RL: Perform voxel-by-voxel division

CMakeLists.txt

@@ -28,7 +28,7 @@ endif()
 
 SET(VERSION_MAJOR 1)
 SET(VERSION_MINOR 2)
-SET(VERSION_PATCH 9)
+SET(VERSION_PATCH 10)
 
 MAKE_DIRECTORY(${PROJECT_BINARY_DIR}/config)
 

lib/petpvcRLPVCImageFilter.h

@@ -29,6 +29,7 @@
 #include <itkExtractImageFilter.h>
 #include <itkMultiplyImageFilter.h>
 #include <itkDivideImageFilter.h>
+#include <itkCastImageFilter.h>
 #include <itkAddImageFilter.h>
 #include <itkSubtractImageFilter.h>
 #include <itkDiscreteGaussianImageFilter.h>
@@ -89,7 +90,6 @@ class RichardsonLucyPVCImageFilter:public ImageToImageFilter< TInputImage, TInpu
         this->m_vecVariance = vec;
     }
 
-
     ITKVectorType GetPSF() {
         return this->m_vecVariance;
     }
@@ -114,6 +114,9 @@ class RichardsonLucyPVCImageFilter:public ImageToImageFilter< TInputImage, TInpu
     /** Does the real work. */
     virtual void GenerateData() ITK_OVERRIDE;
 
+    // Calculate threshold at which a value is zeroed
+    float GetZeroThreshold( typename TInputImage::ConstPointer img );
+
     ITKVectorType m_vecVariance;
     unsigned int m_nIterations;
     float m_fStopCriterion;

lib/petpvcRLPVCImageFilter.txx

@@ -42,6 +42,25 @@ RichardsonLucyPVCImageFilter< TInputImage >
     this->m_fStopCriterion = -3e+06;
 }
 
+template< class TInputImage >
+float RichardsonLucyPVCImageFilter< TInputImage >
+::GetZeroThreshold( typename TInputImage::ConstPointer img )
+{
+    // Default threshold to zero values at:
+    const float fZeroThreshold = 1e-4f;
+
+    // Calculate image statistics
+    typename StatisticsFilterType::Pointer statsFilter = StatisticsFilterType::New();
+    statsFilter->SetInput( img );
+    statsFilter->Update();
+
+    // Return default or (image max * default ) as new threshold
+    const float fNewThreshold = std::min( fZeroThreshold, statsFilter->GetMaximum() * fZeroThreshold );
+
+    // Set to zero if somehow new threshold is negative!
+    return std::max( fNewThreshold , 0.0f );
+}
+
 template< class TInputImage >
 void RichardsonLucyPVCImageFilter< TInputImage >
 ::GenerateData()
@@ -79,50 +98,98 @@ void RichardsonLucyPVCImageFilter< TInputImage >
     typename TInputImage::Pointer imageEstimate;
     //Set image estimate to the original non-negative PET data for the first iteration.
     imageEstimate = duplicator->GetOutput();
+    imageEstimate->DisconnectPipeline();
+
+    // Create a blank image with same dimensions as input to use for division
+    duplicator->SetInputImage(imageEstimate);
+    duplicator->Update();
+
+    typedef itk::CastImageFilter< TInputImage, TInputImage > CastImageFilterType;
+    typename CastImageFilterType::Pointer castFilter = CastImageFilterType::New();
+    castFilter->SetInput( duplicator->GetOutput() );
+    castFilter->Update();
+    typename TInputImage::Pointer blankImage = castFilter->GetOutput();
+
+    typedef typename TInputImage::PixelType PixelType;
+    blankImage->FillBuffer( itk::NumericTraits< PixelType >::Zero );
+
+    // Image to store result of division
+    typename TInputImage::Pointer dividedImage;
 
     int nMaxNumOfIters =  this->m_nIterations;
     int n=1;
 
     bool bStopped = false;
-
 
     while ( ( n <= nMaxNumOfIters ) && ( !bStopped ) ) {
-         float fLog = 0.0;   
+         float fLog = 0.0;  
+            // f_k * h
             blurFilter->SetInput( imageEstimate );
-			divideFilter->SetInput1( thresholdFilter->GetOutput() );
-			divideFilter->SetInput2( blurFilter->GetOutput() );
+            blurFilter->Update();
+
+            // Perform f(x) / [ f_k(x) * h ] voxel-by-voxel as itk::DivideFilterType sets image
+            // to max if denominator is 0.
+            ConstImageIterator numeratorIt( thresholdFilter->GetOutput(),  thresholdFilter->GetOutput()->GetLargestPossibleRegion() );
+            ConstImageIterator denomIt( blurFilter->GetOutput(),  blurFilter->GetOutput()->GetLargestPossibleRegion() );
+
+            duplicator->SetInputImage(blankImage);
+            duplicator->Update();
+            dividedImage = duplicator->GetOutput();
+            dividedImage->DisconnectPipeline();
+
+            ImageRegionIterator<TInputImage> divIt( dividedImage, dividedImage->GetLargestPossibleRegion() );
+
+            numeratorIt.GoToBegin();
+            denomIt.GoToBegin();
+            divIt.GoToBegin();
+
+            // Get zero threshold
+            const float fSmallNum = this->GetZeroThreshold( blurFilter->GetOutput() );
+
+            while ( !divIt.IsAtEnd() )
+    	    {
+				if ( denomIt.Get() > fSmallNum )
+                    divIt.Set(numeratorIt.Get() / denomIt.Get());
+                else
+                    divIt.Set(itk::NumericTraits< PixelType >::Zero);
+
+        		++numeratorIt;
+                ++denomIt;
+                ++divIt;
+        	}
+
+            // Reblur correction factors            
+            blurFilter2->SetInput( dividedImage );
+            blurFilter2->Update();
+
+            // Multiply current image estimate by reblurred correction factors
 			multiplyFilter->SetInput1( imageEstimate );
-			multiplyFilter->SetInput2( divideFilter->GetOutput() );
+			multiplyFilter->SetInput2( blurFilter2->GetOutput() );
             multiplyFilter->Update();
-
+
+            // Update image estimate 
             imageEstimate = multiplyFilter->GetOutput();
             imageEstimate->DisconnectPipeline();
-
-			blurFilter2->SetInput( imageEstimate );
-			blurFilter2->Update();
 
             ConstImageIterator origIt( thresholdFilter->GetOutput(),  thresholdFilter->GetOutput()->GetLargestPossibleRegion() );
-            ConstImageIterator currIt(  blurFilter2->GetOutput(), blurFilter2->GetOutput()->GetLargestPossibleRegion() );
+            ConstImageIterator currIt( imageEstimate, imageEstimate->GetLargestPossibleRegion() );
 
             origIt.GoToBegin();
             currIt.GoToBegin();
 
             while ( !currIt.IsAtEnd() )
     	    {
 				if (currIt.Get() > 0.0) {
-                	float diff = origIt.Get() * log(currIt.Get()) - currIt.Get();
-	                fLog += diff;
+                    fLog += origIt.Get() * log(currIt.Get()) - currIt.Get();
 				}
+
         		++currIt;
                 ++origIt;
         	}
 
             float fCurrentEval = fLog;
             std::cout << n << "\t" << fCurrentEval << std::endl;         
 			n++;
-
-            //if ( fCurrentEval < this->m_fStopCriterion )
-            //   bStopped = true;
     }
     std::cout << std::endl;