Warp projectors

applications/rtkfourdrooster/rtkfourdrooster.cxx

@@ -114,6 +114,13 @@ main(int argc, char * argv[])
   // Also set the forward and back projection filters to be used
   using ROOSTERFilterType = rtk::FourDROOSTERConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
   auto rooster = ROOSTERFilterType::New();
+  if (args_info.fp_arg == fp_arg_CudaWarp || args_info.bp_arg == bp_arg_CudaWarp)
+  {
+    itkGenericExceptionMacro(
+      << "CudaWarp projectors are not currently supported by FourDROOSTER. "
+      << "They require a separate 3D projector DVF via --warp-dvf, while FourDROOSTER uses 4D motion DVFs via "
+      << "--dvf/--idvf.");
+  }
   SetForwardProjectionFromGgo(args_info, rooster.GetPointer());
   SetBackProjectionFromGgo(args_info, rooster.GetPointer());
   rooster->SetInputVolumeSeries(inputFilter->GetOutput());

applications/rtkfourdrooster/rtkfourdrooster.py

@@ -257,6 +257,12 @@ def process(args_info: argparse.Namespace):
     # Positivity
     rooster.SetPerformPositivity(not args_info.nopositivity)
 
+    if args_info.fp == "CudaWarp" or args_info.bp == "CudaWarp":
+        raise RuntimeError(
+            "CudaWarp projectors are not currently supported by FourDROOSTER. "
+            "They require a separate 3D projector DVF via --warp-dvf, while FourDROOSTER uses 4D motion DVFs via --dvf/--idvf."
+        )
+
     rtk.SetForwardProjectionFromArgParse(args_info, rooster)
     rtk.SetBackProjectionFromArgParse(args_info, rooster)
     # Motion mask

applications/rtkiterativefdk/rtkiterativefdk.ggo

@@ -15,7 +15,8 @@ option "hann"      - "Cut frequency for hann window in ]0,1] (0.0 disables it)"
 option "hannY"     - "Cut frequency for hann window in ]0,1] (0.0 disables it)"  double                       no   default="0.0"
 
 section "Projectors"
-option "fp"    f "Forward projection method" values="Joseph","CudaRayCast","JosephAttenuated","Zeng" enum no default="Joseph"
+option "fp"    f "Forward projection method" values="Joseph","CudaRayCast","JosephAttenuated","Zeng","CudaWarp" enum no default="Joseph"
+option "warp-dvf" - "3D DVF for CudaWarp projectors." string no
 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

applications/rtkprojectors_group.py

@@ -12,7 +12,7 @@ def add_rtkprojectors_group(parser):
         "--fp",
         "-f",
         help="Forward projection method",
-        choices=["Joseph", "CudaRayCast", "JosephAttenuated", "Zeng"],
+        choices=["Joseph", "CudaRayCast", "JosephAttenuated", "Zeng", "CudaWarp"],
         default="Joseph",
     )
     rtkprojectors_group.add_argument(
@@ -26,6 +26,7 @@ def add_rtkprojectors_group(parser):
             "CudaRayCast",
             "JosephAttenuated",
             "Zeng",
+            "CudaWarp",
         ],
         default="VoxelBasedBackProjection",
     )
@@ -57,6 +58,10 @@ def add_rtkprojectors_group(parser):
         "--superiorclipimage",
         help="Superior clip of the ray for each pixel of the projections (Joseph only)",
     )
+    rtkprojectors_group.add_argument(
+        "--warp-dvf",
+        help="3D DVF for CudaWarp projectors. This is distinct from the 4D motion DVF (--dvf/--idvf)",
+    )
 
 
 # Mimicks SetBackProjectionFromGgo
@@ -97,6 +102,16 @@ def SetBackProjectionFromArgParse(args_info, recon):
             recon.SetAlphaPSF(args_info.alphapsf)
         if args_info.attenuationmap is not None:
             recon.SetAttenuationMap(attenuation_map)
