ENH: set division threshold for SART

My attempt at fixing #151

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>, OutputImageType)
 
   TRY_AND_EXIT_ON_ITK_EXCEPTION(sart->Update())

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

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<ProjectionType, ProjectionType>;
@@ -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.
+   */
+  void
+  SetDivisionThreshold(ProjectionPixelType threshold);
+
 protected:
   SARTConeBeamReconstructionFilter();
   ~SARTConeBeamReconstructionFilter() override = default;

include/rtkSARTConeBeamReconstructionFilter.hxx

@@ -119,6 +119,13 @@ SARTConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::SetGatingWeigh
   m_IsGated = true;
 }
 
+template <class TVolumeImage, class TProjectionImage>
+void
+SARTConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::SetDivisionThreshold(ProjectionPixelType threshold)
+{
+  m_DivideProjectionFilter->SetThreshold(threshold);
+}
+
 template <class TVolumeImage, class TProjectionImage>
 void
 SARTConeBeamReconstructionFilter<TVolumeImage, TProjectionImage>::VerifyPreconditions() ITKv5_CONST