From a44f65bade5694c8b25da0e8078ce5eeb101d307 Mon Sep 17 00:00:00 2001 From: Ben Thomas Date: Sun, 11 Jul 2021 13:07:48 +0100 Subject: [PATCH] RL: Perform voxel-by-voxel division (#82) * RL: Perform voxel-by-voxel division * Set zero clipping based on image max * Bump version to 1.2.10 --- CMakeLists.txt | 2 +- lib/petpvcRLPVCImageFilter.h | 5 +- lib/petpvcRLPVCImageFilter.txx | 97 ++++++++++++++++++++++++++++------ 3 files changed, 87 insertions(+), 17 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 60727c6..fd2d134 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/lib/petpvcRLPVCImageFilter.h b/lib/petpvcRLPVCImageFilter.h index 4a36a94..3eb690f 100644 --- a/lib/petpvcRLPVCImageFilter.h +++ b/lib/petpvcRLPVCImageFilter.h @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -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; diff --git a/lib/petpvcRLPVCImageFilter.txx b/lib/petpvcRLPVCImageFilter.txx index d2430c1..0c37660 100644 --- a/lib/petpvcRLPVCImageFilter.txx +++ b/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,30 +98,81 @@ 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 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(); @@ -110,9 +180,9 @@ void RichardsonLucyPVCImageFilter< TInputImage > 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; } @@ -120,9 +190,6 @@ void RichardsonLucyPVCImageFilter< TInputImage > float fCurrentEval = fLog; std::cout << n << "\t" << fCurrentEval << std::endl; n++; - - //if ( fCurrentEval < this->m_fStopCriterion ) - // bStopped = true; } std::cout << std::endl;