ENH: Add new inputs for the iterative reconstruction filters

Move the inputs for the clip images and attenuation map from OSEM filter to the base
class of the iterative reconstruction filters so all the filters can uses them.

applications/rtkforwardprojections/rtkforwardprojections.ggo

@@ -13,6 +13,6 @@ section "Projectors"
 option "fp"    f "Forward projection method" values="Joseph","JosephAttenuated","CudaRayCast","Zeng" enum no default="Joseph"
 option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction"   string  no
 option "sigmazero" - "PSF value at a distance of 0 meter of the detector"   double  no
-option "alphapsf" - "Slope of the PSF against the detector distance"   double  no 
+option "alphapsf" - "Slope of the PSF against the detector distance"   double  no
 option "inferiorclipimage" - "Value of the inferior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no
 option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no

applications/rtkiterativefdk/rtkiterativefdk.ggo

@@ -19,3 +19,8 @@ option "hannY"     - "Cut frequency for hann window in ]0,1] (0.0 disables it)"
 
 section "Projectors"
 option "fp"    f "Forward projection method" values="Joseph","CudaRayCast","JosephAttenuated","Zeng" enum no default="Joseph"
+option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction"   string  no
+option "sigmazero" - "PSF value at a distance of 0 meter of the detector"   double  no
+option "alphapsf" - "Slope of the PSF against the detector distance"   double  no
+option "inferiorclipimage" - "Value of the inferior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no
+option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no

applications/rtkosem/rtkosem.cxx

@@ -75,13 +75,6 @@ main(int argc, char * argv[])
     inputFilter = constantImageSource;
   }
 
-  OutputImageType::Pointer attenuationMap;
-  if (args_info.attenuationmap_given)
-  {
-    // Read an existing image to initialize the attenuation map
-    attenuationMap = itk::ReadImage<OutputImageType>(args_info.attenuationmap_arg);
-  }
-
   // OSEM reconstruction filter
   rtk::OSEMConeBeamReconstructionFilter<OutputImageType>::Pointer osem =
     rtk::OSEMConeBeamReconstructionFilter<OutputImageType>::New();
@@ -91,15 +84,9 @@ main(int argc, char * argv[])
   SetBackProjectionFromGgo(args_info, osem.GetPointer());
   osem->SetInput(inputFilter->GetOutput());
   osem->SetInput(1, reader->GetOutput());
-  if (args_info.attenuationmap_given)
-    osem->SetInput(2, attenuationMap);
-  if (args_info.sigmazero_given)
-    osem->SetSigmaZero(args_info.sigmazero_arg);
-  if (args_info.alphapsf_given)
-    osem->SetAlpha(args_info.alphapsf_arg);
+  osem->SetGeometry(geometry);
   if (args_info.betaregularization_given)
     osem->SetBetaRegularization(args_info.betaregularization_arg);
-  osem->SetGeometry(geometry);
 
   osem->SetNumberOfIterations(args_info.niterations_arg);
   osem->SetNumberOfProjectionsPerSubset(args_info.nprojpersubset_arg);

applications/rtkprojectors_section.ggo

@@ -4,3 +4,5 @@ option "bp"    b "Back projection method" values="VoxelBasedBackProjection","Jos
 option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction"   string  no 
 option "sigmazero" - "PSF value at a distance of 0 meter of the detector"   double  no
 option "alphapsf" - "Slope of the PSF against the detector distance"   double  no 
+option "inferiorclipimage" - "Value of the inferior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no
+option "superiorclipimage" - "Value of the superior clip of the ray for each pixel of the projections (only with Joseph-based projector)" string no

include/rtkGgoFunctions.h

@@ -281,6 +281,14 @@ template <class TArgsInfo, class TIterativeReconstructionFilter>
 void
 SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFilter * recon)
 {
+  typename TIterativeReconstructionFilter::VolumeType::Pointer attenuationMap;
+  using VolumeType = typename TIterativeReconstructionFilter::VolumeType;
+  if (args_info.attenuationmap_given)
+  {
+    // Read an existing image to initialize the attenuation map
+    attenuationMap = itk::ReadImage<VolumeType>(args_info.attenuationmap_arg);
+  }
+
   switch (args_info.bp_arg)
   {
     case (-1): // bp__NULL, keep default
@@ -299,9 +307,17 @@ SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFi
       break;
     case (4): // bp_arg_JosephAttenuated
       recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPHATTENUATED);
+      if (args_info.attenuationmap_given)
+        recon->SetAttenuationMap(attenuationMap);
       break;
     case (5): // bp_arg_RotationBased
       recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_ZENG);
