diff --git a/applications/rtkfourdrooster/rtkfourdrooster.cxx b/applications/rtkfourdrooster/rtkfourdrooster.cxx index efe454f3e..b1dfa4fcd 100644 --- a/applications/rtkfourdrooster/rtkfourdrooster.cxx +++ b/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; 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()); diff --git a/applications/rtkfourdrooster/rtkfourdrooster.py b/applications/rtkfourdrooster/rtkfourdrooster.py index 963bec023..577a24e48 100644 --- a/applications/rtkfourdrooster/rtkfourdrooster.py +++ b/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 diff --git a/applications/rtkiterativefdk/rtkiterativefdk.ggo b/applications/rtkiterativefdk/rtkiterativefdk.ggo index aea59e452..36a5fed8f 100644 --- a/applications/rtkiterativefdk/rtkiterativefdk.ggo +++ b/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 diff --git a/applications/rtkprojectors_group.py b/applications/rtkprojectors_group.py index f1a7eca94..da4d23290 100644 --- a/applications/rtkprojectors_group.py +++ b/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) diff --git a/applications/rtkprojectors_section.ggo b/applications/rtkprojectors_section.ggo index cdab60355..1b76336bc 100644 --- a/applications/rtkprojectors_section.ggo +++ b/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 diff --git a/applications/rtkspectralrooster/rtkspectralrooster.cxx b/applications/rtkspectralrooster/rtkspectralrooster.cxx index d5bde5748..8d99cbe9d 100644 --- a/applications/rtkspectralrooster/rtkspectralrooster.cxx +++ b/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::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); diff --git a/include/rtkFourDROOSTERConeBeamReconstructionFilter.h b/include/rtkFourDROOSTERConeBeamReconstructionFilter.h index e40b01aaf..8e2e70731 100644 --- a/include/rtkFourDROOSTERConeBeamReconstructionFilter.h +++ b/include/rtkFourDROOSTERConeBeamReconstructionFilter.h @@ -233,6 +233,8 @@ class ITK_TEMPLATE_EXPORT FourDROOSTERConeBeamReconstructionFilter typename VolumeSeriesType::template RebindImageType; using DVFImageType = typename VolumeSeriesType::template RebindImageType; + // Alias 4D DVF to base-class DisplacementFieldType + using DisplacementFieldType = DVFSequenceImageType; #ifdef RTK_USE_CUDA using AverageOutOfROIFilterType = std::conditional_t, diff --git a/include/rtkGgoFunctions.h b/include/rtkGgoFunctions.h index db5ca0702..de2dafaff 100644 --- a/include/rtkGgoFunctions.h +++ b/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(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(args_info.warp_dvf_arg)); + } + break; } } diff --git a/include/rtkIterativeConeBeamReconstructionFilter.h b/include/rtkIterativeConeBeamReconstructionFilter.h index afdaa1b87..1f806f8a7 100644 --- a/include/rtkIterativeConeBeamReconstructionFilter.h +++ b/include/rtkIterativeConeBeamReconstructionFilter.h @@ -19,6 +19,7 @@ #ifndef rtkIterativeConeBeamReconstructionFilter_h #define rtkIterativeConeBeamReconstructionFilter_h +#include #include // Forward projection filters @@ -70,6 +71,10 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter /** Convenient type alias */ using VolumeType = ProjectionStackType; + using VolumeValueType = typename itk::PixelTraits::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, VolumeType::ImageDimension>; using TClipImageType = itk::Image; using ForwardProjectionType = enum { FP_JOSEPH = 0, @@ -160,6 +165,18 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter return static_cast(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(dvf)); + } + typename DisplacementFieldType::ConstPointer + GetDisplacementField() + { + return static_cast(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; } diff --git a/test/rtkconjugategradientreconstructiontest.cxx b/test/rtkconjugategradientreconstructiontest.cxx index 2e3902e6e..4c93a56fe 100644 --- a/test/rtkconjugategradientreconstructiontest.cxx +++ b/test/rtkconjugategradientreconstructiontest.cxx @@ -7,6 +7,7 @@ #ifdef USE_CUDA # include "itkCudaImage.h" #endif +#include /** * \file rtkconjugategradientreconstructiontest.cxx @@ -139,6 +140,26 @@ rtkconjugategradientreconstructiontest(int, char *[]) CheckImageQuality(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; + using DVFImageType = itk::CudaImage; + auto dvfSource = rtk::ConstantImageSource::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(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; diff --git a/test/rtkfourdconjugategradienttest.cxx b/test/rtkfourdconjugategradienttest.cxx index 2eabc9acc..56414ded1 100644 --- a/test/rtkfourdconjugategradienttest.cxx +++ b/test/rtkfourdconjugategradienttest.cxx @@ -1,6 +1,7 @@ #include #include #include +#include #include #include "rtkConstantImageSource.h" @@ -252,6 +253,29 @@ rtkfourdconjugategradienttest(int, char *[]) CheckImageQuality(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; + using DVFImageType = itk::CudaImage; + auto dvfSource = rtk::ConstantImageSource::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(conjugategradient->GetOutput(), join->GetOutput(), 0.4, 12, 2.0); + std::cout << "\n\nTest PASSED! " << std::endl; #endif itksys::SystemTools::RemoveFile(signalFileName);