+    elif args_info.bp == "CudaWarp":
+        recon.SetBackProjectionFilter(ReconType.BackProjectionType_BP_CUDAWARP)
+        if getattr(args_info, "warp_dvf", None) is None:
+            raise RuntimeError(
+                "A displacement field is required with CudaWarp projectors (use --warp-dvf)"
+            )
+        displacement_field = itk.imread(args_info.warp_dvf)
+        if hasattr(itk, "CudaImage"):
+            displacement_field = itk.cuda_image_from_image(displacement_field)
+        recon.SetDisplacementField(displacement_field)
 
 
 # Mimicks SetForwardProjectionFromGgo
@@ -127,3 +142,15 @@ def SetForwardProjectionFromArgParse(args_info, recon):
             recon.SetAlphaPSF(args_info.alphapsf)
         if args_info.attenuationmap is not None:
             recon.SetAttenuationMap(attenuation_map)
+    elif args_info.fp == "CudaWarp":
+        recon.SetForwardProjectionFilter(ReconType.ForwardProjectionType_FP_CUDAWARP)
+        if args_info.step is not None:
+            recon.SetStepSize(args_info.step)
+        if getattr(args_info, "warp_dvf", None) is None:
+            raise RuntimeError(
+                "A displacement field is required with CudaWarp projectors (use --warp-dvf)"
+            )
+        displacement_field = itk.imread(args_info.warp_dvf)
+        if hasattr(itk, "CudaImage"):
+            displacement_field = itk.cuda_image_from_image(displacement_field)
+        recon.SetDisplacementField(displacement_field)

applications/rtkprojectors_section.ggo

@@ -1,9 +1,10 @@
 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 "fp"    f "Forward projection method" values="Joseph","CudaRayCast","JosephAttenuated","Zeng","CudaWarp" enum no default="Joseph"
+option "bp"    b "Back projection method" values="VoxelBasedBackProjection","Joseph","CudaVoxelBased","CudaRayCast","JosephAttenuated", "Zeng","CudaWarp" enum no default="VoxelBasedBackProjection"
 option "step" - "Step size along ray (for CudaRayCast only, default to the minimum voxel spacing of the input volume if set to 0)"                double   no   default="0"
 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
 option "inferiorclipimage" - "Inferior clip of the ray for each pixel of the projections (Joseph only)" string no
 option "superiorclipimage" - "Superior clip of the ray for each pixel of the projections (Joseph only)" string no
+option "warp-dvf" - "3D DVF for CudaWarp projectors" string no

applications/rtkspectralrooster/rtkspectralrooster.cxx

@@ -197,6 +197,13 @@ main(int argc, char * argv[])
 
   // Set the forward and back projection filters to be used
   auto rooster = rtk::FourDROOSTERConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>::New();
+  if (args_info.fp_arg == fp_arg_CudaWarp || args_info.bp_arg == bp_arg_CudaWarp)
+  {
+    itkGenericExceptionMacro(
+      << "CudaWarp projectors are not currently supported by SpectralROOSTER. "
+      << "They require a separate 3D projector DVF via --warp-dvf, while ROOSTER-based motion compensation uses 4D "
+      << "motion DVFs.");
+  }
   SetForwardProjectionFromGgo(args_info, rooster.GetPointer());
   SetBackProjectionFromGgo(args_info, rooster.GetPointer());
   rooster->SetInputVolumeSeries(input);

include/rtkFourDROOSTERConeBeamReconstructionFilter.h

@@ -233,6 +233,8 @@ class ITK_TEMPLATE_EXPORT FourDROOSTERConeBeamReconstructionFilter
     typename VolumeSeriesType::template RebindImageType<DVFVectorType, VolumeSeriesType::ImageDimension>;
   using DVFImageType =
     typename VolumeSeriesType::template RebindImageType<DVFVectorType, VolumeSeriesType::ImageDimension - 1>;
+  // Alias 4D DVF to base-class DisplacementFieldType
+  using DisplacementFieldType = DVFSequenceImageType;
 
 #ifdef RTK_USE_CUDA
   using AverageOutOfROIFilterType = std::conditional_t<std::is_same_v<VolumeSeriesType, CPUVolumeSeriesType>,

include/rtkGgoFunctions.h

@@ -412,6 +412,14 @@ SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFi
       if (args_info.attenuationmap_given)
         recon->SetAttenuationMap(attenuationMap);
       break;