+      if (args_info.sigmazero_given)
+        recon->SetSigmaZero(args_info.sigmazero_arg);
+      if (args_info.alphapsf_given)
+        recon->SetAlphaPSF(args_info.alphapsf_arg);
+      if (args_info.attenuationmap_given)
+        recon->SetAttenuationMap(attenuationMap);
       break;
   }
 }
@@ -318,21 +334,57 @@ template <class TArgsInfo, class TIterativeReconstructionFilter>
 void
 SetForwardProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFilter * recon)
 {
+  typename TIterativeReconstructionFilter::VolumeType::Pointer attenuationMap;
+  using VolumeType = typename TIterativeReconstructionFilter::VolumeType;
+  if (args_info.attenuationmap_given)
+  {
+    // Read an existing image to initialize the attenuation map
+    attenuationMap = itk::ReadImage<VolumeType>(args_info.attenuationmap_arg);
+  }
+  using ClipImageType = itk::Image<double, TIterativeReconstructionFilter::VolumeType::ImageDimension>;
+  typename ClipImageType::Pointer inferiorClipImage, superiorClipImage;
+  if (args_info.inferiorclipimage_given)
+  {
+    // Read an existing image to initialize the attenuation map
+    inferiorClipImage = itk::ReadImage<ClipImageType>(args_info.inferiorclipimage_arg);
+  }
+  if (args_info.superiorclipimage_given)
+  {
+    // Read an existing image to initialize the attenuation map
+    superiorClipImage = itk::ReadImage<ClipImageType>(args_info.superiorclipimage_arg);
+  }
+
   switch (args_info.fp_arg)
   {
     case (-1): // fp__NULL, keep default
       break;
     case (0): // fp_arg_Joseph
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPH);
+      if (args_info.superiorclipimage_given)
+        recon->SetSuperiorClipImage(superiorClipImage);
+      if (args_info.inferiorclipimage_given)
+        recon->SetInferiorClipImage(inferiorClipImage);
       break;
     case (1): // fp_arg_CudaRayCast
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDARAYCAST);
       break;
     case (2): // fp_arg_JosephAttenuated
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPHATTENUATED);
+      if (args_info.superiorclipimage_given)
+        recon->SetSuperiorClipImage(superiorClipImage);
+      if (args_info.inferiorclipimage_given)
+        recon->SetInferiorClipImage(inferiorClipImage);
+      if (args_info.attenuationmap_given)
+        recon->SetAttenuationMap(attenuationMap);
       break;
     case (3): // fp_arg_RotationBased
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_ZENG);
+      if (args_info.sigmazero_given)
+        recon->SetSigmaZero(args_info.sigmazero_arg);
+      if (args_info.alphapsf_given)
+        recon->SetAlphaPSF(args_info.alphapsf_arg);
+      if (args_info.attenuationmap_given)
+        recon->SetAttenuationMap(attenuationMap);
       break;
   }
 }

include/rtkIterativeConeBeamReconstructionFilter.h

@@ -68,6 +68,7 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
 
   /** Convenient type alias */
   using VolumeType = ProjectionStackType;
+  using TClipImageType = itk::Image<double, VolumeType::ImageDimension>;
   typedef enum
   {
     FP_JOSEPH = 0,
@@ -113,6 +114,58 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
     return m_CurrentBackProjectionConfiguration;
   }
 
