CUDA Zeng forward and back projectors in RTK

applications/rtkbackprojections/rtkbackprojections.cxx

@@ -29,6 +29,7 @@
 #  include "rtkCudaBackProjectionImageFilter.h"
 #  include "rtkCudaFDKBackProjectionImageFilter.h"
 #  include "rtkCudaRayCastBackProjectionImageFilter.h"
+#  include "rtkCudaZengBackProjectionImageFilter.h"
 #endif
 #include "rtkCyclicDeformationImageFilter.h"
 
@@ -125,6 +126,21 @@ main(int argc, char * argv[])
       bp = zeng;
       break;
     }
+    case (bp_arg_CudaZeng):
+    {
+#ifdef RTK_USE_CUDA
+      auto zeng = rtk::CudaZengBackProjectionImageFilter::New();
+      if (args_info.sigmazero_given)
+        zeng->SetSigmaZero(args_info.sigmazero_arg);
+      if (args_info.alphapsf_given)
+        zeng->SetAlpha(args_info.alphapsf_arg);
+      bp = zeng;
+#else
+      std::cerr << "The program has not been compiled with cuda option" << std::endl;
+      return EXIT_FAILURE;
+#endif
+      break;
+    }
     case (bp_arg_CudaFDKBackProjection):
 #ifdef RTK_USE_CUDA
       bp = rtk::CudaFDKBackProjectionImageFilter::New();

applications/rtkbackprojections/rtkbackprojections.ggo

@@ -4,7 +4,7 @@ option "geometry"  g  "XML geometry file name"
 option "output"    o "Output volume file name"                                   string   yes
 
 section "Projectors"
-option "bp"    - "Backprojection method" values="VoxelBasedBackProjection","FDKBackProjection","FDKWarpBackProjection","Joseph","JosephAttenuated", "Zeng", "CudaFDKBackProjection","CudaBackProjection","CudaRayCast"  enum no default="VoxelBasedBackProjection"
+option "bp"    - "Backprojection method" values="VoxelBasedBackProjection","FDKBackProjection","FDKWarpBackProjection","Joseph","JosephAttenuated", "Zeng", "CudaZeng", "CudaFDKBackProjection","CudaBackProjection","CudaRayCast"  enum no default="VoxelBasedBackProjection"
 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/rtkbackprojections/rtkbackprojections.py

@@ -31,6 +31,7 @@ def build_parser():
             "Joseph",
             "JosephAttenuated",
             "Zeng",
+            "CudaZeng",
             "CudaFDKBackProjection",
             "CudaBackProjection",
             "CudaRayCast",
@@ -130,11 +131,19 @@ def process(args_info: argparse.Namespace):
         bp = rtk.JosephBackAttenuatedProjectionImageFilter[
             OutputImageType, OutputImageType
         ].New()
-    elif args_info.bp == "Zeng":
-        bp = rtk.ZengBackProjectionImageFilter[OutputImageType, OutputImageType].New()
-        if args_info.sigmazero:
+    elif args_info.bp in ("Zeng", "CudaZeng"):
+        if args_info.bp == "CudaZeng":
+            if not hasattr(itk, "CudaImage"):
+                print("The program has not been compiled with cuda option")
+                sys.exit(1)
+            bp = rtk.CudaZengBackProjectionImageFilter.New()
+        else:
+            bp = rtk.ZengBackProjectionImageFilter[
+                OutputImageType, OutputImageType
+            ].New()
+        if args_info.sigmazero is not None:
             bp.SetSigmaZero(args_info.sigmazero)
-        if args_info.alphapsf:
+        if args_info.alphapsf is not None:
             bp.SetAlpha(args_info.alphapsf)
 
     elif args_info.bp == "CudaFDKBackProjection":
@@ -156,7 +165,12 @@ def process(args_info: argparse.Namespace):
             print("The program has not been compiled with cuda option")
             sys.exit(1)
 
-    if args_info.bp in ["CudaFDKBackProjection", "CudaBackProjection", "CudaRayCast"]:
+    if args_info.bp in [
+        "CudaFDKBackProjection",
+        "CudaBackProjection",
+        "CudaRayCast",
+        "CudaZeng",
+    ]:
         bp.SetInput(itk.cuda_image_from_image(constantImageSource.GetOutput()))
         bp.SetInput(1, itk.cuda_image_from_image(reader.GetOutput()))
         if attenuation_map:

applications/rtkforwardprojections/rtkforwardprojections.cxx

@@ -26,6 +26,7 @@
 #include "rtkZengForwardProjectionImageFilter.h"
 #ifdef RTK_USE_CUDA
 #  include "rtkCudaForwardProjectionImageFilter.h"
+#  include "rtkCudaZengForwardProjectionImageFilter.h"
 #endif
 
 #include <itkImageFileReader.h>
@@ -153,6 +154,21 @@ main(int argc, char * argv[])
 #else
       std::cerr << "The program has not been compiled with cuda option" << std::endl;
       return EXIT_FAILURE;
+#endif
+      break;
+    }
+    case (fp_arg_CudaZeng):
+    {
+#ifdef RTK_USE_CUDA
+      auto zeng = rtk::CudaZengForwardProjectionImageFilter::New();
+      if (args_info.sigmazero_given)
+        zeng->SetSigmaZero(args_info.sigmazero_arg);
+      if (args_info.alphapsf_given)
+        zeng->SetAlpha(args_info.alphapsf_arg);
+      forwardProjection = zeng;
+#else
+      std::cerr << "The program has not been compiled with cuda option" << std::endl;
+      return EXIT_FAILURE;
 #endif
       break;
     }

applications/rtkforwardprojections/rtkforwardprojections.ggo

@@ -7,7 +7,7 @@ option "step"      s "Step size along ray (for CudaRayCast only, default to the
 option "lowmem"    l "Compute only one projection at a time"                     flag     off
 
 section "Projectors"
-option "fp"    f "Forward projection method" values="Joseph","JosephAttenuated","CudaRayCast","Zeng","MIP" enum no default="Joseph"
+option "fp"    f "Forward projection method" values="Joseph","JosephAttenuated","CudaRayCast","Zeng","CudaZeng","MIP" 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

applications/rtkforwardprojections/rtkforwardprojections.py

@@ -44,6 +44,7 @@ def build_parser():
             "JosephAttenuated",
             "CudaRayCast",
             "Zeng",
+            "CudaZeng",
             "MIP",
             "CudaWrapRayCast",
         ],
@@ -157,10 +158,16 @@ def process(args_info):
                 "Please provide it using --attenuationmap."
             )
             sys.exit(1)
-    elif args_info.fp == "Zeng":
-        forwardProjection = rtk.ZengForwardProjectionImageFilter[
-            OutputImageType, OutputImageType
-        ].New()
+    elif args_info.fp in ("Zeng", "CudaZeng"):
+        if args_info.fp == "CudaZeng":
+            if not hasattr(itk, "CudaImage"):
+                print("The program has not been compiled with cuda option")
+                sys.exit(1)
+            forwardProjection = rtk.CudaZengForwardProjectionImageFilter.New()
+        else:
+            forwardProjection = rtk.ZengForwardProjectionImageFilter[
+                OutputImageType, OutputImageType
+            ].New()
     elif args_info.fp == "MIP":
         forwardProjection = rtk.MaximumIntensityProjectionImageFilter[
             OutputImageType, OutputImageType
@@ -191,7 +198,7 @@ def process(args_info):
             print("The program has not been compiled with cuda option")
             sys.exit(1)
 
-    if args_info.fp in ["CudaWrapRayCast", "CudaRayCast"]:
+    if args_info.fp in ["CudaWrapRayCast", "CudaRayCast", "CudaZeng"]:
         forwardProjection.SetInput(
             itk.cuda_image_from_image(constantImageSource.GetOutput())
         )
@@ -218,9 +225,9 @@ def process(args_info):
             forwardProjection.SetSuperiorClipImage(superiorClipImage)
         elif args_info.fp == "MIP":
             forwardProjection.SetSuperiorClipImage(superiorClipImage)
-    if args_info.sigmazero and args_info.fp == "Zeng":
+    if args_info.sigmazero and args_info.fp in ["Zeng", "CudaZeng"]:
         forwardProjection.SetSigmaZero(args_info.sigmazero)
-    if args_info.alphapsf and args_info.fp == "Zeng":
+    if args_info.alphapsf and args_info.fp in ["Zeng", "CudaZeng"]:
         forwardProjection.SetAlpha(args_info.alphapsf)
 
     forwardProjection.SetGeometry(geometry)

applications/rtkiterativefdk/rtkiterativefdk.ggo

@@ -15,7 +15,7 @@ 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","CudaZeng" 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

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", "CudaZeng"],
         default="Joseph",
     )
     rtkprojectors_group.add_argument(
@@ -26,6 +26,7 @@ def add_rtkprojectors_group(parser):
             "CudaRayCast",
             "JosephAttenuated",
             "Zeng",
+            "CudaZeng",
         ],
         default="VoxelBasedBackProjection",
     )