+    case (6): // bp_arg_CudaWarp
+      recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDAWARP);
+      if (args_info.warp_dvf_given)
+      {
+        using DisplacementFieldType = typename TIterativeReconstructionFilter::DisplacementFieldType;
+        recon->SetDisplacementField(itk::ReadImage<DisplacementFieldType>(args_info.warp_dvf_arg));
+      }
+      break;
   }
 }
 
@@ -481,6 +489,16 @@ SetForwardProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructio
       if (args_info.attenuationmap_given)
         recon->SetAttenuationMap(attenuationMap);
       break;
+    case (4): // fp_arg_CudaWarp
+      recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDAWARP);
+      if (args_info.step_given)
+        recon->SetStepSize(args_info.step_arg);
+      if (args_info.warp_dvf_given)
+      {
+        using DisplacementFieldType = typename TIterativeReconstructionFilter::DisplacementFieldType;
+        recon->SetDisplacementField(itk::ReadImage<DisplacementFieldType>(args_info.warp_dvf_arg));
+      }
+      break;
   }
 }
 

include/rtkIterativeConeBeamReconstructionFilter.h

@@ -19,6 +19,7 @@
 #ifndef rtkIterativeConeBeamReconstructionFilter_h
 #define rtkIterativeConeBeamReconstructionFilter_h
 
+#include <itkCovariantVector.h>
 #include <itkPixelTraits.h>
 
 // Forward projection filters
@@ -70,6 +71,10 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
 
   /** Convenient type alias */
   using VolumeType = ProjectionStackType;
+  using VolumeValueType = typename itk::PixelTraits<typename VolumeType::PixelType>::ValueType;
+  // Keep the same image backend as the volume/projection stack while changing the pixel type to a 3D DVF vector.
+  using DisplacementFieldType =
+    typename VolumeType::template RebindImageType<itk::CovariantVector<VolumeValueType, 3>, VolumeType::ImageDimension>;
   using TClipImageType = itk::Image<double, VolumeType::ImageDimension>;
   using ForwardProjectionType = enum {
     FP_JOSEPH = 0,
@@ -160,6 +165,18 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
     return static_cast<const TClipImageType *>(this->itk::ProcessObject::GetInput("SuperiorClipImage"));
   }
 
+  /** Set/Get the displacement field (DVF) required by warp forward/back projectors. */
+  void
+  SetDisplacementField(const DisplacementFieldType * dvf)
+  {
+    this->SetInput("DisplacementField", const_cast<DisplacementFieldType *>(dvf));
+  }
+  typename DisplacementFieldType::ConstPointer
+  GetDisplacementField()
+  {
+    return static_cast<const DisplacementFieldType *>(this->itk::ProcessObject::GetInput("DisplacementField"));
+  }
+
   /** Get / Set the sigma zero of the PSF. Default is 1.5417233052142099 */
   itkGetMacro(SigmaZero, double);
   itkSetMacro(SigmaZero, double);
@@ -266,7 +283,16 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
 #ifdef RTK_USE_CUDA
     auto cudaWarpFw = CudaWarpForwardProjectionImageFilter::New();
     cudaWarpFw->SetStepSize(m_StepSize);
-    fw = cudaWarpFw;
+    if (this->GetDisplacementField().IsNull())
+    {
+      itkExceptionMacro(
+        << "CudaWarpForwardProjectionImageFilter requires a displacement field. Call SetDisplacementField().");
+    }
+    else
+    {
+      cudaWarpFw->SetDisplacementField(this->GetDisplacementField().GetPointer());
+    }
+    fw = cudaWarpFw.GetPointer();
 #endif
     return fw;
   }
@@ -367,7 +393,17 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter
   {
     BackProjectionPointerType bp;
 #ifdef RTK_USE_CUDA
-    bp = CudaWarpBackProjectionImageFilter::New();
+    auto cudaWarpBp = CudaWarpBackProjectionImageFilter::New();
+    if (this->GetDisplacementField().IsNull())
+    {
+      itkExceptionMacro(
+        << "CudaWarpBackProjectionImageFilter requires a displacement field. Call SetDisplacementField().");
+    }
+    else
+    {
+      cudaWarpBp->SetDisplacementField(this->GetDisplacementField().GetPointer());
+    }
+    bp = cudaWarpBp.GetPointer();
 #endif
     return bp;
   }