+  /** Set/Get the attenuation map for SPECT reconstruction.
+   * */
+  void
+  SetAttenuationMap(const VolumeType * attenuationMap)
+  {
+    // Process object is not const-correct so the const casting is required.
+    this->SetNthInput(2, const_cast<VolumeType *>(attenuationMap));
+  }
+  typename VolumeType::ConstPointer
+  GetAttenuationMap()
+  {
+    return static_cast<const VolumeType *>(this->itk::ProcessObject::GetInput(2));
+  }
+
+  /** Set/Get the inferior clip image. Each pixel of the image
+   ** corresponds to the value of the inferior clip of the ray
+   ** emitted from that pixel. */
+  void
+  SetInferiorClipImage(const TClipImageType * inferiorClipImage)
+  {
+    // Process object is not const-correct so the const casting is required.
+    this->SetInput("InferiorClipImage", const_cast<TClipImageType *>(inferiorClipImage));
+  }
+  typename TClipImageType::ConstPointer
+  GetInferiorClipImage()
+  {
+    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("InferiorClipImage"));
+  }
+
+  /** Set/Get the superior clip image. Each pixel of the image
+   ** corresponds to the value of the superior clip of the ray
+   ** emitted from that pixel. */
+  void
+  SetSuperiorClipImage(const TClipImageType * superiorClipImage)
+  {
+    // Process object is not const-correct so the const casting is required.
+    this->SetInput("SuperiorClipImage", const_cast<TClipImageType *>(superiorClipImage));
+  }
+  typename TClipImageType::ConstPointer
+  GetSuperiorClipImage()
+  {
+    return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("SuperiorClipImage"));
+  }
+
+  /** Get / Set the sigma zero of the PSF. Default is 1.5417233052142099 */
+  itkGetMacro(SigmaZero, double);
+  itkSetMacro(SigmaZero, double);
+
+  /** Get / Set the alpha of the PSF. Default is 0.016241189545787734 */
+  itkGetMacro(AlphaPSF, double);
+  itkSetMacro(AlphaPSF, double);
+
 protected:
   IterativeConeBeamReconstructionFilter();
   ~IterativeConeBeamReconstructionFilter() override = default;
@@ -136,6 +189,10 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
    */
   std::default_random_engine m_DefaultRandomEngine = std::default_random_engine{};
 
+  /** PSF correction coefficients */
+  double m_SigmaZero{ 1.5417233052142099 };
+  double m_AlphaPSF{ 0.016241189545787734 };
+
   /** Instantiate forward and back projectors using SFINAE. */
   using CPUImageType =
     typename itk::Image<typename ProjectionStackType::PixelType, ProjectionStackType::ImageDimension>;
@@ -207,6 +264,27 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   {
     ForwardProjectionPointerType fw;
     fw = JosephForwardAttenuatedProjectionImageFilter<VolumeType, ProjectionStackType>::New();
+    if (this->GetAttenuationMap().IsNotNull())
+    {
+      fw->SetInput(2, this->GetAttenuationMap());
+    }
+    else
+    {
+      itkExceptionMacro(<< "Set Joseph attenuated forward projection filter but no attenuation map is given");
+      return nullptr;
+    }
+    if (this->GetSuperiorClipImage().IsNotNull())
+    {
+      dynamic_cast<rtk::JosephForwardAttenuatedProjectionImageFilter<VolumeType, ProjectionStackType> *>(
+        fw.GetPointer())
+        ->SetSuperiorClipImage(this->GetSuperiorClipImage());
+    }
+    if (this->GetInferiorClipImage().IsNotNull())
+    {
+      dynamic_cast<rtk::JosephForwardAttenuatedProjectionImageFilter<VolumeType, ProjectionStackType> *>(
+        fw.GetPointer())
+        ->SetInferiorClipImage(this->GetInferiorClipImage());
+    }
     return fw;
   }
 
@@ -225,6 +303,14 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   {
     ForwardProjectionPointerType fw;
     fw = ZengForwardProjectionImageFilter<VolumeType, ProjectionStackType>::New();
+    if (this->GetAttenuationMap().IsNotNull())
+    {
+      fw->SetInput(2, this->GetAttenuationMap());
+    }
+    dynamic_cast<rtk::ZengForwardProjectionImageFilter<VolumeType, ProjectionStackType> *>(fw.GetPointer())
+      ->SetSigmaZero(m_SigmaZero);
+    dynamic_cast<rtk::ZengForwardProjectionImageFilter<VolumeType, ProjectionStackType> *>(fw.GetPointer())
+      ->SetAlpha(m_AlphaPSF);
     return fw;
   }
 
@@ -286,7 +372,16 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   {
     BackProjectionPointerType bp;
     bp = JosephBackAttenuatedProjectionImageFilter<ImageType, ImageType>::New();
-    return bp;
+    if (this->GetAttenuationMap().IsNotNull())
+    {
+      bp->SetInput(2, this->GetAttenuationMap());
+      return bp;
+    }
+    else
+    {
+      itkExceptionMacro(<< "Set Joseph attenuated backprojection filter but no attenuation map is given");
+      return nullptr;
+    }
   }
 
   template <typename ImageType, EnableVectorType<ImageType> * = nullptr>
