ENH: Step size accessors for Cuda ray tracing in iterative recon

This commit fixes the issue RTKConsortium#564.
arobert01 · Nov 24, 2023 · a04d0c0 · a04d0c0
1 parent 0f218d3
commit a04d0c0
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 0 deletions.
diff --git a/applications/rtkiterativefdk/rtkiterativefdk.ggo b/applications/rtkiterativefdk/rtkiterativefdk.ggo
@@ -24,3 +24,4 @@ option "sigmazero" - "PSF value at a distance of 0 meter of the detector"   doub
 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
+option "step" - "Step size along ray (for CudaRayCast only)"                double   no   default="1"
diff --git a/applications/rtkprojectors_section.ggo b/applications/rtkprojectors_section.ggo
@@ -1,6 +1,7 @@
 section "Projectors"
 option "fp"    f "Forward projection method" values="Joseph","CudaRayCast","JosephAttenuated","Zeng" enum no default="Joseph"
 option "bp"    b "Back projection method" values="VoxelBasedBackProjection","Joseph","CudaVoxelBased","CudaRayCast","JosephAttenuated", "Zeng" enum no default="VoxelBasedBackProjection"
+option "step" - "Step size along ray (for CudaRayCast only)"                double   no   default="1"
 option "attenuationmap" - "Attenuation map relative to the volume to perfom the attenuation correction (JosephAttenuated and Zeng)"   string  no
 option "sigmazero" - "PSF value at a distance of 0 meter of the detector (Zeng only)"   double  no
 option "alphapsf" - "Slope of the PSF against the detector distance (Zeng only)"   double  no

diff --git a/include/rtkGgoFunctions.h b/include/rtkGgoFunctions.h
@@ -304,6 +304,8 @@ SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFi
       break;
     case (3): // bp_arg_CudaRayCast
       recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDARAYCAST);
+      if (args_info.step_given)
+        recon->SetStepSize(args_info.step_arg);
       break;
     case (4): // bp_arg_JosephAttenuated
       recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPHATTENUATED);
@@ -367,6 +369,8 @@ SetForwardProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructio
       break;
     case (1): // fp_arg_CudaRayCast
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDARAYCAST);
+      if (args_info.step_given)
+        recon->SetStepSize(args_info.step_arg);
       break;
     case (2): // fp_arg_JosephAttenuated
       recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPHATTENUATED);

diff --git a/include/rtkIterativeConeBeamReconstructionFilter.h b/include/rtkIterativeConeBeamReconstructionFilter.h
@@ -166,6 +166,10 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   itkGetMacro(AlphaPSF, double);
   itkSetMacro(AlphaPSF, double);
 
+  /** Get / Set step size along ray (in mm). Default is 1 mm. */
+  itkGetConstMacro(StepSize, double);
+  itkSetMacro(StepSize, double);
+
 protected:
   IterativeConeBeamReconstructionFilter();
   ~IterativeConeBeamReconstructionFilter() override = default;
@@ -193,6 +197,9 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   double m_SigmaZero{ 1.5417233052142099 };
   double m_AlphaPSF{ 0.016241189545787734 };
 
+  /** Step size along ray (in mm). */
+  double m_StepSize{ 1.0 };
+
   /** Instantiate forward and back projectors using SFINAE. */
   using CPUImageType =
     typename itk::Image<typename ProjectionStackType::PixelType, ProjectionStackType::ImageDimension>;
@@ -234,6 +241,8 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
     ForwardProjectionPointerType fw;
 #ifdef RTK_USE_CUDA
     fw = CudaForwardProjectionImageFilter<ImageType, ImageType>::New();
+    dynamic_cast<rtk::CudaForwardProjectionImageFilter<ImageType, ImageType> *>(fw.GetPointer())
+      ->SetStepSize(m_StepSize);
 #endif
     return fw;
   }
@@ -343,6 +352,7 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
     BackProjectionPointerType bp;
 #ifdef RTK_USE_CUDA
     bp = CudaRayCastBackProjectionImageFilter::New();
+    dynamic_cast<rtk::CudaRayCastBackProjectionImageFilter *>(bp.GetPointer())->SetStepSize(m_StepSize);
 #endif
     return bp;
   }