@@ -37,16 +38,16 @@ def add_rtkprojectors_group(parser):
     )
     rtkprojectors_group.add_argument(
         "--attenuationmap",
-        help="Attenuation map relative to the volume to perfom the attenuation correction (JosephAttenuated and Zeng)",
+        help="Attenuation map relative to the volume to perform the attenuation correction (JosephAttenuated, Zeng and CudaZeng)",
     )
     rtkprojectors_group.add_argument(
         "--sigmazero",
-        help="PSF value at a distance of 0 meter of the detector (Zeng only)",
+        help="PSF value at a distance of 0 meter of the detector (Zeng and CudaZeng only)",
         type=float,
     )
     rtkprojectors_group.add_argument(
         "--alphapsf",
-        help="Slope of the PSF against the detector distance (Zeng only)",
+        help="Slope of the PSF against the detector distance (Zeng and CudaZeng only)",
         type=float,
     )
     rtkprojectors_group.add_argument(
@@ -89,8 +90,13 @@ def SetBackProjectionFromArgParse(args_info, recon):
         if args_info.superiorclipimage is not None:
             recon.SetSuperiorClipImage(superior_clip_image)
         recon.SetAttenuationMap(attenuation_map)
-    elif args_info.bp == "Zeng":  # bp_arg_RotationBased
-        recon.SetBackProjectionFilter(ReconType.BackProjectionType_BP_ZENG)
+    elif args_info.bp in ("Zeng", "CudaZeng"):
+        projector = (
+            ReconType.BackProjectionType_BP_ZENG
+            if args_info.bp == "Zeng"
+            else ReconType.BackProjectionType_BP_CUDAZENG
+        )
+        recon.SetBackProjectionFilter(projector)
         if args_info.sigmazero is not None:
             recon.SetSigmaZero(args_info.sigmazero)
         if args_info.alphapsf is not None:
@@ -119,8 +125,13 @@ def SetForwardProjectionFromArgParse(args_info, recon):
             ReconType.ForwardProjectionType_FP_JOSEPHATTENUATED
         )
         recon.SetAttenuationMap(attenuation_map)
-    elif args_info.fp == "Zeng":  # fp_arg_RotationBased
-        recon.SetForwardProjectionFilter(ReconType.ForwardProjectionType_FP_ZENG)
+    elif args_info.fp in ("Zeng", "CudaZeng"):
+        projector = (
+            ReconType.ForwardProjectionType_FP_ZENG
+            if args_info.fp == "Zeng"
+            else ReconType.ForwardProjectionType_FP_CUDAZENG
+        )
+        recon.SetForwardProjectionFilter(projector)
         if args_info.sigmazero is not None:
             recon.SetSigmaZero(args_info.sigmazero)
         if args_info.alphapsf is not None:

applications/rtkprojectors_section.ggo

@@ -1,9 +1,9 @@
 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","CudaZeng" enum no default="Joseph"
+option "bp"    b "Back projection method" values="VoxelBasedBackProjection","Joseph","CudaVoxelBased","CudaRayCast","JosephAttenuated", "Zeng", "CudaZeng" 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 "attenuationmap" - "Attenuation map relative to the volume to perform the attenuation correction (JosephAttenuated, Zeng and CudaZeng)"   string  no
+option "sigmazero" - "PSF value at a distance of 0 meter of the detector (Zeng and CudaZeng only)"   double  no
+option "alphapsf" - "Slope of the PSF against the detector distance (Zeng and CudaZeng 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

