From a0e604c77d727c05795e13f7623b53e9b3b2fc1b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Tue, 14 Apr 2020 11:06:31 +0200 Subject: [PATCH] ENH: set division threshold for SART My attempt at fixing #151 --- applications/rtksart/rtksart.cxx | 5 +++++ applications/rtksart/rtksart.ggo | 1 + include/rtkSARTConeBeamReconstructionFilter.h | 10 ++++++++++ include/rtkSARTConeBeamReconstructionFilter.hxx | 1 + 4 files changed, 17 insertions(+) diff --git a/applications/rtksart/rtksart.cxx b/applications/rtksart/rtksart.cxx index c73d30e4c..f2477237d 100644 --- a/applications/rtksart/rtksart.cxx +++ b/applications/rtksart/rtksart.cxx @@ -119,6 +119,11 @@ main(int argc, char * argv[]) sart->SetEnforcePositivity(true); } + if (args_info.divisionthreshold_given) + { + sart->SetDivisionThreshold(args_info.divisionthreshold_arg); + } + REPORT_ITERATIONS(sart, rtk::SARTConeBeamReconstructionFilter, OutputImageType) TRY_AND_EXIT_ON_ITK_EXCEPTION(sart->Update()) diff --git a/applications/rtksart/rtksart.ggo b/applications/rtksart/rtksart.ggo index d3456d89d..80867d10f 100644 --- a/applications/rtksart/rtksart.ggo +++ b/applications/rtksart/rtksart.ggo @@ -11,6 +11,7 @@ option "positivity" - "Enforces positivity during the reconstruction" f option "input" i "Input volume" string no option "nprojpersubset" - "Number of projections processed between each update of the reconstructed volume (1 for SART, several for OSSART, all for SIRT)" int no default="1" option "nodisplaced" - "Disable the displaced detector filter" flag off +option "divisionthreshold" - "Threshold below which pixels in the denominator in the projection space are considered zero" double no section "Phase gating" option "signal" - "File containing the phase of each projection" string no diff --git a/include/rtkSARTConeBeamReconstructionFilter.h b/include/rtkSARTConeBeamReconstructionFilter.h index c3dd1d981..2b782239b 100644 --- a/include/rtkSARTConeBeamReconstructionFilter.h +++ b/include/rtkSARTConeBeamReconstructionFilter.h @@ -147,6 +147,7 @@ class ITK_EXPORT SARTConeBeamReconstructionFilter /** Some convenient type alias. */ using VolumeType = TVolumeImage; using ProjectionType = TProjectionImage; + using ProjectionPixelType = typename ProjectionType::PixelType; /** Typedefs of each subfilter of this composite filter */ using ExtractFilterType = itk::ExtractImageFilter; @@ -209,6 +210,13 @@ class ITK_EXPORT SARTConeBeamReconstructionFilter itkSetMacro(DisableDisplacedDetectorFilter, bool); itkGetMacro(DisableDisplacedDetectorFilter, bool); + /** Set the threshold below which pixels in the denominator in the projection space are considered zero. The division + * by zero will then be evaluated at zero. Avoid noise magnification from low projections values when working with + * noisy and/or simulated data. + */ + itkSetMacro(DivisionThreshold, ProjectionPixelType); + itkGetMacro(DivisionThreshold, ProjectionPixelType); + protected: SARTConeBeamReconstructionFilter(); ~SARTConeBeamReconstructionFilter() override = default; @@ -251,6 +259,8 @@ class ITK_EXPORT SARTConeBeamReconstructionFilter typename DisplacedDetectorFilterType::Pointer m_DisplacedDetectorFilter; typename GatingWeightsFilterType::Pointer m_GatingWeightsFilter; + ProjectionPixelType m_DivisionThreshold; + bool m_EnforcePositivity; bool m_DisableDisplacedDetectorFilter; diff --git a/include/rtkSARTConeBeamReconstructionFilter.hxx b/include/rtkSARTConeBeamReconstructionFilter.hxx index 0c1c7b2d9..a1a0e5c99 100644 --- a/include/rtkSARTConeBeamReconstructionFilter.hxx +++ b/include/rtkSARTConeBeamReconstructionFilter.hxx @@ -155,6 +155,7 @@ void SARTConeBeamReconstructionFilter::GenerateOutputInformation() { m_DisplacedDetectorFilter->SetDisable(m_DisableDisplacedDetectorFilter); + m_DivideProjectionFilter->SetThreshold(m_DivisionThreshold); // We only set the first sub-stack at that point, the rest will be // requested in the GenerateData function