test/rtkconjugategradientreconstructiontest.cxx

@@ -7,6 +7,7 @@
 #ifdef USE_CUDA
 #  include "itkCudaImage.h"
 #endif
+#include <itkCovariantVector.h>
 
 /**
  * \file rtkconjugategradientreconstructiontest.cxx
@@ -139,6 +140,26 @@ rtkconjugategradientreconstructiontest(int, char *[])
 
   CheckImageQuality<OutputImageType>(conjugategradient->GetOutput(), dsl->GetOutput(), 0.08, 23, 2.0);
   std::cout << "\n\nTest PASSED! " << std::endl;
+
+  std::cout << "\n\n****** Case 3b: CUDA Warp Forward/Back Projectors ******" << std::endl;
+  using DVFPixelType = itk::CovariantVector<OutputPixelType, 3>;
+  using DVFImageType = itk::CudaImage<DVFPixelType, 3>;
+  auto dvfSource = rtk::ConstantImageSource<DVFImageType>::New();
+  dvfSource->SetOrigin(tomographySource->GetOutput()->GetOrigin());
+  dvfSource->SetSpacing(tomographySource->GetOutput()->GetSpacing());
+  dvfSource->SetSize(tomographySource->GetOutput()->GetLargestPossibleRegion().GetSize());
+  DVFPixelType dvf;
+  dvf.Fill(0.0);
+  dvfSource->SetConstant(dvf);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(dvfSource->Update())
+
+  conjugategradient->SetDisplacementField(dvfSource->GetOutput());
+  conjugategradient->SetForwardProjectionFilter(ConjugateGradientType::FP_CUDAWARP);
+  conjugategradient->SetBackProjectionFilter(ConjugateGradientType::BP_CUDAWARP);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(conjugategradient->Update());
+
+  CheckImageQuality<OutputImageType>(conjugategradient->GetOutput(), dsl->GetOutput(), 0.08, 23, 2.0);
+  std::cout << "\n\nTest PASSED! " << std::endl;
 #endif
 
   std::cout << "\n\n****** Case 4: Joseph Backprojector, weighted least squares  ******" << std::endl;

test/rtkfourdconjugategradienttest.cxx

@@ -1,6 +1,7 @@
 #include <itkImageRegionConstIterator.h>
 #include <itkJoinSeriesImageFilter.h>
 #include <itkPasteImageFilter.h>
+#include <itkCovariantVector.h>
 #include <itksys/SystemTools.hxx>
 
 #include "rtkConstantImageSource.h"
@@ -252,6 +253,29 @@ rtkfourdconjugategradienttest(int, char *[])
 
   CheckImageQuality<VolumeSeriesType>(conjugategradient->GetOutput(), join->GetOutput(), 0.4, 12, 2.0);
   std::cout << "\n\nTest PASSED! " << std::endl;
+
+  std::cout << "\n\n****** Case 3: CUDA Warp forward projector, CUDA Warp back projector, GPU interpolation and splat "
+               "******"
+            << std::endl;
+
+  using DVFPixelType = itk::CovariantVector<OutputPixelType, 3>;
+  using DVFImageType = itk::CudaImage<DVFPixelType, 3>;
+  auto dvfSource = rtk::ConstantImageSource<DVFImageType>::New();
+  dvfSource->SetOrigin(tomographySource->GetOutput()->GetOrigin());
+  dvfSource->SetSpacing(tomographySource->GetOutput()->GetSpacing());
+  dvfSource->SetSize(tomographySource->GetOutput()->GetLargestPossibleRegion().GetSize());
+  DVFPixelType dvf;
+  dvf.Fill(0.0);
+  dvfSource->SetConstant(dvf);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(dvfSource->Update())
+
+  conjugategradient->SetDisplacementField(dvfSource->GetOutput());
+  conjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAWARP);
+  conjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDAWARP);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(conjugategradient->Update());
+
+  CheckImageQuality<VolumeSeriesType>(conjugategradient->GetOutput(), join->GetOutput(), 0.4, 12, 2.0);
+  std::cout << "\n\nTest PASSED! " << std::endl;
 #endif
 
   itksys::SystemTools::RemoveFile(signalFileName);