documentation/docs/Projectors.md

@@ -57,6 +57,7 @@ RTK supports the following forward projector implementations :
 Cuda based projector:
 - [CudaRayCast](https://www.openrtk.org/Doxygen/classrtk_1_1CudaForwardProjectionImageFilter.html)
 - [CudaWarpRayCast](https://www.openrtk.org/Doxygen/classrtk_1_1CudaWarpForwardProjectionImageFilter.html)
+- [CudaZeng](https://www.openrtk.org/Doxygen/classrtk_1_1CudaZengForwardProjectionImageFilter.html)
 
 ## Back projectors
 
@@ -95,3 +96,4 @@ Cuda based projector:
 - [CudaFDKBackProjection](https://www.openrtk.org/Doxygen/classrtk_1_1CudaFDKBackProjectionImageFilter.html)
 - [CudaRayCast](https://www.openrtk.org/Doxygen/classrtk_1_1CudaRayCastBackProjectionImageFilter.html)
 - [CudaWarpBackProjection](https://www.openrtk.org/Doxygen/classrtk_1_1CudaWarpBackProjectionImageFilter.html)
+- [CudaZeng](https://www.openrtk.org/Doxygen/classrtk_1_1CudaZengBackProjectionImageFilter.html)

include/rtkCudaZengBackProjectionImageFilter.h

@@ -0,0 +1,82 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef rtkCudaZengBackProjectionImageFilter_h
+#define rtkCudaZengBackProjectionImageFilter_h
+
+#include "rtkConfiguration.h"
+#ifdef RTK_USE_CUDA
+
+#  include "rtkZengBackProjectionImageFilter.h"
+#  include "RTKExport.h"
+#  include <itkCudaImage.h>
+#  include <itkCudaInPlaceImageFilter.h>
+#  include <itkImageBase.h>
+
+namespace rtk
+{
+
+/** \class CudaZengBackProjectionImageFilter
+ * \brief CUDA implementation of the rotation-based Zeng backprojector.
+ *
+ * The implementation reproduces the slice recursion of
+ * ZengBackProjectionImageFilter, including the depth-dependent Gaussian PSF
+ * and the optional attenuation map (input 2).
+ * See Zeng et al., IEEE Transactions on Medical Imaging, 1999,
+ * doi:10.1109/42.796285.
+ *
+ * \ingroup RTK Projector CudaImageToImageFilter
+ */
+class RTK_EXPORT CudaZengBackProjectionImageFilter
+  : public itk::CudaInPlaceImageFilter<
+      itk::CudaImage<float, 3>,
+      itk::CudaImage<float, 3>,
+      ZengBackProjectionImageFilter<itk::CudaImage<float, 3>, itk::CudaImage<float, 3>>>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(CudaZengBackProjectionImageFilter);
+  using ImageType = itk::CudaImage<float, 3>;
+  /** Metadata-only image used to build the rotated-grid transforms. Using
+   * ImageBase avoids the CUDA image factory replacing itk::Image::New(). */
+  using CPUImageType = itk::ImageBase<3>;
+  using ProjectorType = ZengBackProjectionImageFilter<ImageType, ImageType>;
+  using Self = CudaZengBackProjectionImageFilter;
+  using Superclass = itk::CudaInPlaceImageFilter<ImageType, ImageType, ProjectorType>;
+  using Pointer = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+
+  itkNewMacro(Self);
+  itkOverrideGetNameOfClassMacro(CudaZengBackProjectionImageFilter);
+
+  /** Maximum number of projections processed together. Zero selects the automatic size (default). */
+  itkSetMacro(BatchSize, unsigned int);
+  itkGetConstMacro(BatchSize, unsigned int);
+
+protected:
+  CudaZengBackProjectionImageFilter();
+  ~CudaZengBackProjectionImageFilter() override;
+  void
+  GPUGenerateData() override;
+
+private:
+  void *       m_CudaWorkspace{ nullptr };
+  unsigned int m_BatchSize{ 0 };
+};
+
+} // namespace rtk
+#endif // RTK_USE_CUDA
+#endif // rtkCudaZengBackProjectionImageFilter_h