Added time driven Monitor to CM Solver

modules/cudaconfidencemaps/core/cudaconfidencemaps_cuda.cu

 #include "cudaconfidencemaps_cuda.h"
 #include "cudautils.h"
+#include "cuspmonitors.h"
 
 #include <cusp/blas.h>
 #include <cusp/dia_matrix.h>
@@ -69,32 +70,6 @@ namespace cuda {
         float systemSolveTime;
     };
 
-    template <typename ValueType>
-    class iteration_monitor : public cusp::default_monitor<ValueType>
-    {
-        typedef typename cusp::norm_type<ValueType>::type Real;
-        typedef cusp::default_monitor<ValueType> super;
-
-    public:
-        template <typename Vector>
-        iteration_monitor(const Vector& b, size_t iteration_limit = 500)
-            : super(b, iteration_limit, 0.0f, 0.0f)
-        { }
-
-        template <typename Vector>
-        bool finished(const Vector& r)
-        {
-            // Only if the maximum iteration count has been reached, actually go ahead and
-            // compute the error
-            if (super::iteration_count() >= super::iteration_limit()) {
-                super::r_norm = cusp::blas::nrm2(r);
-                return true;
-            }
-
-            return false;
-        }
-    };
-
     CudaConfidenceMapsSystemSolver::CudaConfidenceMapsSystemSolver()
         : _gpuData(new CudaConfidenceMapsSystemGPUData())
     {
@@ -173,12 +148,12 @@ namespace cuda {
     }
 
     // FIXME: Remove errorTolerance parameter
-    void CudaConfidenceMapsSystemSolver::solve(int maximumIterations, float errorTolerance) {
+    void CudaConfidenceMapsSystemSolver::solve(float millisecondBudget) {
         // Measure execution time and record it in the _gpuData datastructure
         CUDAClock clock; clock.start();
 
         // The solution is computed using Conjugate Gradient with a Diagonal (Jacobi) preconditioner
-        iteration_monitor<float> monitor(_gpuData->b_d, maximumIterations);
+        deadline_monitor<float> monitor(_gpuData->b_d, millisecondBudget);
         cusp::precond::diagonal<float, cusp::device_memory> M(_gpuData->L_d);
         cusp::krylov::cg(_gpuData->L_d, _gpuData->x_d, _gpuData->b_d, monitor, M);
         _gpuData->solutionResidualNorm = monitor.residual_norm();

modules/cudaconfidencemaps/core/cudaconfidencemaps_cuda.h

@@ -58,10 +58,9 @@ namespace cuda {
         /**
          * After calling \see uploadImage(), this functions launches a solver on the GPU that will solve
          * the diffusion problem.
-         * \param   maximumIterations maximum number of iterations the solver will preform
-         * \param   errorTolerance    if the solution error sinks below this value, the solver stops early
+         * \param   millisecondBudget the time budget the solver has to come up with a solution.
          */
-        void solve(int maximumIterations, float errorTolerance);
+        void solve(float millisecondBudget);
 
         /**
          * Returns a host buffer of the last solution computed by the solver. The pointer is guaranteed to

modules/cudaconfidencemaps/core/cuspmonitors.h

+#ifndef CUSPMONITORS_H__
+#define CUSPMONITORS_H__
+
+#include <cusp/monitor.h>
+#include <tbb/tick_count.h>
+
+namespace campvis {
+namespace cuda {
+
+    /**
+     * This class allows to execute a fixed number of conjugate gradient iterations
+     * in CUSP. Unlike the default_monitor, this class only calculates the residual
+     * norm when the iteration count is reached.
+     */
+    template <typename ValueType>
+    class iteration_monitor : public cusp::default_monitor<ValueType>
+    {
+        typedef typename cusp::norm_type<ValueType>::type Real;
+        typedef cusp::default_monitor<ValueType> super;
+
+    public:
+        template <typename Vector>
+        iteration_monitor(const Vector& b, size_t iteration_limit = 500)
+            : super(b, iteration_limit, 0.0f, 0.0f)
+        { }
+
+        template <typename Vector>
+        bool finished(const Vector& r)
+        {
+            // Only if the maximum iteration count has been reached, actually go ahead and
+            // compute the error
+            if (super::iteration_count() >= super::iteration_limit()) {
+                super::r_norm = cusp::blas::nrm2(r);
+                return true;
+            }
+
+            return false;
+        }
+    };
+
+    /**
+     * This monitor allows to set a deadline, after which the computation has to stop.
+     */
+    template <typename ValueType>
+    class deadline_monitor : public cusp::default_monitor<ValueType>
+    {
+        typedef typename cusp::norm_type<ValueType>::type Real;
+        typedef cusp::default_monitor<ValueType> super;
+
+    public:
+        template <typename Vector>
+        deadline_monitor(const Vector& b, float milliseconds)
+            : super(b, 0, 0.0f, 0.0f), _startTime(tbb::tick_count::now()), _seconds(milliseconds / 1000.0f)
+        { }
+
+        template <typename Vector>
+        bool finished(const Vector& r)
+        {
+            // Only if the deadline is reached, stop and compute the error
+            if ((tbb::tick_count::now() - _startTime).seconds() > _seconds) {
+                super::r_norm = cusp::blas::nrm2(r);
+                return true;
+            }
+
+            return false;
+        }
+
+    private:
+        tbb::tick_count _startTime;
+        float _seconds;
+    };
+}
+}
+
+#endif // CUSPMONITORS_H__
\ No newline at end of file

modules/cudaconfidencemaps/cudaconfidencemaps.cmake

@@ -47,7 +47,10 @@ IF(ModuleEnabled)
         ENDIF()
 
         # CUSP Include directory
-        CUDA_INCLUDE_DIRECTORIES(${ThisModDir}/ext/cusplibrary-0.4.0)
+        CUDA_INCLUDE_DIRECTORIES(
+            ${ThisModDir}/ext/cusplibrary-0.4.0
+            ${TBB_INCLUDE_DIR}
+            )
 
         # Build CUDA portion of the code (STATICALLY!?)
         FILE(GLOB cuda_SOURCES modules/cudaconfidencemaps/core/*.cu)
@@ -55,10 +58,6 @@ IF(ModuleEnabled)
             ${cuda_SOURCES}
             )        
 
-        # Make sure code can find the CUSP include files included with this module
-        #SET(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS};"-I ${ThisModDir}/ext/cusplibrary-0.4.0")
-        set(CUDA_NVCC_FLAGS "ajkladjfl" CACHE STRING "adsf")
-
         # Link CUDA code to module
         LIST(APPEND ThisModExternalLibs cudaconfidencemaps-cuda)
     ELSE()

modules/cudaconfidencemaps/pipelines/cudaconfidencemapsdemo.cpp

@@ -30,18 +30,17 @@
 
 namespace campvis {
 
-	CudaConfidenceMapsDemo::CudaConfidenceMapsDemo(DataContainer* dc)
-		: AutoEvaluationPipeline(dc)
-		, _usIgtlReader()
-		, _usCropFilter(&_canvasSize)
+    CudaConfidenceMapsDemo::CudaConfidenceMapsDemo(DataContainer* dc)
+        : AutoEvaluationPipeline(dc)
+        , _usIgtlReader()
+        , _usCropFilter(&_canvasSize)
         , _usBlurFilter(&_canvasSize)
         , _usResampler(&_canvasSize)
         , _usMapsSolver()
         , _usFusion(&_canvasSize)
         , _usFanRenderer(&_canvasSize)
-        , p_iterations("Iterations", "Number of CG Iterations", 150, 1, 500)
         , p_autoIterationCount("AutoIterationCount", "Estimate iteration count based on time slot", true)
-        , p_timeSlot("TimeSlot", "Milliseconds per frame", 32.0f, 10.0f, 250.0f)
+        , p_millisecondBudget("MillisecondBudget", "Milliseconds per frame", 32.0f, 10.0f, 1000.0f)
         , p_connectToIGTLinkServer("ConnectToIGTLink", "Connect/Disconnect to IGTLink")
         , p_gaussianFilterSize("GaussianSigma", "Blur amount", 2.5f, 1.0f, 10.0f)
         , p_resamplingScale("ResampleScale", "Resample Scale", 0.25f, 0.01f, 1.0f)
@@ -52,21 +51,18 @@ namespace campvis {
         , p_useAlphaBetaFilter("UseAlphaBetaFilter", "Alpha-Beta-Filter", true)
         , p_fanHalfAngle("FanHalfAngle", "Fan Half Angle", 37.0f, 1.0f, 90.0f)
         , p_fanInnerRadius("FanInnerRadius", "Fan Inner Radius", 0.222f, 0.001f, 0.999f)
-        , _cgIterationsPerMsRunningAverage(1.0f)
-        , _cgTimeslotRunningAverage(1.0f)
         , _statisticsLastUpdateTime()
     {
-		addProcessor(&_usIgtlReader);
-		addProcessor(&_usCropFilter);
+        addProcessor(&_usIgtlReader);
+        addProcessor(&_usCropFilter);
         addProcessor(&_usBlurFilter);
         addProcessor(&_usResampler);
         addProcessor(&_usMapsSolver);
         addProcessor(&_usFusion);
         addProcessor(&_usFanRenderer);
 
-        addProperty(p_iterations);
         addProperty(p_autoIterationCount);
-        addProperty(p_timeSlot);
+        addProperty(p_millisecondBudget);
         addProperty(p_connectToIGTLinkServer);
         addProperty(p_gaussianFilterSize);
         addProperty(p_resamplingScale);
@@ -98,8 +94,8 @@ namespace campvis {
         // Create connectors
         _usIgtlReader.p_targetImagePrefix.setValue("us.igtl.");
 
-		_usCropFilter.p_inputImage.setValue("us.igtl.CAMPUS");
-		_usCropFilter.p_outputImage.setValue("us");
+        _usCropFilter.p_inputImage.setValue("us.igtl.CAMPUS");
+        _usCropFilter.p_outputImage.setValue("us");
 
         _usBlurFilter.p_inputImage.setValue("us");
         _usBlurFilter.p_outputImage.setValue("us.blurred");
@@ -125,9 +121,6 @@ namespace campvis {
         _renderTargetID.setValue("us.fused_fan");
 
         // Bind pipeline proeprties to processor properties
-        p_iterations.addSharedProperty(&_usMapsSolver.p_iterations);
-        //p_autoIterationCount.addSharedProperty();
-        //p_timeSlot.addSharedProperty();
         p_connectToIGTLinkServer.addSharedProperty(&_usIgtlReader.p_connect);
         p_gaussianFilterSize.addSharedProperty(&_usBlurFilter.p_sigma);
         p_resamplingScale.addSharedProperty(&_usResampler.p_resampleScale);
@@ -151,6 +144,7 @@ namespace campvis {
         else {
             // Only launch the pipeline if the IgtlReader has recieved new data
             if (!_usIgtlReader.isValid()) {
+                float millisecondBudget = p_millisecondBudget.getValue();
                 auto startTime = tbb::tick_count::now();
 
                 // Make sure that the whole pipeline gets invalidated
@@ -160,12 +154,13 @@ namespace campvis {
                 _usMapsSolver.invalidate(AbstractProcessor::INVALID_RESULT);
                 _usFusion.invalidate(AbstractProcessor::INVALID_RESULT);
 
-				executeProcessorAndCheckOpenGLState(&_usIgtlReader);
-				executeProcessorAndCheckOpenGLState(&_usCropFilter);
-                executeProcessorAndCheckOpenGLState(&_usBlurFilter);
+                executeProcessorAndCheckOpenGLState(&_usIgtlReader);
+                executeProcessorAndCheckOpenGLState(&_usCropFilter);
+                executeProcessorAndCheckOpenGLState(&_usBlurFilter) ;
                 executeProcessorAndCheckOpenGLState(&_usResampler);
 
                 auto solverStartTime = tbb::tick_count::now();
+                _usMapsSolver.p_millisecondBudget.setValue(millisecondBudget);
                 executeProcessorAndCheckOpenGLState(&_usMapsSolver);
                 auto solverEndTime = tbb::tick_count::now();
 
@@ -186,24 +181,6 @@ namespace campvis {
                     string << "Error: " << _usMapsSolver.getResidualNorm() << std::endl;
                     _usFanRenderer.p_text.setValue(string.str());
                 }
-
-
-                auto ms = (solverEndTime - solverStartTime).seconds() * 1000.0f;
-                auto iterationsPerMs = _usMapsSolver.getActualConjugentGradientIterations() / ms;
-
-
-                // Factor out the time needed for pre- and postprocessing the image from the time slot
-                float timeSlot = (p_timeSlot.getValue() -
-                                 (solverStartTime - startTime).seconds() * 1000.0f -
-                                 (endTime - solverEndTime).seconds() * 1000.0f);
-
-                // Compute the running average using an exponential filter
-                float expAlpha = 0.05f;
-                _cgIterationsPerMsRunningAverage = _cgIterationsPerMsRunningAverage * (1.0f - expAlpha) + iterationsPerMs * expAlpha;
-                _cgTimeslotRunningAverage = _cgTimeslotRunningAverage * (1.0f - expAlpha) + timeSlot * expAlpha;
-
-                int iterations = static_cast<int>(_cgTimeslotRunningAverage * _cgIterationsPerMsRunningAverage);
-				p_iterations.setValue(iterations);
             }
         }
     }

modules/cudaconfidencemaps/pipelines/cudaconfidencemapsdemo.h

@@ -72,7 +72,7 @@ namespace campvis {
 
     protected:
         OpenIGTLinkClient _usIgtlReader;
-		GlImageCrop      _usCropFilter;
+        GlImageCrop      _usCropFilter;
         GlGaussianFilter _usBlurFilter;
         GlImageResampler _usResampler;
         CudaConfidenceMapsSolver _usMapsSolver;
@@ -80,9 +80,8 @@ namespace campvis {
         UsFanRenderer _usFanRenderer;
 
 
-        NumericProperty<int> p_iterations;
         BoolProperty   p_autoIterationCount;
-        FloatProperty  p_timeSlot;
+        FloatProperty  p_millisecondBudget;
 
         ButtonProperty p_connectToIGTLinkServer;
 
@@ -98,9 +97,6 @@ namespace campvis {
         FloatProperty  p_fanHalfAngle;
         FloatProperty  p_fanInnerRadius;
 
-
-        float _cgIterationsPerMsRunningAverage;
-        float _cgTimeslotRunningAverage;
         tbb::tick_count _statisticsLastUpdateTime;
     };
 

modules/cudaconfidencemaps/processors/cudaconfidencemapssolver.cpp

@@ -41,7 +41,7 @@ namespace campvis {
         , p_outputConfidenceMap("OutputConfidenceMap", "Output Confidence Map", "us.confidence", DataNameProperty::WRITE)
         , p_resetResult("ResetSolution", "Reset solution vector")
         , p_use8Neighbourhood("Use8Neighbourhood", "Use 8 Neighbourhood (otherwise 4)", true)
-        , p_iterations("IterationCount", "Conjugate Gradient Iterations", 200, 1, 500)
+        , p_millisecondBudget("IterationCount", "Conjugate Gradient Iterations", 25.0f, 1.0f, 1000.0f)
         , p_gradientScaling("GradientScaling", "Scaling factor for gradients", 2.0f, 0.001f, 10.0f)
         , p_paramAlpha("Alpha", "Alpha (TGC)", 2.0f, 0.001f, 10.0f)
         , p_paramBeta("Beta", "Beta (Weight mapping)", 20.0f, 0.001f, 200.0f)
@@ -57,7 +57,7 @@ namespace campvis {
 
         addProperty(p_resetResult);
         addProperty(p_use8Neighbourhood);
-        addProperty(p_iterations);
+        addProperty(p_millisecondBudget);
         addProperty(p_gradientScaling);
         addProperty(p_paramAlpha);
         addProperty(p_paramBeta);
@@ -84,7 +84,7 @@ namespace campvis {
         ImageRepresentationLocal::ScopedRepresentation img(data, p_inputImage.getValue());
         if (img != 0) {
             bool use8Neighbourhood = p_use8Neighbourhood.getValue();
-            int iterations = p_iterations.getValue();
+            float millisecondBudget = p_millisecondBudget.getValue();
             float gradientScaling = p_gradientScaling.getValue();
             float alpha = p_paramAlpha.getValue();
             float beta = p_paramBeta.getValue();
@@ -99,7 +99,7 @@ namespace campvis {
             auto image = (unsigned char*)img->getWeaklyTypedPointer()._pointer;
 
             _solver.uploadImage(image, size.x, size.y, gradientScaling, alpha, beta, gamma, use8Neighbourhood);
-            _solver.solve(iterations, 1e-10f);
+            _solver.solve(millisecondBudget);
 
             const float *solution = _solver.getSolution(size.x, size.y);
 

modules/cudaconfidencemaps/processors/cudaconfidencemapssolver.h

@@ -90,7 +90,7 @@ namespace campvis {
 
         BoolProperty  p_use8Neighbourhood;          ///< Wether to use 8- or 4-neighbourhood
 
-        NumericProperty<int> p_iterations;          ///< Number of CG-Iterations to do
+        FloatProperty p_millisecondBudget;          ///< Maximum number of ms the solver can run
 
         FloatProperty p_gradientScaling;
         FloatProperty p_paramAlpha;

modules/cudaconfidencemaps/processors/usfanrenderer.cpp

@@ -71,7 +71,7 @@ namespace campvis {
         // Creates the grid, with the origin at the center of the top edge, with the +y axis representing depth
         _grid = GeometryDataFactory::createGrid(cgt::vec3(-0.5f, 1.0f, 0.0f), cgt::vec3(0.5f, 0.0f, 0.0f),
                                                 cgt::vec3(0.0f, 1.0f, 0.0f), cgt::vec3(1.0f, 0.0f, 0.0f),
-                                                16, 4);
+                                                32, 32);
 
         // Initialize font rendering
         updateFontAtlas();

modules/openigtlink/processors/openigtlinkclient.cpp

@@ -160,7 +160,7 @@ namespace campvis {
                 // If the voxel size boundled with the packet is practically 0.0f, make it 1.0f
                 // this makes sure we don't get invalid mapping informations (non invertable matrix)
                 if (minElem(voxelSize) <= 1e-10f) {
-                    voxelSize = 1.0f;
+                    voxelSize = cgt::vec3(1.0f);
                 }
 
                 size_t dimensionality = (size_i[2] == 1) ? ((size_i[1] == 1) ? 1 : 2) : 3;

modules/preprocessing/glsl/glimageresampler.frag

@@ -35,7 +35,7 @@ uniform sampler2D _texture;
 #endif
 
 void main() {
-	#ifdef GLRESAMPLER_3D
+    #ifdef GLRESAMPLER_3D
     vec4 sample = texture(_texture, vec3(ex_TexCoord.xy, _zTexCoord));
     #endif