ConfidenceMapsSolver now supports fixed iterations

added support for fixed iteration count to ConfidenceMapsSolver again.
It is now possible to chose wether to use a time budget or a fixed
iteration count.

modules/cudaconfidencemaps/core/cudaconfidencemaps_cuda.cu

@@ -153,7 +153,7 @@ namespace cuda {
     }
 
     // FIXME: Remove errorTolerance parameter
-    void CudaConfidenceMapsSystemSolver::solve(float millisecondBudget) {
+    void CudaConfidenceMapsSystemSolver::solveWithFixedTimeBudget(float millisecondBudget) {
         // Measure execution time and record it in the _gpuData datastructure
         CUDAClock clock; clock.start();
 
@@ -164,24 +164,25 @@ namespace cuda {
         _gpuData->solutionResidualNorm = monitor.residual_norm();
         _gpuData->iterationCount = static_cast<int>(monitor.iteration_count());
 
-        if (alphaBetaFilterEnabled()) {
-            // X' = X' + V'
-            // R' = X - X'
-            cusp::blas::axpy(_gpuData->abFilterV_d, _gpuData->abFilterX_d, 1.0f);
-            cusp::blas::axpby(_gpuData->x_d, _gpuData->abFilterX_d, _gpuData->abFilterR_d, 1.0f, -1.0f);
+        // Applies alpha-beta-filtering (if enabled) and downloads solution
+        performOutputFiltering();
 
-            // X' = X' + alpha * R'
-            // V' = V' + beta * R'
-            cusp::blas::axpy(_gpuData->abFilterR_d, _gpuData->abFilterX_d, _gpuData->abFilterAlpha);
-            cusp::blas::axpy(_gpuData->abFilterR_d, _gpuData->abFilterV_d, _gpuData->abFilterBeta);
+        _gpuData->systemSolveTime = clock.getElapsedMilliseconds();
+    }
 
-            // Download the smoothed solution to the host
-            _gpuData->x_h = _gpuData->abFilterX_d;
-        }
-        else {
-            // Downlaod the actual solution, which has been computed to the host
-            _gpuData->x_h = _gpuData->x_d;
-        }
+    void CudaConfidenceMapsSystemSolver::solveWithFixedIterationCount(int iterations) {
+        // 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, iterations);
+        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();
+        _gpuData->iterationCount = static_cast<int>(monitor.iteration_count());
+
+        // Applies alpha-beta-filtering (if enabled) and downloads solution
+        performOutputFiltering();
 
         _gpuData->systemSolveTime = clock.getElapsedMilliseconds();
     }
@@ -209,6 +210,26 @@ namespace cuda {
         return _gpuData->systemSolveTime;
     }
 
+    void CudaConfidenceMapsSystemSolver::performOutputFiltering() {
+        if (alphaBetaFilterEnabled()) {
+            // X' = X' + V'
+            // R' = X - X'
+            cusp::blas::axpy(_gpuData->abFilterV_d, _gpuData->abFilterX_d, 1.0f);
+            cusp::blas::axpby(_gpuData->x_d, _gpuData->abFilterX_d, _gpuData->abFilterR_d, 1.0f, -1.0f);
+
+            // X' = X' + alpha * R'
+            // V' = V' + beta * R'
+            cusp::blas::axpy(_gpuData->abFilterR_d, _gpuData->abFilterX_d, _gpuData->abFilterAlpha);
+            cusp::blas::axpy(_gpuData->abFilterR_d, _gpuData->abFilterV_d, _gpuData->abFilterBeta);
+
+            // Download the smoothed solution to the host
+            _gpuData->x_h = _gpuData->abFilterX_d;
+        }
+        else {
+            // Downlaod the actual solution, which has been computed to the host
+            _gpuData->x_h = _gpuData->x_d;
+        }
+    }
 
     void CudaConfidenceMapsSystemSolver::resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown, bool use8Neighbourhood) {
         // If the problem size changed, reset the solution vector, as well as all

modules/cudaconfidencemaps/core/cudaconfidencemaps_cuda.h

@@ -56,11 +56,18 @@ namespace cuda {
         void setAlphaBetaFilterParameters(float alpha, float beta);
 
         /**
-         * After calling \see uploadImage(), this functions launches a solver on the GPU that will solve
+         * After calling \see uploadImage(), this function launches a solver on the GPU that will solve
          * the diffusion problem.
          * \param   millisecondBudget the time budget the solver has to come up with a solution.
          */
-        void solve(float millisecondBudget);
+        void solveWithFixedTimeBudget(float millisecondBudget);
+
+        /**
+         * After calling \see uploadImage(), this function lanuches a solver on the GPU that will solve
+         * the diffusion problem.
+         * \param  iterations the number of conjugate iterations the solver is allowed to take.
+         */ 
+        void solveWithFixedIterationCount(int iterations);
 
         /**
          * Returns a host buffer of the last solution computed by the solver. The pointer is guaranteed to
@@ -95,6 +102,7 @@ namespace cuda {
         float getSystemSolveTime() const;
 
     private:
+        void performOutputFiltering();
         void resizeDataStructures(int imageWidth, int imageHeight, bool isUpsideDown, bool use8Neighbourhood);
         void createSystemGPU(const unsigned char* imageData, int imageWidth, int imageHeight,
                              float gradientScaling, float alpha, float beta, float gamma, bool isUpsideDown=true);

modules/cudaconfidencemaps/processors/cudaconfidencemapssolver.cpp

@@ -41,7 +41,9 @@ 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_millisecondBudget("IterationCount", "Conjugate Gradient Iterations", 25.0f, 1.0f, 1000.0f)
+        , p_useFixedIterationCount("UseFixedIterationCount", "Use Fixed Iteration Count", false)
+        , p_millisecondBudget("MillisecondBudget", "(P)CG Solver Time Budget", 25.0f, 1.0f, 1000.0f)
+        , p_iterationBudget("IterationBudget", "(P)CG Solver Iteration Count", 100, 0, 1000)
         , 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 +59,9 @@ namespace campvis {
 
         addProperty(p_resetResult);
         addProperty(p_use8Neighbourhood);
+        addProperty(p_useFixedIterationCount);
         addProperty(p_millisecondBudget);
+        addProperty(p_iterationBudget);
         addProperty(p_gradientScaling);
         addProperty(p_paramAlpha);
         addProperty(p_paramBeta);
@@ -66,6 +70,8 @@ namespace campvis {
         addProperty(p_useAlphaBetaFilter);
         addProperty(p_filterAlpha);
         addProperty(p_filterBeta);
+
+        updatePropertyVisibility();
     }
 
     CudaConfidenceMapsSolver::~CudaConfidenceMapsSolver() { }
@@ -84,7 +90,6 @@ namespace campvis {
         ImageRepresentationLocal::ScopedRepresentation img(data, p_inputImage.getValue());
         if (img != 0) {
             bool use8Neighbourhood = p_use8Neighbourhood.getValue();
-            float millisecondBudget = p_millisecondBudget.getValue();
             float gradientScaling = p_gradientScaling.getValue();
             float alpha = p_paramAlpha.getValue();
             float beta = p_paramBeta.getValue();
@@ -98,8 +103,16 @@ namespace campvis {
             size_t elementCount = cgt::hmul(size);
             auto image = (unsigned char*)img->getWeaklyTypedPointer()._pointer;
 
+            // Copy the image on the GPU and generate the equation system
             _solver.uploadImage(image, size.x, size.y, gradientScaling, alpha, beta, gamma, use8Neighbourhood);
-            _solver.solve(millisecondBudget);
+
+            // Solve the equation system using Conjugate Gradient
+            if (p_useFixedIterationCount.getValue()) {
+                _solver.solveWithFixedIterationCount(p_iterationBudget.getValue());
+            }
+            else {
+                _solver.solveWithFixedTimeBudget(p_millisecondBudget.getValue());
+            }
 
             const float *solution = _solver.getSolution(size.x, size.y);
 
@@ -115,15 +128,15 @@ namespace campvis {
         }
     }
 
-    void CudaConfidenceMapsSolver::updateProperties(DataContainer& dataContainer) { }
+    void CudaConfidenceMapsSolver::onPropertyChanged(const AbstractProperty* prop) {
+        updatePropertyVisibility();
+    }
 
-    int CudaConfidenceMapsSolver::getActualConjugentGradientIterations() const
-    {
+    int CudaConfidenceMapsSolver::getActualConjugentGradientIterations() const {
         return _solver.getSolutionIterationCount();
     }
 
-    float CudaConfidenceMapsSolver::getResidualNorm() const
-    {
+    float CudaConfidenceMapsSolver::getResidualNorm() const {
         return _solver.getSolutionResidualNorm();
     }
 
@@ -131,4 +144,11 @@ namespace campvis {
         // Create a linear gradient image of the same size as the input image
         _solver.resetSolution();
     }
+
+    void CudaConfidenceMapsSolver::updatePropertyVisibility() {
+        // Hide properties that currently do not affect the processor
+        bool useFixedIterationCount = p_useFixedIterationCount.getValue();
+        p_millisecondBudget.setVisible(!useFixedIterationCount);
+        p_iterationBudget.setVisible(useFixedIterationCount);
+    }
 }
\ No newline at end of file

modules/cudaconfidencemaps/processors/cudaconfidencemapssolver.h

@@ -77,8 +77,7 @@ namespace campvis {
          */
         virtual void updateResult(DataContainer& dataContainer);
 
-        /// \see AbstractProcessor::updateProperties
-        virtual void updateProperties(DataContainer& dataContainer);
+        virtual void onPropertyChanged(const AbstractProperty* prop);
 
         int getActualConjugentGradientIterations() const;
         float getResidualNorm() const;
@@ -86,13 +85,15 @@ namespace campvis {
         DataNameProperty p_inputImage;              ///< ID for input volume
         DataNameProperty p_outputConfidenceMap;     ///< ID for output gradient volume
 
-        ButtonProperty p_resetResult;
+        ButtonProperty p_resetResult;               ///< Resets solution vector to a linear confidence ramp
 
         BoolProperty  p_use8Neighbourhood;          ///< Wether to use 8- or 4-neighbourhood
+        BoolProperty  p_useFixedIterationCount;     ///< If set to true the CG iterations are fixed, and not determined by time budget
 
         FloatProperty p_millisecondBudget;          ///< Maximum number of ms the solver can run
+        IntProperty   p_iterationBudget;            ///< Number of CG iterations to perform per frame
 
-        FloatProperty p_gradientScaling;
+        FloatProperty p_gradientScaling;            ///< Multiplicaton factor for gradients
         FloatProperty p_paramAlpha;
         FloatProperty p_paramBeta;
         FloatProperty p_paramGamma;
@@ -103,6 +104,7 @@ namespace campvis {
         
     protected:
         void resetSolutionVector();
+        void updatePropertyVisibility();
 
         cuda::CudaConfidenceMapsSystemSolver _solver;