@@ -304,6 +399,14 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   {
     BackProjectionPointerType bp;
     bp = ZengBackProjectionImageFilter<ImageType, ImageType>::New();
+    if (this->GetAttenuationMap().IsNotNull())
+    {
+      bp->SetInput(2, this->GetAttenuationMap());
+    }
+    dynamic_cast<rtk::ZengBackProjectionImageFilter<VolumeType, ProjectionStackType> *>(bp.GetPointer())
+      ->SetSigmaZero(m_SigmaZero);
+    dynamic_cast<rtk::ZengBackProjectionImageFilter<VolumeType, ProjectionStackType> *>(bp.GetPointer())
+      ->SetAlpha(m_AlphaPSF);
     return bp;
   }
 

include/rtkIterativeConeBeamReconstructionFilter.hxx

@@ -39,6 +39,16 @@ IterativeConeBeamReconstructionFilter<TOutputImage, ProjectionStackType>::Instan
   {
     case (FP_JOSEPH):
       fw = JosephForwardProjectionImageFilter<VolumeType, ProjectionStackType>::New();
+      if (this->GetSuperiorClipImage().IsNotNull())
+      {
+        dynamic_cast<rtk::JosephForwardProjectionImageFilter<VolumeType, ProjectionStackType> *>(fw.GetPointer())
+          ->SetSuperiorClipImage(this->GetSuperiorClipImage());
+      }
+      if (this->GetInferiorClipImage().IsNotNull())
+      {
+        dynamic_cast<rtk::JosephForwardProjectionImageFilter<VolumeType, ProjectionStackType> *>(fw.GetPointer())
+          ->SetInferiorClipImage(this->GetInferiorClipImage());
+      }
       break;
     case (FP_CUDARAYCAST):
       fw = InstantiateCudaForwardProjection<ProjectionStackType>();

include/rtkJosephForwardProjectionImageFilter.h

@@ -306,7 +306,7 @@ class ITK_TEMPLATE_EXPORT JosephForwardProjectionImageFilter
   ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override;
 
   /** If a third input is given, it should be in the same physical space
-   * than the first one. */
+   * as the first one. */
   void
   VerifyInputInformation() const override;
 

include/rtkJosephForwardProjectionImageFilter.hxx

@@ -243,9 +243,8 @@ JosephForwardProjectionImageFilter<TInputImage,
   typename BoxShape::VectorType stepMM, np, fp;
   for (unsigned int pix = 0; pix < outputRegionForThread.GetNumberOfPixels(); pix++, itIn->Next(), ++itOut)
   {
-    typename InputRegionIterator::PointType  pixelPosition = itIn->GetPixelPosition();
-    typename InputRegionIterator::PointType  dirVox = -itIn->GetSourceToPixel();
-    typename InputRegionIterator ::IndexType pixelIndex = itIn->GetIndex();
+    typename InputRegionIterator::PointType pixelPosition = itIn->GetPixelPosition();
+    typename InputRegionIterator::PointType dirVox = -itIn->GetSourceToPixel();
 
     // Select main direction
     unsigned int                  mainDir = 0;

include/rtkOSEMConeBeamReconstructionFilter.h

@@ -163,14 +163,6 @@ class ITK_TEMPLATE_EXPORT OSEMConeBeamReconstructionFilter
   itkGetMacro(NumberOfProjectionsPerSubset, unsigned int);
   itkSetMacro(NumberOfProjectionsPerSubset, unsigned int);
 
-  /** Get / Set the sigma zero of the PSF. Default is 1.5417233052142099 */
-  itkGetMacro(SigmaZero, double);
-  itkSetMacro(SigmaZero, double);
-
-  /** Get / Set the alpha of the PSF. Default is 0.016241189545787734 */
-  itkGetMacro(Alpha, double);
-  itkSetMacro(Alpha, double);
-
   /** Get / Set the hyperparameter for the regularization. Default is 0.01 */
   itkGetMacro(BetaRegularization, double);
   itkSetMacro(BetaRegularization, double);
@@ -229,10 +221,6 @@ class ITK_TEMPLATE_EXPORT OSEMConeBeamReconstructionFilter
   /** Number of iterations */
   unsigned int m_NumberOfIterations{ 3 };
 
-  /** PSF correction coefficients */
-  double m_SigmaZero{ -1. };
-  double m_Alpha{ -1. };
-
   /** Hyperparameter for the regularization */
   double m_BetaRegularization{ 0. };