From 9070a962c9d5b507b8305a7b3bc9bc38f7f9280f Mon Sep 17 00:00:00 2001 From: "C. Constanza" Date: Wed, 29 Jul 2026 22:09:32 +0200 Subject: [PATCH 1/4] ENH: Add CUDA Zeng projectors --- .../rtkbackprojections/rtkbackprojections.cxx | 16 + .../rtkbackprojections/rtkbackprojections.ggo | 2 +- .../rtkbackprojections/rtkbackprojections.py | 19 +- .../rtkforwardprojections.cxx | 16 + .../rtkforwardprojections.ggo | 2 +- .../rtkforwardprojections.py | 13 +- .../rtkiterativefdk/rtkiterativefdk.ggo | 2 +- applications/rtkprojectors_group.py | 19 +- applications/rtkprojectors_section.ggo | 10 +- documentation/docs/Projectors.md | 2 + .../rtkCudaZengBackProjectionImageFilter.h | 81 ++ .../rtkCudaZengForwardProjectionImageFilter.h | 81 ++ include/rtkCudaZengProjectionImageFilter.hcu | 57 ++ include/rtkGgoFunctions.h | 18 + ...rtkIterativeConeBeamReconstructionFilter.h | 56 +- ...kIterativeConeBeamReconstructionFilter.hxx | 6 + src/CMakeLists.txt | 3 + src/rtkCudaZengBackProjectionImageFilter.cxx | 172 +++++ ...tkCudaZengForwardProjectionImageFilter.cxx | 170 +++++ src/rtkCudaZengProjectionImageFilter.cu | 709 ++++++++++++++++++ test/CMakeLists.txt | 2 + test/rtkadjointoperatorstest.cxx | 42 +- test/rtkzengforwardprojectiontest.cxx | 15 + test/rtkzengprojectioncomparisoncudatest.cxx | 202 +++++ .../rtkCudaZengBackProjectionImageFilter.wrap | 3 + ...kCudaZengForwardProjectionImageFilter.wrap | 3 + 26 files changed, 1702 insertions(+), 19 deletions(-) create mode 100644 include/rtkCudaZengBackProjectionImageFilter.h create mode 100644 include/rtkCudaZengForwardProjectionImageFilter.h create mode 100644 include/rtkCudaZengProjectionImageFilter.hcu create mode 100644 src/rtkCudaZengBackProjectionImageFilter.cxx create mode 100644 src/rtkCudaZengForwardProjectionImageFilter.cxx create mode 100644 src/rtkCudaZengProjectionImageFilter.cu create mode 100644 test/rtkzengprojectioncomparisoncudatest.cxx create mode 100644 wrapping/rtkCudaZengBackProjectionImageFilter.wrap create mode 100644 wrapping/rtkCudaZengForwardProjectionImageFilter.wrap diff --git a/applications/rtkbackprojections/rtkbackprojections.cxx b/applications/rtkbackprojections/rtkbackprojections.cxx index d229f9b8e..caa2f8ecc 100644 --- a/applications/rtkbackprojections/rtkbackprojections.cxx +++ b/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(); diff --git a/applications/rtkbackprojections/rtkbackprojections.ggo b/applications/rtkbackprojections/rtkbackprojections.ggo index 1289bf722..b7fc6f561 100644 --- a/applications/rtkbackprojections/rtkbackprojections.ggo +++ b/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 diff --git a/applications/rtkbackprojections/rtkbackprojections.py b/applications/rtkbackprojections/rtkbackprojections.py index 46dd39ff2..11d1d26f0 100644 --- a/applications/rtkbackprojections/rtkbackprojections.py +++ b/applications/rtkbackprojections/rtkbackprojections.py @@ -31,6 +31,7 @@ def build_parser(): "Joseph", "JosephAttenuated", "Zeng", + "CudaZeng", "CudaFDKBackProjection", "CudaBackProjection", "CudaRayCast", @@ -137,6 +138,17 @@ def process(args_info: argparse.Namespace): if args_info.alphapsf: bp.SetAlpha(args_info.alphapsf) + elif args_info.bp == "CudaZeng": + if hasattr(itk, "CudaImage"): + bp = rtk.CudaZengBackProjectionImageFilter.New() + if args_info.sigmazero: + bp.SetSigmaZero(args_info.sigmazero) + if args_info.alphapsf: + bp.SetAlpha(args_info.alphapsf) + else: + print("The program has not been compiled with cuda option") + sys.exit(1) + elif args_info.bp == "CudaFDKBackProjection": if hasattr(itk, "CudaImage"): bp = rtk.CudaFDKBackProjectionImageFilter.New() @@ -156,7 +168,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: diff --git a/applications/rtkforwardprojections/rtkforwardprojections.cxx b/applications/rtkforwardprojections/rtkforwardprojections.cxx index 536e33874..55ad3d29f 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.cxx +++ b/applications/rtkforwardprojections/rtkforwardprojections.cxx @@ -26,6 +26,7 @@ #include "rtkZengForwardProjectionImageFilter.h" #ifdef RTK_USE_CUDA # include "rtkCudaForwardProjectionImageFilter.h" +# include "rtkCudaZengForwardProjectionImageFilter.h" #endif #include @@ -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; } diff --git a/applications/rtkforwardprojections/rtkforwardprojections.ggo b/applications/rtkforwardprojections/rtkforwardprojections.ggo index 096d19d62..4b7a69dd4 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.ggo +++ b/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 diff --git a/applications/rtkforwardprojections/rtkforwardprojections.py b/applications/rtkforwardprojections/rtkforwardprojections.py index 85ca23df9..789de1664 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.py +++ b/applications/rtkforwardprojections/rtkforwardprojections.py @@ -44,6 +44,7 @@ def build_parser(): "JosephAttenuated", "CudaRayCast", "Zeng", + "CudaZeng", "MIP", "CudaWrapRayCast", ], @@ -161,6 +162,12 @@ def process(args_info): forwardProjection = rtk.ZengForwardProjectionImageFilter[ OutputImageType, OutputImageType ].New() + elif args_info.fp == "CudaZeng": + if hasattr(itk, "CudaImage"): + forwardProjection = rtk.CudaZengForwardProjectionImageFilter.New() + else: + print("The program has not been compiled with cuda option") + sys.exit(1) 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) diff --git a/applications/rtkiterativefdk/rtkiterativefdk.ggo b/applications/rtkiterativefdk/rtkiterativefdk.ggo index aea59e452..ebc4be584 100644 --- a/applications/rtkiterativefdk/rtkiterativefdk.ggo +++ b/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 diff --git a/applications/rtkprojectors_group.py b/applications/rtkprojectors_group.py index f1a7eca94..85ba0d887 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", "CudaZeng"], default="Joseph", ) rtkprojectors_group.add_argument( @@ -26,6 +26,7 @@ def add_rtkprojectors_group(parser): "CudaRayCast", "JosephAttenuated", "Zeng", + "CudaZeng", ], default="VoxelBasedBackProjection", ) @@ -97,6 +98,14 @@ 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 == "CudaZeng": + recon.SetBackProjectionFilter(ReconType.BackProjectionType_BP_CUDAZENG) + if args_info.sigmazero is not None: + recon.SetSigmaZero(args_info.sigmazero) + if args_info.alphapsf is not None: + recon.SetAlphaPSF(args_info.alphapsf) + if args_info.attenuationmap is not None: + recon.SetAttenuationMap(attenuation_map) # Mimicks SetForwardProjectionFromGgo @@ -127,3 +136,11 @@ 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 == "CudaZeng": + recon.SetForwardProjectionFilter(ReconType.ForwardProjectionType_FP_CUDAZENG) + if args_info.sigmazero is not None: + recon.SetSigmaZero(args_info.sigmazero) + if args_info.alphapsf is not None: + recon.SetAlphaPSF(args_info.alphapsf) + if args_info.attenuationmap is not None: + recon.SetAttenuationMap(attenuation_map) diff --git a/applications/rtkprojectors_section.ggo b/applications/rtkprojectors_section.ggo index cdab60355..6d4de1ff7 100644 --- a/applications/rtkprojectors_section.ggo +++ b/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 diff --git a/documentation/docs/Projectors.md b/documentation/docs/Projectors.md index 488e2a107..9714b7e0a 100644 --- a/documentation/docs/Projectors.md +++ b/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) diff --git a/include/rtkCudaZengBackProjectionImageFilter.h b/include/rtkCudaZengBackProjectionImageFilter.h new file mode 100644 index 000000000..6de7519f6 --- /dev/null +++ b/include/rtkCudaZengBackProjectionImageFilter.h @@ -0,0 +1,81 @@ +/*========================================================================= + * + * 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 "rtkBackProjectionImageFilter.h" +# include "RTKExport.h" +# include +# include +# include + +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). + * + * \ingroup RTK Projector CudaImageToImageFilter + */ +class RTK_EXPORT CudaZengBackProjectionImageFilter + : public itk::CudaInPlaceImageFilter, + itk::CudaImage, + BackProjectionImageFilter, itk::CudaImage>> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(CudaZengBackProjectionImageFilter); + using ImageType = itk::CudaImage; + /** 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 = BackProjectionImageFilter; + using Self = CudaZengBackProjectionImageFilter; + using Superclass = itk::CudaInPlaceImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + itkNewMacro(Self); + itkOverrideGetNameOfClassMacro(CudaZengBackProjectionImageFilter); + + itkGetConstMacro(SigmaZero, double); + itkSetMacro(SigmaZero, double); + itkGetConstMacro(Alpha, double); + itkSetMacro(Alpha, double); + +protected: + CudaZengBackProjectionImageFilter(); + ~CudaZengBackProjectionImageFilter() override; + void + GPUGenerateData() override; + +private: + double m_SigmaZero{ 1.5417233052142099 }; + double m_Alpha{ 0.016241189545787734 }; + void * m_CudaWorkspace{ nullptr }; +}; + +} // namespace rtk +#endif // RTK_USE_CUDA +#endif // rtkCudaZengBackProjectionImageFilter_h diff --git a/include/rtkCudaZengForwardProjectionImageFilter.h b/include/rtkCudaZengForwardProjectionImageFilter.h new file mode 100644 index 000000000..8b4d2fa15 --- /dev/null +++ b/include/rtkCudaZengForwardProjectionImageFilter.h @@ -0,0 +1,81 @@ +/*========================================================================= + * + * 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 rtkCudaZengForwardProjectionImageFilter_h +#define rtkCudaZengForwardProjectionImageFilter_h + +#include "rtkConfiguration.h" +#ifdef RTK_USE_CUDA + +# include "rtkForwardProjectionImageFilter.h" +# include "RTKExport.h" +# include +# include +# include + +namespace rtk +{ + +/** \class CudaZengForwardProjectionImageFilter + * \brief CUDA implementation of the rotation-based Zeng forward projector. + * + * The implementation reproduces the slice recursion of + * ZengForwardProjectionImageFilter, including the depth-dependent Gaussian + * PSF and the optional attenuation map (input 2). + * + * \ingroup RTK Projector CudaImageToImageFilter + */ +class RTK_EXPORT CudaZengForwardProjectionImageFilter + : public itk::CudaInPlaceImageFilter, + itk::CudaImage, + ForwardProjectionImageFilter, itk::CudaImage>> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(CudaZengForwardProjectionImageFilter); + using ImageType = itk::CudaImage; + /** 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 = ForwardProjectionImageFilter; + using Self = CudaZengForwardProjectionImageFilter; + using Superclass = itk::CudaInPlaceImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + itkNewMacro(Self); + itkOverrideGetNameOfClassMacro(CudaZengForwardProjectionImageFilter); + + itkGetConstMacro(SigmaZero, double); + itkSetMacro(SigmaZero, double); + itkGetConstMacro(Alpha, double); + itkSetMacro(Alpha, double); + +protected: + CudaZengForwardProjectionImageFilter(); + ~CudaZengForwardProjectionImageFilter() override; + void + GPUGenerateData() override; + +private: + double m_SigmaZero{ 1.5417233052142099 }; + double m_Alpha{ 0.016241189545787734 }; + void * m_CudaWorkspace{ nullptr }; +}; + +} // namespace rtk +#endif // RTK_USE_CUDA +#endif // rtkCudaZengForwardProjectionImageFilter_h diff --git a/include/rtkCudaZengProjectionImageFilter.hcu b/include/rtkCudaZengProjectionImageFilter.hcu new file mode 100644 index 000000000..6148f1dd3 --- /dev/null +++ b/include/rtkCudaZengProjectionImageFilter.hcu @@ -0,0 +1,57 @@ +/*========================================================================= + * + * 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 rtkCudaZengProjectionImageFilter_hcu +#define rtkCudaZengProjectionImageFilter_hcu + +#include "RTKExport.h" + +void RTK_EXPORT +CUDA_zeng_forward_project(const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + const float *rotatedToVolumeMatrices, + const float *farDistances, + const float *devProjectionIn, + float * devProjectionOut, + const float *devVolume, + const float *devAttenuation, + float sigmaZero, + float alpha, + void ** workspace); + +void RTK_EXPORT +CUDA_zeng_back_project(const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + const float *volumeToRotatedMatrices, + const float *nearDistances, + const int * firstSlices, + const float *devVolumeIn, + float * devVolumeOut, + const float *devProjections, + const float *devAttenuation, + float sigmaZero, + float alpha, + void ** workspace); + +void RTK_EXPORT +CUDA_zeng_release_workspace(void * workspace); + +#endif // rtkCudaZengProjectionImageFilter_hcu diff --git a/include/rtkGgoFunctions.h b/include/rtkGgoFunctions.h index db5ca0702..a1b4fec5e 100644 --- a/include/rtkGgoFunctions.h +++ b/include/rtkGgoFunctions.h @@ -412,6 +412,15 @@ SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFi if (args_info.attenuationmap_given) recon->SetAttenuationMap(attenuationMap); break; + case (6): // bp_arg_CudaZeng + recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDAZENG); + 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; } } @@ -481,6 +490,15 @@ SetForwardProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructio if (args_info.attenuationmap_given) recon->SetAttenuationMap(attenuationMap); break; + case (4): // fp_arg_CudaZeng + recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDAZENG); + 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; } } diff --git a/include/rtkIterativeConeBeamReconstructionFilter.h b/include/rtkIterativeConeBeamReconstructionFilter.h index afdaa1b87..6402fbe9c 100644 --- a/include/rtkIterativeConeBeamReconstructionFilter.h +++ b/include/rtkIterativeConeBeamReconstructionFilter.h @@ -35,6 +35,8 @@ # include "rtkCudaRayCastBackProjectionImageFilter.h" # include "rtkCudaWarpBackProjectionImageFilter.h" # include "rtkCudaWarpForwardProjectionImageFilter.h" +# include "rtkCudaZengBackProjectionImageFilter.h" +# include "rtkCudaZengForwardProjectionImageFilter.h" #endif #include @@ -76,7 +78,8 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter FP_CUDARAYCAST = 2, FP_JOSEPHATTENUATED = 3, FP_ZENG = 4, - FP_CUDAWARP = 5 + FP_CUDAWARP = 5, + FP_CUDAZENG = 6 }; using BackProjectionType = enum { BP_VOXELBASED = 0, @@ -85,7 +88,8 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter BP_CUDARAYCAST = 4, BP_JOSEPHATTENUATED = 5, BP_ZENG = 6, - BP_CUDAWARP = 7 + BP_CUDAWARP = 7, + BP_CUDAZENG = 8 }; /** Typedefs of each subfilter of this composite filter */ @@ -324,6 +328,30 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter return nullptr; } + template * = nullptr> + ForwardProjectionPointerType + InstantiateCudaZengForwardProjection() + { + ForwardProjectionPointerType fw; +#ifdef RTK_USE_CUDA + auto zeng = CudaZengForwardProjectionImageFilter::New(); + if (this->GetAttenuationMap().IsNotNull()) + zeng->SetInput(2, this->GetAttenuationMap()); + zeng->SetSigmaZero(m_SigmaZero); + zeng->SetAlpha(m_AlphaPSF); + fw = zeng; +#endif + return fw; + } + + template * = nullptr> + ForwardProjectionPointerType + InstantiateCudaZengForwardProjection() + { + itkGenericExceptionMacro(<< "CudaZengForwardProjectionImageFilter only available with 3D CudaImage."); + return nullptr; + } + template * = nullptr> ForwardProjectionPointerType @@ -440,6 +468,30 @@ class ITK_TEMPLATE_EXPORT IterativeConeBeamReconstructionFilter return nullptr; } + template * = nullptr> + BackProjectionPointerType + InstantiateCudaZengBackProjection() + { + BackProjectionPointerType bp; +#ifdef RTK_USE_CUDA + auto zeng = CudaZengBackProjectionImageFilter::New(); + if (this->GetAttenuationMap().IsNotNull()) + zeng->SetInput(2, this->GetAttenuationMap()); + zeng->SetSigmaZero(m_SigmaZero); + zeng->SetAlpha(m_AlphaPSF); + bp = zeng; +#endif + return bp; + } + + template * = nullptr> + BackProjectionPointerType + InstantiateCudaZengBackProjection() + { + itkGenericExceptionMacro(<< "CudaZengBackProjectionImageFilter only available with 3D CudaImage."); + return nullptr; + } + template * = nullptr> BackProjectionPointerType diff --git a/include/rtkIterativeConeBeamReconstructionFilter.hxx b/include/rtkIterativeConeBeamReconstructionFilter.hxx index ae210811e..bcf080a49 100644 --- a/include/rtkIterativeConeBeamReconstructionFilter.hxx +++ b/include/rtkIterativeConeBeamReconstructionFilter.hxx @@ -62,6 +62,9 @@ IterativeConeBeamReconstructionFilter::Instan case (FP_CUDAWARP): fw = InstantiateCudaWarpForwardProjection(); break; + case (FP_CUDAZENG): + fw = InstantiateCudaZengForwardProjection(); + break; default: itkGenericExceptionMacro(<< "Unhandled --fp value."); } @@ -96,6 +99,9 @@ IterativeConeBeamReconstructionFilter::Instan case (BP_CUDAWARP): bp = InstantiateCudaWarpBackProjection(); break; + case (BP_CUDAZENG): + bp = InstantiateCudaZengBackProjection(); + break; default: itkGenericExceptionMacro(<< "Unhandled --bp value."); } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 21a70b865..18bc565fc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -96,6 +96,8 @@ if(RTK_USE_CUDA) rtkCudaParkerShortScanImageFilter.cxx rtkCudaPolynomialGainCorrectionImageFilter.cxx rtkCudaRayCastBackProjectionImageFilter.cxx + rtkCudaZengBackProjectionImageFilter.cxx + rtkCudaZengForwardProjectionImageFilter.cxx rtkCudaScatterGlareCorrectionImageFilter.cxx rtkCudaSplatImageFilter.cxx rtkCudaTotalVariationDenoisingBPDQImageFilter.cxx @@ -128,6 +130,7 @@ if(RTK_USE_CUDA) rtkCudaParkerShortScanImageFilter.cu rtkCudaPolynomialGainCorrectionImageFilter.cu rtkCudaRayCastBackProjectionImageFilter.cu + rtkCudaZengProjectionImageFilter.cu rtkCudaSplatImageFilter.cu rtkCudaTotalVariationDenoisingBPDQImageFilter.cu rtkCudaUtilities.cu diff --git a/src/rtkCudaZengBackProjectionImageFilter.cxx b/src/rtkCudaZengBackProjectionImageFilter.cxx new file mode 100644 index 000000000..40a7cbcf6 --- /dev/null +++ b/src/rtkCudaZengBackProjectionImageFilter.cxx @@ -0,0 +1,172 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +#include "rtkCudaZengBackProjectionImageFilter.h" + +#ifdef RTK_USE_CUDA +# include "rtkCudaZengProjectionImageFilter.hcu" +# include "rtkHomogeneousMatrix.h" +# include + +namespace rtk +{ +namespace +{ +using MatrixType = ThreeDCircularProjectionGeometry::ThreeDHomogeneousMatrixType; + +MatrixType +TransformMatrix(const itk::CenteredEuler3DTransform * transform) +{ + MatrixType matrix; + matrix.SetIdentity(); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 3; ++c) + matrix[r][c] = transform->GetMatrix()[r][c]; + matrix[r][3] = transform->GetOffset()[r]; + } + return matrix; +} +void +StoreMatrix(const MatrixType & matrix, float * destination) +{ + for (unsigned int r = 0; r < 3; ++r) + for (unsigned int c = 0; c < 4; ++c) + destination[4 * r + c] = static_cast(matrix[r][c]); +} +} // namespace + +CudaZengBackProjectionImageFilter::CudaZengBackProjectionImageFilter() { this->InPlaceOff(); } + +CudaZengBackProjectionImageFilter::~CudaZengBackProjectionImageFilter() +{ + CUDA_zeng_release_workspace(m_CudaWorkspace); +} + +void +CudaZengBackProjectionImageFilter::GPUGenerateData() +{ + const auto * geometry = this->GetGeometry(); + if (!geometry) + itkGenericExceptionMacro(<< "CudaZengBackProjectionImageFilter requires a projection geometry."); + + const auto & volumeRegion = this->GetOutput()->GetBufferedRegion(); + const auto & projectionRegion = this->GetInput(1)->GetBufferedRegion(); + const unsigned int firstProjection = this->GetInput(1)->GetRequestedRegion().GetIndex(2); + const unsigned int numberOfProjections = this->GetInput(1)->GetRequestedRegion().GetSize(2); + int projectionSize[3] = { static_cast(projectionRegion.GetSize(0)), + static_cast(projectionRegion.GetSize(1)), + static_cast(numberOfProjections) }; + int volumeSize[3] = { static_cast(volumeRegion.GetSize(0)), + static_cast(volumeRegion.GetSize(1)), + static_cast(volumeRegion.GetSize(2)) }; + int rotatedSize[3] = { projectionSize[0], + projectionSize[1], + static_cast(std::ceil(volumeRegion.GetSize(2) * std::sqrt(2.0))) }; + float rotatedSpacing[3] = { static_cast(this->GetInput(1)->GetSpacing()[0]), + static_cast(this->GetInput(1)->GetSpacing()[1]), + static_cast(this->GetInput(0)->GetSpacing()[2]) }; + + itk::ContinuousIndex centerIndex; + for (unsigned int d = 0; d < 3; ++d) + centerIndex[d] = volumeRegion.GetIndex(d) + (volumeRegion.GetSize(d) - 1.0) / 2.0; + CPUImageType::PointType volumeCenter; + this->GetInput(0)->TransformContinuousIndexToPhysicalPoint(centerIndex, volumeCenter); + + auto rotatedImage = CPUImageType::New(); + CPUImageType::RegionType rotatedRegion; + CPUImageType::IndexType zeroIndex{}; + CPUImageType::SizeType rotatedImageSize; + CPUImageType::SpacingType spacing; + for (unsigned int d = 0; d < 3; ++d) + { + rotatedImageSize[d] = rotatedSize[d]; + spacing[d] = rotatedSpacing[d]; + } + rotatedRegion.SetIndex(zeroIndex); + rotatedRegion.SetSize(rotatedImageSize); + rotatedImage->SetRegions(rotatedRegion); + rotatedImage->SetSpacing(spacing); + rotatedImage->SetDirection(this->GetInput(1)->GetDirection()); + + std::vector matrices(12 * numberOfProjections); + std::vector nearDistances(numberOfProjections); + std::vector firstSlices(numberOfProjections); + MatrixType volumeIndexTranslation; + volumeIndexTranslation.SetIdentity(); + for (unsigned int d = 0; d < 3; ++d) + volumeIndexTranslation[d][3] = volumeRegion.GetIndex(d); + using TransformType = itk::CenteredEuler3DTransform; + auto transform = TransformType::New(); + TransformType::InputPointType zero{}; + transform->SetCenter(zero); + + for (unsigned int local = 0; local < numberOfProjections; ++local) + { + const unsigned int projection = firstProjection + local; + const double angle = geometry->GetGantryAngles()[projection]; + transform->SetRotation(0., angle, 0.); + transform->SetTranslation(itk::MakeVector(geometry->GetProjectionOffsetsX()[projection] * std::cos(-angle), + geometry->GetProjectionOffsetsY()[projection], + geometry->GetProjectionOffsetsX()[projection] * std::sin(-angle))); + const auto rotatedCenter = transform->GetMatrix() * volumeCenter; + CPUImageType::PointType origin; + origin[0] = this->GetInput(1)->GetOrigin()[0]; + origin[1] = this->GetInput(1)->GetOrigin()[1]; + origin[2] = rotatedCenter[2] - spacing[2] * (rotatedSize[2] - 1.0) / 2.0; + rotatedImage->SetOrigin(origin); + + auto inverse = TransformType::New(); + if (!transform->GetInverse(inverse)) + itkGenericExceptionMacro(<< "Could not invert Zeng rotation transform."); + const MatrixType matrix = GetPhysicalPointToIndexMatrix(rotatedImage.GetPointer()).GetVnlMatrix() * + TransformMatrix(inverse).GetVnlMatrix() * + GetIndexToPhysicalPointMatrix(this->GetInput(0)).GetVnlMatrix() * + volumeIndexTranslation.GetVnlMatrix(); + StoreMatrix(matrix, matrices.data() + 12 * local); + + const double firstDistance = geometry->GetSourceToIsocenterDistances()[projection] + origin[2]; + firstSlices[local] = std::max(0, static_cast(std::ceil(-firstDistance / spacing[2]))); + nearDistances[local] = static_cast(firstDistance + firstSlices[local] * spacing[2]); + } + + const auto projectionOffset = firstProjection - projectionRegion.GetIndex(2); + const auto pixelsPerProjection = projectionSize[0] * projectionSize[1]; + const float * projections = static_cast(this->GetInput(1)->GetCudaDataManager()->GetGPUBufferPointer()) + + projectionOffset * pixelsPerProjection; + const float * volumeIn = static_cast(this->GetInput(0)->GetCudaDataManager()->GetGPUBufferPointer()); + float * volumeOut = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); + const float * attenuation = + this->GetInput(2) ? static_cast(this->GetInput(2)->GetCudaDataManager()->GetGPUBufferPointer()) : nullptr; + CUDA_zeng_back_project(projectionSize, + volumeSize, + rotatedSize, + rotatedSpacing, + matrices.data(), + nearDistances.data(), + firstSlices.data(), + volumeIn, + volumeOut, + projections, + attenuation, + static_cast(m_SigmaZero), + static_cast(m_Alpha), + &m_CudaWorkspace); +} + +} // namespace rtk +#endif diff --git a/src/rtkCudaZengForwardProjectionImageFilter.cxx b/src/rtkCudaZengForwardProjectionImageFilter.cxx new file mode 100644 index 000000000..9d763f1cc --- /dev/null +++ b/src/rtkCudaZengForwardProjectionImageFilter.cxx @@ -0,0 +1,170 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +#include "rtkCudaZengForwardProjectionImageFilter.h" + +#ifdef RTK_USE_CUDA +# include "rtkCudaZengProjectionImageFilter.hcu" +# include "rtkHomogeneousMatrix.h" +# include + +namespace rtk +{ +namespace +{ +using MatrixType = ThreeDCircularProjectionGeometry::ThreeDHomogeneousMatrixType; + +MatrixType +TransformMatrix(const itk::CenteredEuler3DTransform * transform) +{ + MatrixType matrix; + matrix.SetIdentity(); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 3; ++c) + matrix[r][c] = transform->GetMatrix()[r][c]; + matrix[r][3] = transform->GetOffset()[r]; + } + return matrix; +} + +void +StoreTextureMatrix(const MatrixType & matrix, const itk::Index<3> & bufferedIndex, float * destination) +{ + for (unsigned int r = 0; r < 3; ++r) + for (unsigned int c = 0; c < 4; ++c) + destination[4 * r + c] = static_cast(matrix[r][c]); + destination[3] -= bufferedIndex[0]; + destination[7] -= bufferedIndex[1]; + destination[11] -= bufferedIndex[2]; +} +} // namespace + +CudaZengForwardProjectionImageFilter::CudaZengForwardProjectionImageFilter() { this->InPlaceOff(); } + +CudaZengForwardProjectionImageFilter::~CudaZengForwardProjectionImageFilter() +{ + CUDA_zeng_release_workspace(m_CudaWorkspace); +} + +void +CudaZengForwardProjectionImageFilter::GPUGenerateData() +{ + const auto * geometry = this->GetGeometry(); + if (!geometry) + itkGenericExceptionMacro(<< "CudaZengForwardProjectionImageFilter requires a projection geometry."); + + constexpr unsigned int Dimension = 3; + const auto & projectionRegion = this->GetOutput()->GetBufferedRegion(); + const auto & volumeRegion = this->GetInput(1)->GetBufferedRegion(); + const unsigned int firstProjection = this->GetOutput()->GetRequestedRegion().GetIndex(2); + const unsigned int numberOfProjections = this->GetOutput()->GetRequestedRegion().GetSize(2); + + int projectionSize[3] = { static_cast(projectionRegion.GetSize(0)), + static_cast(projectionRegion.GetSize(1)), + static_cast(numberOfProjections) }; + int volumeSize[3] = { static_cast(volumeRegion.GetSize(0)), + static_cast(volumeRegion.GetSize(1)), + static_cast(volumeRegion.GetSize(2)) }; + int rotatedSize[3] = { projectionSize[0], + projectionSize[1], + static_cast(std::ceil(volumeRegion.GetSize(2) * std::sqrt(2.0))) }; + float rotatedSpacing[3] = { static_cast(this->GetInput(0)->GetSpacing()[0]), + static_cast(this->GetInput(0)->GetSpacing()[1]), + static_cast(this->GetInput(1)->GetSpacing()[2]) }; + + CPUImageType::PointType volumeCenter; + itk::ContinuousIndex centerIndex; + for (unsigned int d = 0; d < Dimension; ++d) + centerIndex[d] = volumeRegion.GetIndex(d) + (volumeRegion.GetSize(d) - 1.0) / 2.0; + this->GetInput(1)->TransformContinuousIndexToPhysicalPoint(centerIndex, volumeCenter); + + std::vector matrices(12 * numberOfProjections); + std::vector farDistances(numberOfProjections); + using TransformType = itk::CenteredEuler3DTransform; + auto transform = TransformType::New(); + TransformType::InputPointType zero{}; + transform->SetCenter(zero); + + auto rotatedImage = CPUImageType::New(); + CPUImageType::SpacingType spacing; + CPUImageType::PointType origin; + CPUImageType::DirectionType direction = this->GetInput(0)->GetDirection(); + CPUImageType::RegionType rotatedRegion; + CPUImageType::IndexType rotatedIndex{}; + CPUImageType::SizeType rotatedImageSize; + for (unsigned int d = 0; d < 3; ++d) + { + spacing[d] = rotatedSpacing[d]; + rotatedImageSize[d] = rotatedSize[d]; + } + rotatedRegion.SetIndex(rotatedIndex); + rotatedRegion.SetSize(rotatedImageSize); + rotatedImage->SetRegions(rotatedRegion); + rotatedImage->SetSpacing(spacing); + rotatedImage->SetDirection(direction); + + for (unsigned int local = 0; local < numberOfProjections; ++local) + { + const unsigned int projection = firstProjection + local; + const double angle = geometry->GetGantryAngles()[projection]; + transform->SetRotation(0., angle, 0.); + transform->SetTranslation(itk::MakeVector(geometry->GetProjectionOffsetsX()[projection] * std::cos(-angle), + geometry->GetProjectionOffsetsY()[projection], + geometry->GetProjectionOffsetsX()[projection] * std::sin(-angle))); + + const auto rotatedCenter = transform->GetMatrix() * volumeCenter; + origin[0] = this->GetInput(0)->GetOrigin()[0]; + origin[1] = this->GetInput(0)->GetOrigin()[1]; + origin[2] = rotatedCenter[2] - spacing[2] * (rotatedSize[2] - 1.0) / 2.0; + rotatedImage->SetOrigin(origin); + + const MatrixType matrix = GetPhysicalPointToIndexMatrix(this->GetInput(1)).GetVnlMatrix() * + TransformMatrix(transform).GetVnlMatrix() * + GetIndexToPhysicalPointMatrix(rotatedImage.GetPointer()).GetVnlMatrix(); + StoreTextureMatrix(matrix, volumeRegion.GetIndex(), matrices.data() + 12 * local); + farDistances[local] = static_cast(geometry->GetSourceToIsocenterDistances()[projection] + origin[2] + + spacing[2] * (rotatedSize[2] - 1)); + } + + const auto projectionOffset = firstProjection - projectionRegion.GetIndex(2); + const auto pixelsPerProjection = projectionSize[0] * projectionSize[1]; + const float * projectionIn = static_cast(this->GetInput(0)->GetCudaDataManager()->GetGPUBufferPointer()) + + projectionOffset * pixelsPerProjection; + float * projectionOut = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()) + + projectionOffset * pixelsPerProjection; + const float * volume = static_cast(this->GetInput(1)->GetCudaDataManager()->GetGPUBufferPointer()); + const float * attenuation = + this->GetInput(2) ? static_cast(this->GetInput(2)->GetCudaDataManager()->GetGPUBufferPointer()) : nullptr; + + CUDA_zeng_forward_project(projectionSize, + volumeSize, + rotatedSize, + rotatedSpacing, + matrices.data(), + farDistances.data(), + projectionIn, + projectionOut, + volume, + attenuation, + static_cast(m_SigmaZero), + static_cast(m_Alpha), + &m_CudaWorkspace); +} + +} // namespace rtk +#endif diff --git a/src/rtkCudaZengProjectionImageFilter.cu b/src/rtkCudaZengProjectionImageFilter.cu new file mode 100644 index 000000000..a09cdc8a6 --- /dev/null +++ b/src/rtkCudaZengProjectionImageFilter.cu @@ -0,0 +1,709 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +#include "rtkCudaZengProjectionImageFilter.hcu" +#include "rtkCudaUtilities.hcu" + +#include +#include +#include +#include +#include +#include + +namespace +{ +inline __device__ float3 +applyMatrix(const float * m, float x, float y, float z) +{ + return make_float3(m[0] * x + m[1] * y + m[2] * z + m[3], + m[4] * x + m[5] * y + m[6] * z + m[7], + m[8] * x + m[9] * y + m[10] * z + m[11]); +} + +__device__ float +trilinearZero(const float * image, int3 size, float3 p) +{ + const int x0 = static_cast(floorf(p.x)); + const int y0 = static_cast(floorf(p.y)); + const int z0 = static_cast(floorf(p.z)); + const float fx = p.x - x0; + const float fy = p.y - y0; + const float fz = p.z - z0; + float result = 0.f; + for (int dz = 0; dz <= 1; ++dz) + for (int dy = 0; dy <= 1; ++dy) + for (int dx = 0; dx <= 1; ++dx) + { + const int x = x0 + dx, y = y0 + dy, z = z0 + dz; + if (x >= 0 && x < size.x && y >= 0 && y < size.y && z >= 0 && z < size.z) + { + const float wx = dx ? fx : 1.f - fx; + const float wy = dy ? fy : 1.f - fy; + const float wz = dz ? fz : 1.f - fz; + result += wx * wy * wz * image[(z * size.y + y) * size.x + x]; + } + } + return result; +} + +__global__ void +sampleForwardSlice(float * slice, + const float * previous, + const float * volume, + const float * attenuation, + int3 volumeSize, + const float * matrix, + int width, + int height, + int z, + float attenuationStep, + bool addPrevious) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= width || y >= height) + return; + const float3 position = applyMatrix(matrix, x, y, z); + float value = trilinearZero(volume, volumeSize, position); + if (addPrevious) + value += previous[y * width + x]; + if (attenuation) + value *= expf(-attenuationStep * trilinearZero(attenuation, volumeSize, position)); + slice[y * width + x] = value; +} + +__global__ void +gaussianXShared(const float * input, float * output, int width, int height, const float * coefficients, int radius) +{ + extern __shared__ float tile[]; + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + const int tileWidth = blockDim.x + 2 * radius; + const int tileElements = tileWidth * blockDim.y; + const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; + const int threadCount = blockDim.x * blockDim.y; + for (int index = threadIndex; index < tileElements; index += threadCount) + { + const int localX = index % tileWidth; + const int localY = index / tileWidth; + const int globalX = blockIdx.x * blockDim.x + localX - radius; + const int globalY = blockIdx.y * blockDim.y + localY; + tile[index] = globalX >= 0 && globalX < width && globalY < height ? input[globalY * width + globalX] : 0.f; + } + __syncthreads(); + if (x >= width || y >= height) + return; + float sum = 0.f; + for (int offset = -radius; offset <= radius; ++offset) + sum += coefficients[abs(offset)] * tile[threadIdx.y * tileWidth + threadIdx.x + radius + offset]; + output[y * width + x] = sum; +} + +__global__ void +gaussianYShared(const float * input, float * output, int width, int height, const float * coefficients, int radius) +{ + extern __shared__ float tile[]; + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + const int tileHeight = blockDim.y + 2 * radius; + const int tileElements = blockDim.x * tileHeight; + const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; + const int threadCount = blockDim.x * blockDim.y; + for (int index = threadIndex; index < tileElements; index += threadCount) + { + const int localX = index % blockDim.x; + const int localY = index / blockDim.x; + const int globalX = blockIdx.x * blockDim.x + localX; + const int globalY = blockIdx.y * blockDim.y + localY - radius; + tile[index] = globalX < width && globalY >= 0 && globalY < height ? input[globalY * width + globalX] : 0.f; + } + __syncthreads(); + if (x >= width || y >= height) + return; + float sum = 0.f; + for (int offset = -radius; offset <= radius; ++offset) + sum += coefficients[abs(offset)] * tile[(threadIdx.y + radius + offset) * blockDim.x + threadIdx.x]; + output[y * width + x] = sum; +} + +std::vector +discreteGaussianCoefficients(double variance) +{ + if (variance <= 1.e-12) + return { 1.f }; + + constexpr double maximumError = 1.e-5; + constexpr size_t maximumKernelWidth = 32; + const double exponential = std::exp(-variance); + std::vector coefficients; + coefficients.push_back(exponential * std::cyl_bessel_i(0., variance)); + double sum = coefficients[0]; + for (int order = 1; sum < 1. - maximumError; ++order) + { + const double coefficient = exponential * std::cyl_bessel_i(static_cast(order), variance); + coefficients.push_back(coefficient); + sum += 2. * coefficient; + if (coefficient <= 0. || coefficients.size() > maximumKernelWidth) + break; + } + std::vector normalized(coefficients.size()); + std::transform(coefficients.begin(), coefficients.end(), normalized.begin(), [sum](double value) { + return static_cast(value / sum); + }); + return normalized; +} + +constexpr size_t GaussianCoefficientStride = 33; + +struct GaussianMetadata +{ + std::vector coefficientsX; + std::vector coefficientsY; + std::vector radiiX; + std::vector radiiY; +}; + +struct MetadataKey +{ + int mode{}; + std::vector integers; + std::vector values; + + bool + operator<(const MetadataKey & other) const + { + return std::tie(mode, integers, values) < std::tie(other.mode, other.integers, other.values); + } +}; + +struct DeviceMetadata +{ + float * coefficientsX{}; + float * coefficientsY{}; + float * matrices{}; + float * inverseMatrices{}; + std::vector radiiX; + std::vector radiiY; + + ~DeviceMetadata() + { + cudaFree(inverseMatrices); + cudaFree(matrices); + cudaFree(coefficientsY); + cudaFree(coefficientsX); + } +}; + +struct CudaZengWorkspace +{ + float * current{}; + float * blurred{}; + float * gaussianScratch{}; + float * rotated{}; + size_t sliceCapacity{}; + size_t rotatedCapacity{}; + std::map metadata; + + ~CudaZengWorkspace() + { + cudaFree(rotated); + cudaFree(gaussianScratch); + cudaFree(blurred); + cudaFree(current); + } +}; + +void +ensureSliceBuffers(CudaZengWorkspace & workspace, size_t requiredElements) +{ + if (requiredElements <= workspace.sliceCapacity) + return; + cudaFree(workspace.gaussianScratch); + cudaFree(workspace.blurred); + cudaFree(workspace.current); + cudaMalloc(&workspace.current, requiredElements * sizeof(float)); + cudaMalloc(&workspace.blurred, requiredElements * sizeof(float)); + cudaMalloc(&workspace.gaussianScratch, requiredElements * sizeof(float)); + workspace.sliceCapacity = requiredElements; +} + +void +ensureRotatedBuffer(CudaZengWorkspace & workspace, size_t requiredElements) +{ + if (requiredElements <= workspace.rotatedCapacity) + return; + cudaFree(workspace.rotated); + cudaMalloc(&workspace.rotated, requiredElements * sizeof(float)); + workspace.rotatedCapacity = requiredElements; +} + +MetadataKey +makeMetadataKey(int mode, + const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + float sigmaZero, + float alpha, + const float * matrices, + const float * distances, + const int * firstSlices) +{ + MetadataKey key; + key.mode = mode; + key.integers.insert(key.integers.end(), projectionSize, projectionSize + 3); + key.integers.insert(key.integers.end(), volumeSize, volumeSize + 3); + key.integers.insert(key.integers.end(), rotatedSize, rotatedSize + 3); + if (firstSlices) + key.integers.insert(key.integers.end(), firstSlices, firstSlices + projectionSize[2]); + key.values.insert(key.values.end(), rotatedSpacing, rotatedSpacing + 3); + key.values.push_back(sigmaZero); + key.values.push_back(alpha); + key.values.insert(key.values.end(), matrices, matrices + 12 * projectionSize[2]); + key.values.insert(key.values.end(), distances, distances + projectionSize[2]); + return key; +} + +void +appendGaussianMetadata(GaussianMetadata & metadata, + std::map> & coefficientCache, + float variance, + float spacingX, + float spacingY) +{ + const double variances[] = { variance / (spacingX * spacingX), variance / (spacingY * spacingY) }; + std::vector * packed[] = { &metadata.coefficientsX, &metadata.coefficientsY }; + std::vector * radii[] = { &metadata.radiiX, &metadata.radiiY }; + for (int dimension = 0; dimension < 2; ++dimension) + { + auto [iterator, inserted] = coefficientCache.try_emplace(variances[dimension]); + if (inserted) + iterator->second = discreteGaussianCoefficients(variances[dimension]); + const auto & coefficients = iterator->second; + radii[dimension]->push_back(static_cast(coefficients.size()) - 1); + packed[dimension]->insert(packed[dimension]->end(), coefficients.begin(), coefficients.end()); + packed[dimension]->resize(packed[dimension]->size() + GaussianCoefficientStride - coefficients.size(), 0.f); + } +} + +void +gaussian2D(const float * input, + float * output, + float * scratch, + const float * deviceCoefficientsX, + const float * deviceCoefficientsY, + int radiusX, + int radiusY, + int width, + int height, + dim3 grid, + dim3 block) +{ + const size_t sharedX = (block.x + 2 * radiusX) * block.y * sizeof(float); + const size_t sharedY = block.x * (block.y + 2 * radiusY) * sizeof(float); + gaussianXShared<<>>(input, scratch, width, height, deviceCoefficientsX, radiusX); + gaussianYShared<<>>(scratch, output, width, height, deviceCoefficientsY, radiusY); +} + +__global__ void +finishForward(const float * input, const float * zeng, float * output, int count, float thickness) +{ + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < count) + output[i] = input[i] + thickness * zeng[i]; +} + +__global__ void +copyProjection(const float * projections, float * slice, int count, int projection) +{ + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < count) + slice[i] = projections[projection * count + i]; +} + +__global__ void +attenuateSlice(float * slice, + cudaTextureObject_t attenuation, + const float * rotatedToVolume, + int width, + int height, + int z, + float attenuationStep) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= width || y >= height) + return; + const float3 position = applyMatrix(rotatedToVolume, x, y, z); + slice[y * width + x] *= expf(-attenuationStep * tex3D(attenuation, position.x, position.y, position.z)); +} + +__global__ void +storeSlice(const float * slice, float * volume, int count, int z) +{ + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < count) + volume[z * count + i] = slice[i]; +} + +__global__ void +addRotatedVolume(const float * rotated, + const float * input, + float * output, + int3 volumeSize, + int3 rotatedSize, + const float * volumeToRotated, + float thickness) +{ + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + const int z = blockIdx.z * blockDim.z + threadIdx.z; + if (x >= volumeSize.x || y >= volumeSize.y || z >= volumeSize.z) + return; + const int index = (z * volumeSize.y + y) * volumeSize.x + x; + const float3 position = applyMatrix(volumeToRotated, x, y, z); + output[index] = input[index] + thickness * trilinearZero(rotated, rotatedSize, position); +} + +void +invertAffine(const float * source, float * inverse) +{ + const double a = source[0], b = source[1], c = source[2]; + const double d = source[4], e = source[5], f = source[6]; + const double g = source[8], h = source[9], i = source[10]; + const double determinant = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g); + const double inv[9] = { (e * i - f * h) / determinant, (c * h - b * i) / determinant, (b * f - c * e) / determinant, + (f * g - d * i) / determinant, (a * i - c * g) / determinant, (c * d - a * f) / determinant, + (d * h - e * g) / determinant, (b * g - a * h) / determinant, (a * e - b * d) / determinant }; + for (int r = 0; r < 3; ++r) + { + for (int col = 0; col < 3; ++col) + inverse[4 * r + col] = static_cast(inv[3 * r + col]); + inverse[4 * r + 3] = + static_cast(-(inv[3 * r] * source[3] + inv[3 * r + 1] * source[7] + inv[3 * r + 2] * source[11]) + 0.5); + } +} +} // namespace + +void +CUDA_zeng_forward_project(const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + const float * rotatedToVolumeMatrices, + const float * farDistances, + const float * devProjectionIn, + float * devProjectionOut, + const float * devVolume, + const float * devAttenuation, + float sigmaZero, + float alpha, + void ** workspacePointer) +{ + const int3 cudaVolumeSize = make_int3(volumeSize[0], volumeSize[1], volumeSize[2]); + if (!*workspacePointer) + *workspacePointer = new CudaZengWorkspace; + auto & workspace = *static_cast(*workspacePointer); + const int pixelCount = rotatedSize[0] * rotatedSize[1]; + ensureSliceBuffers(workspace, pixelCount); + + const auto key = makeMetadataKey(0, + projectionSize, + volumeSize, + rotatedSize, + rotatedSpacing, + sigmaZero, + alpha, + rotatedToVolumeMatrices, + farDistances, + nullptr); + auto [metadataIterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = metadataIterator->second; + if (inserted) + { + GaussianMetadata metadata; + std::map> coefficientCache; + for (int projection = 0; projection < projectionSize[2]; ++projection) + { + float distance = farDistances[projection]; + for (int z = rotatedSize[2] - 2; z >= 0 && distance - rotatedSpacing[2] >= 0.f; --z) + { + const float variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + + 2.f * rotatedSpacing[2] * alpha * sigmaZero - + alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; + appendGaussianMetadata( + metadata, coefficientCache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); + distance -= rotatedSpacing[2]; + } + const float finalVariance = (alpha * distance + sigmaZero) * (alpha * distance + sigmaZero); + appendGaussianMetadata(metadata, coefficientCache, finalVariance, rotatedSpacing[0], rotatedSpacing[1]); + } + deviceMetadata.radiiX = std::move(metadata.radiiX); + deviceMetadata.radiiY = std::move(metadata.radiiY); + cudaMalloc(&deviceMetadata.coefficientsX, metadata.coefficientsX.size() * sizeof(float)); + cudaMalloc(&deviceMetadata.coefficientsY, metadata.coefficientsY.size() * sizeof(float)); + cudaMalloc(&deviceMetadata.matrices, 12 * projectionSize[2] * sizeof(float)); + cudaMemcpy(deviceMetadata.coefficientsX, + metadata.coefficientsX.data(), + metadata.coefficientsX.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy(deviceMetadata.coefficientsY, + metadata.coefficientsY.data(), + metadata.coefficientsY.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy( + deviceMetadata.matrices, rotatedToVolumeMatrices, 12 * projectionSize[2] * sizeof(float), cudaMemcpyHostToDevice); + } + + float * current = workspace.current; + float * blurred = workspace.blurred; + float * gaussianScratch = workspace.gaussianScratch; + const dim3 block(16, 16); + const dim3 grid(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16)); + + size_t gaussianIndex = 0; + for (int projection = 0; projection < projectionSize[2]; ++projection) + { + const float * matrix = deviceMetadata.matrices + 12 * projection; + float distance = farDistances[projection]; + int z = rotatedSize[2] - 1; + sampleForwardSlice<<>>(current, + nullptr, + devVolume, + devAttenuation, + cudaVolumeSize, + matrix, + rotatedSize[0], + rotatedSize[1], + z, + rotatedSpacing[2], + false); + for (--z; z >= 0 && distance - rotatedSpacing[2] >= 0.f; --z) + { + gaussian2D(current, + blurred, + gaussianScratch, + deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.radiiX[gaussianIndex], + deviceMetadata.radiiY[gaussianIndex], + rotatedSize[0], + rotatedSize[1], + grid, + block); + ++gaussianIndex; + sampleForwardSlice<<>>(current, + blurred, + devVolume, + devAttenuation, + cudaVolumeSize, + matrix, + rotatedSize[0], + rotatedSize[1], + z, + rotatedSpacing[2], + true); + distance -= rotatedSpacing[2]; + } + gaussian2D(current, + blurred, + gaussianScratch, + deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.radiiX[gaussianIndex], + deviceMetadata.radiiY[gaussianIndex], + rotatedSize[0], + rotatedSize[1], + grid, + block); + ++gaussianIndex; + finishForward<<>>(devProjectionIn + projection * pixelCount, + blurred, + devProjectionOut + projection * pixelCount, + pixelCount, + rotatedSpacing[2]); + } + CUDA_CHECK_ERROR; +} + +void +CUDA_zeng_back_project(const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + const float * volumeToRotatedMatrices, + const float * nearDistances, + const int * firstSlices, + const float * devVolumeIn, + float * devVolumeOut, + const float * devProjections, + const float * devAttenuation, + float sigmaZero, + float alpha, + void ** workspacePointer) +{ + const size_t volumeBytes = static_cast(volumeSize[0]) * volumeSize[1] * volumeSize[2] * sizeof(float); + if (!*workspacePointer) + *workspacePointer = new CudaZengWorkspace; + auto & workspace = *static_cast(*workspacePointer); + const int slicePixels = rotatedSize[0] * rotatedSize[1]; + const size_t rotatedElements = static_cast(slicePixels) * rotatedSize[2]; + const size_t rotatedBytes = rotatedElements * sizeof(float); + ensureSliceBuffers(workspace, slicePixels); + ensureRotatedBuffer(workspace, rotatedElements); + + const auto key = makeMetadataKey(1, + projectionSize, + volumeSize, + rotatedSize, + rotatedSpacing, + sigmaZero, + alpha, + volumeToRotatedMatrices, + nearDistances, + firstSlices); + auto [metadataIterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = metadataIterator->second; + if (inserted) + { + GaussianMetadata metadata; + std::map> coefficientCache; + std::vector inverseMatrices(12 * projectionSize[2]); + for (int projection = 0; projection < projectionSize[2]; ++projection) + { + float distance = nearDistances[projection]; + float variance = (alpha * distance + sigmaZero) * (alpha * distance + sigmaZero); + appendGaussianMetadata(metadata, coefficientCache, variance, rotatedSpacing[0], rotatedSpacing[1]); + for (int z = firstSlices[projection]; z + 1 < rotatedSize[2]; ++z) + { + distance += rotatedSpacing[2]; + variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + 2.f * rotatedSpacing[2] * alpha * sigmaZero - + alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; + appendGaussianMetadata( + metadata, coefficientCache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); + } + invertAffine(volumeToRotatedMatrices + 12 * projection, inverseMatrices.data() + 12 * projection); + } + deviceMetadata.radiiX = std::move(metadata.radiiX); + deviceMetadata.radiiY = std::move(metadata.radiiY); + cudaMalloc(&deviceMetadata.coefficientsX, metadata.coefficientsX.size() * sizeof(float)); + cudaMalloc(&deviceMetadata.coefficientsY, metadata.coefficientsY.size() * sizeof(float)); + cudaMalloc(&deviceMetadata.matrices, 12 * projectionSize[2] * sizeof(float)); + cudaMalloc(&deviceMetadata.inverseMatrices, 12 * projectionSize[2] * sizeof(float)); + cudaMemcpy(deviceMetadata.coefficientsX, + metadata.coefficientsX.data(), + metadata.coefficientsX.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy(deviceMetadata.coefficientsY, + metadata.coefficientsY.data(), + metadata.coefficientsY.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy( + deviceMetadata.matrices, volumeToRotatedMatrices, 12 * projectionSize[2] * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(deviceMetadata.inverseMatrices, + inverseMatrices.data(), + inverseMatrices.size() * sizeof(float), + cudaMemcpyHostToDevice); + } + + float * rotated = workspace.rotated; + float * current = workspace.current; + float * blurred = workspace.blurred; + float * gaussianScratch = workspace.gaussianScratch; + cudaMemcpy(devVolumeOut, devVolumeIn, volumeBytes, cudaMemcpyDeviceToDevice); + + cudaArray * attenuationArray = nullptr; + cudaTextureObject_t attenuationTexture = 0; + if (devAttenuation) + prepareScalarTextureObject(const_cast(volumeSize), + const_cast(devAttenuation), + attenuationArray, + attenuationTexture, + false, + true, + cudaAddressModeBorder); + + const dim3 block2(16, 16); + const dim3 grid2(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16)); + const dim3 block3(8, 8, 4); + const dim3 grid3(iDivUp(volumeSize[0], 8), iDivUp(volumeSize[1], 8), iDivUp(volumeSize[2], 4)); + const int3 cudaVolumeSize = make_int3(volumeSize[0], volumeSize[1], volumeSize[2]); + const int3 cudaRotatedSize = make_int3(rotatedSize[0], rotatedSize[1], rotatedSize[2]); + + size_t gaussianIndex = 0; + for (int projection = 0; projection < projectionSize[2]; ++projection) + { + const float * matrix = deviceMetadata.matrices + 12 * projection; + const float * inverseMatrix = deviceMetadata.inverseMatrices + 12 * projection; + cudaMemset(rotated, 0, rotatedBytes); + copyProjection<<>>(devProjections, current, slicePixels, projection); + float distance = nearDistances[projection]; + gaussian2D(current, + blurred, + gaussianScratch, + deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.radiiX[gaussianIndex], + deviceMetadata.radiiY[gaussianIndex], + rotatedSize[0], + rotatedSize[1], + grid2, + block2); + ++gaussianIndex; + std::swap(current, blurred); + + for (int z = firstSlices[projection]; z < rotatedSize[2]; ++z) + { + storeSlice<<>>(current, rotated, slicePixels, z); + if (z + 1 == rotatedSize[2]) + break; + if (attenuationTexture) + attenuateSlice<<>>( + current, attenuationTexture, inverseMatrix, rotatedSize[0], rotatedSize[1], z + 1, rotatedSpacing[2]); + distance += rotatedSpacing[2]; + gaussian2D(current, + blurred, + gaussianScratch, + deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, + deviceMetadata.radiiX[gaussianIndex], + deviceMetadata.radiiY[gaussianIndex], + rotatedSize[0], + rotatedSize[1], + grid2, + block2); + ++gaussianIndex; + std::swap(current, blurred); + } + addRotatedVolume<<>>( + rotated, devVolumeOut, devVolumeOut, cudaVolumeSize, cudaRotatedSize, matrix, rotatedSpacing[2]); + } + CUDA_CHECK_ERROR; + if (attenuationArray) + { + cudaDestroyTextureObject(attenuationTexture); + cudaFreeArray(attenuationArray); + } + CUDA_CHECK_ERROR; +} + +void +CUDA_zeng_release_workspace(void * workspace) +{ + delete static_cast(workspace); +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 18da9b27c..3794b4f4e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -154,6 +154,8 @@ rtk_add_test(rtkForwardProjectionTest rtkforwardprojectiontest.cxx) rtk_add_cuda_test(rtkForwardProjectionCudaTest rtkforwardprojectiontest.cxx) rtk_add_test(rtkForwardAttenuatedProjectionTest rtkforwardattenuatedprojectiontest.cxx) rtk_add_test(rtkZengProjectionTest rtkzengforwardprojectiontest.cxx) +rtk_add_cuda_test(rtkZengProjectionCudaTest rtkzengforwardprojectiontest.cxx) +rtk_add_cuda_test(rtkZengProjectionComparisonCudaTest rtkzengprojectioncomparisoncudatest.cxx) rtk_add_test(rtkMaximumIntensityProjectionTest rtkmaximumintensityprojectiontest.cxx) rtk_add_test(rtkGeometryFileTest rtkgeometryfiletest.cxx) diff --git a/test/rtkadjointoperatorstest.cxx b/test/rtkadjointoperatorstest.cxx index ce1af3eb0..8152047cb 100644 --- a/test/rtkadjointoperatorstest.cxx +++ b/test/rtkadjointoperatorstest.cxx @@ -11,6 +11,8 @@ #ifdef USE_CUDA # include "rtkCudaForwardProjectionImageFilter.h" # include "rtkCudaRayCastBackProjectionImageFilter.h" +# include "rtkCudaZengBackProjectionImageFilter.h" +# include "rtkCudaZengForwardProjectionImageFilter.h" #endif /** @@ -20,8 +22,8 @@ * * This test generates a random volume "v" and a random set of projections "p", * and compares the scalar products and , where R is either the - * Joseph forward projector or the Cuda ray cast forward projector, - * and R* is either the Joseph back projector or the Cuda ray cast back projector. + * Joseph or Zeng forward projector, including their CUDA implementations, + * and R* is the corresponding backprojector. * If R* is indeed the adjoint of R, these scalar products are equal. * * \author Cyril Mory @@ -299,48 +301,80 @@ rtkadjointoperatorstest(int, char *[]) std::cout << "\n\n****** Zeng Forward projector ******" << std::endl; +#ifdef USE_CUDA + using ZengForwardProjectorType = rtk::CudaZengForwardProjectionImageFilter; +#else using ZengForwardProjectorType = rtk::ZengForwardProjectionImageFilter; +#endif auto zfw = ZengForwardProjectorType::New(); zfw->SetInput(0, constantProjectionsSource->GetOutput()); zfw->SetInput(1, randomVolumeSource->GetOutput()); zfw->SetGeometry(geometry_parallel); + zfw->SetSigmaZero(1.5); + zfw->SetAlpha(0.016); TRY_AND_EXIT_ON_ITK_EXCEPTION(zfw->Update()); std::cout << "\n\n****** Zeng Back projector ******" << std::endl; +#ifdef USE_CUDA + using ZengBackProjectorType = rtk::CudaZengBackProjectionImageFilter; +#else using ZengBackProjectorType = rtk::ZengBackProjectionImageFilter; +#endif auto zbp = ZengBackProjectorType::New(); zbp->SetInput(0, constantVolumeSource->GetOutput()); zbp->SetInput(1, randomProjectionsSource->GetOutput()); zbp->SetGeometry(geometry_parallel); + zbp->SetSigmaZero(1.5); + zbp->SetAlpha(0.016); TRY_AND_EXIT_ON_ITK_EXCEPTION(zbp->Update()); +#ifdef USE_CUDA + if (!zfw->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty() || + !zbp->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "CUDA Zeng output was unexpectedly synchronized back to the CPU." << std::endl; + return EXIT_FAILURE; + } +#endif + CheckScalarProducts( randomVolumeSource->GetOutput(), zbp->GetOutput(), randomProjectionsSource->GetOutput(), zfw->GetOutput()); std::cout << "\n\nTest PASSED! " << std::endl; std::cout << "\n\n****** Attenuated Zeng Forward projector ******" << std::endl; - using ZengForwardProjectorType = rtk::ZengForwardProjectionImageFilter; auto attzfw = ZengForwardProjectorType::New(); attzfw->SetInput(0, constantProjectionsSource->GetOutput()); attzfw->SetInput(1, randomVolumeSource->GetOutput()); attzfw->SetInput(2, constantAttenuationSource->GetOutput()); attzfw->SetGeometry(geometry_parallel); + attzfw->SetSigmaZero(1.5); + attzfw->SetAlpha(0.016); TRY_AND_EXIT_ON_ITK_EXCEPTION(attzfw->Update()); std::cout << "\n\n****** Attenuated Zeng Back projector ******" << std::endl; - using ZengBackProjectorType = rtk::ZengBackProjectionImageFilter; auto attzbp = ZengBackProjectorType::New(); attzbp->SetInput(0, constantVolumeSource->GetOutput()); attzbp->SetInput(1, randomProjectionsSource->GetOutput()); attzbp->SetInput(2, constantAttenuationSource->GetOutput()); attzbp->SetGeometry(geometry_parallel); + attzbp->SetSigmaZero(1.5); + attzbp->SetAlpha(0.016); TRY_AND_EXIT_ON_ITK_EXCEPTION(attzbp->Update()); +#ifdef USE_CUDA + if (!attzfw->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty() || + !attzbp->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "Attenuated CUDA Zeng output was unexpectedly synchronized back to the CPU." << std::endl; + return EXIT_FAILURE; + } +#endif + CheckScalarProducts( randomVolumeSource->GetOutput(), attzbp->GetOutput(), randomProjectionsSource->GetOutput(), attzfw->GetOutput()); std::cout << "\n\nTest PASSED! " << std::endl; diff --git a/test/rtkzengforwardprojectiontest.cxx b/test/rtkzengforwardprojectiontest.cxx index 8b4d46992..43b86bda4 100644 --- a/test/rtkzengforwardprojectiontest.cxx +++ b/test/rtkzengforwardprojectiontest.cxx @@ -4,6 +4,9 @@ #include "rtkJosephForwardAttenuatedProjectionImageFilter.h" #include "rtkThreeDCircularProjectionGeometryXMLFile.h" #include "rtkZengForwardProjectionImageFilter.h" +#ifdef USE_CUDA +# include "rtkCudaZengForwardProjectionImageFilter.h" +#endif #include #include #include @@ -85,7 +88,11 @@ rtkzengforwardprojectiontest(int, char *[]) // Zeng Forward Projection filter +#ifdef USE_CUDA + auto jfp = rtk::CudaZengForwardProjectionImageFilter::New(); +#else auto jfp = rtk::ZengForwardProjectionImageFilter::New(); +#endif jfp->InPlaceOff(); jfp->SetInput(projInput->GetOutput()); jfp->SetInput(1, volInput->GetOutput()); @@ -112,6 +119,14 @@ rtkzengforwardprojectiontest(int, char *[]) jfp->SetSigmaZero(0.); jfp->Update(); +#ifdef USE_CUDA + if (!jfp->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "CUDA Zeng output was unexpectedly synchronized back to the CPU." << std::endl; + return EXIT_FAILURE; + } +#endif + CheckImageQuality(jfp->GetOutput(), attjfp->GetOutput(), 0.1, 44.0, 255.0); std::cout << "\n\nTest PASSED! " << std::endl; diff --git a/test/rtkzengprojectioncomparisoncudatest.cxx b/test/rtkzengprojectioncomparisoncudatest.cxx new file mode 100644 index 000000000..592d5c4d3 --- /dev/null +++ b/test/rtkzengprojectioncomparisoncudatest.cxx @@ -0,0 +1,202 @@ +#include "rtkConstantImageSource.h" +#include "rtkCudaZengBackProjectionImageFilter.h" +#include "rtkCudaZengForwardProjectionImageFilter.h" +#include "rtkDrawEllipsoidImageFilter.h" +#include "rtkMacro.h" +#include "rtkThreeDCircularProjectionGeometry.h" +#include "rtkZengBackProjectionImageFilter.h" +#include "rtkZengForwardProjectionImageFilter.h" + +#include +#include +#include +#include +#include +#include + +/** + * \file rtkzengprojectioncomparisoncudatest.cxx + * + * \brief Compares the CPU and CUDA Zeng projectors. + * + * The test checks forward projection and backprojection, with and without an + * attenuation map. It also verifies that CUDA outputs remain GPU-resident. + */ + +namespace +{ +template +double +RelativeL2Error(const TReferenceImage * reference, const TResultImage * result) +{ + itk::ImageRegionConstIterator referenceIterator(reference, reference->GetBufferedRegion()); + itk::ImageRegionConstIterator resultIterator(result, result->GetBufferedRegion()); + double squaredError = 0.; + double squaredReference = 0.; + for (referenceIterator.GoToBegin(), resultIterator.GoToBegin(); !referenceIterator.IsAtEnd(); + ++referenceIterator, ++resultIterator) + { + const double referenceValue = referenceIterator.Get(); + const double resultValue = resultIterator.Get(); + if (!std::isfinite(resultValue)) + { + std::cerr << "CUDA Zeng produced a non-finite value." << std::endl; + return std::numeric_limits::infinity(); + } + const double difference = resultValue - referenceValue; + squaredError += difference * difference; + squaredReference += referenceValue * referenceValue; + } + return std::sqrt(squaredError / squaredReference); +} + +void +CheckError(const char * label, double value, double tolerance) +{ + std::cout << label << " relative L2 error: " << value << std::endl; + if (!(value < tolerance)) + { + std::cerr << label << " relative L2 error exceeds " << tolerance << std::endl; + std::exit(EXIT_FAILURE); + } +} +} // namespace + +int +rtkzengprojectioncomparisoncudatest(int, char *[]) +{ + constexpr unsigned int Dimension = 3; + using CPUImageType = itk::Image; + using CudaImageType = itk::CudaImage; + using CPUConstantSourceType = rtk::ConstantImageSource; + using CudaConstantSourceType = rtk::ConstantImageSource; + + auto cpuVolumeSource = CPUConstantSourceType::New(); + auto cudaVolumeSource = CudaConstantSourceType::New(); + cpuVolumeSource->SetOrigin(itk::MakePoint(-30., -30., -30.)); + cpuVolumeSource->SetSpacing(itk::MakeVector(4., 4., 4.)); + cpuVolumeSource->SetSize(itk::MakeSize(16, 16, 16)); + cpuVolumeSource->SetConstant(0.); + cudaVolumeSource->SetOrigin(itk::MakePoint(-30., -30., -30.)); + cudaVolumeSource->SetSpacing(itk::MakeVector(4., 4., 4.)); + cudaVolumeSource->SetSize(itk::MakeSize(16, 16, 16)); + cudaVolumeSource->SetConstant(0.); + + auto cpuVolume = rtk::DrawEllipsoidImageFilter::New(); + cpuVolume->SetInput(cpuVolumeSource->GetOutput()); + cpuVolume->SetCenter(itk::MakePoint(0., 0., 0.)); + cpuVolume->SetAxis(itk::MakeVector(20., 16., 12.)); + cpuVolume->SetDensity(1.); + auto cudaVolume = rtk::DrawEllipsoidImageFilter::New(); + cudaVolume->SetInput(cudaVolumeSource->GetOutput()); + cudaVolume->SetCenter(itk::MakePoint(0., 0., 0.)); + cudaVolume->SetAxis(itk::MakeVector(20., 16., 12.)); + cudaVolume->SetDensity(1.); + + auto cpuProjections = CPUConstantSourceType::New(); + cpuProjections->SetOrigin(itk::MakePoint(-30., -30., 0.)); + cpuProjections->SetSpacing(itk::MakeVector(4., 4., 1.)); + cpuProjections->SetSize(itk::MakeSize(16, 16, 4)); + cpuProjections->SetConstant(0.); + auto cudaProjections = CudaConstantSourceType::New(); + cudaProjections->SetOrigin(itk::MakePoint(-30., -30., 0.)); + cudaProjections->SetSpacing(itk::MakeVector(4., 4., 1.)); + cudaProjections->SetSize(itk::MakeSize(16, 16, 4)); + cudaProjections->SetConstant(0.); + + auto cpuAttenuation = CPUConstantSourceType::New(); + cpuAttenuation->SetOrigin(itk::MakePoint(-30., -30., -30.)); + cpuAttenuation->SetSpacing(itk::MakeVector(4., 4., 4.)); + cpuAttenuation->SetSize(itk::MakeSize(16, 16, 16)); + cpuAttenuation->SetConstant(0.01); + auto cudaAttenuation = CudaConstantSourceType::New(); + cudaAttenuation->SetOrigin(itk::MakePoint(-30., -30., -30.)); + cudaAttenuation->SetSpacing(itk::MakeVector(4., 4., 4.)); + cudaAttenuation->SetSize(itk::MakeSize(16, 16, 16)); + cudaAttenuation->SetConstant(0.01); + + auto geometry = rtk::ThreeDCircularProjectionGeometry::New(); + for (unsigned int projection = 0; projection < 4; ++projection) + geometry->AddProjection(100., 0., projection * 90.); + + using CPUForwardType = rtk::ZengForwardProjectionImageFilter; + auto cpuForward = CPUForwardType::New(); + cpuForward->InPlaceOff(); + cpuForward->SetInput(0, cpuProjections->GetOutput()); + cpuForward->SetInput(1, cpuVolume->GetOutput()); + cpuForward->SetGeometry(geometry); + cpuForward->SetSigmaZero(1.5); + cpuForward->SetAlpha(0.016); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cpuForward->Update()); + + auto cudaForward = rtk::CudaZengForwardProjectionImageFilter::New(); + cudaForward->InPlaceOff(); + cudaForward->SetInput(0, cudaProjections->GetOutput()); + cudaForward->SetInput(1, cudaVolume->GetOutput()); + cudaForward->SetGeometry(geometry); + cudaForward->SetSigmaZero(1.5); + cudaForward->SetAlpha(0.016); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cudaForward->Update()); + if (!cudaForward->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "CUDA Zeng forward output was synchronized to the CPU." << std::endl; + return EXIT_FAILURE; + } + CheckError("Forward", RelativeL2Error(cpuForward->GetOutput(), cudaForward->GetOutput()), 5.e-5); + + using CPUBackType = rtk::ZengBackProjectionImageFilter; + auto cpuBack = CPUBackType::New(); + cpuBack->InPlaceOff(); + cpuBack->SetInput(0, cpuVolumeSource->GetOutput()); + cpuBack->SetInput(1, cpuForward->GetOutput()); + cpuBack->SetGeometry(geometry); + cpuBack->SetSigmaZero(1.5); + cpuBack->SetAlpha(0.016); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cpuBack->Update()); + + auto cudaBack = rtk::CudaZengBackProjectionImageFilter::New(); + cudaBack->InPlaceOff(); + cudaBack->SetInput(0, cudaVolumeSource->GetOutput()); + cudaBack->SetInput(1, cudaForward->GetOutput()); + cudaBack->SetGeometry(geometry); + cudaBack->SetSigmaZero(1.5); + cudaBack->SetAlpha(0.016); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cudaBack->Update()); + if (!cudaBack->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "CUDA Zeng backprojection output was synchronized to the CPU." << std::endl; + return EXIT_FAILURE; + } + CheckError("Backprojection", RelativeL2Error(cpuBack->GetOutput(), cudaBack->GetOutput()), 5.e-5); + + cpuForward->SetInput(2, cpuAttenuation->GetOutput()); + cpuForward->Modified(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cpuForward->Update()); + cudaForward->SetInput(2, cudaAttenuation->GetOutput()); + cudaForward->Modified(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cudaForward->Update()); + if (!cudaForward->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "Attenuated CUDA Zeng forward output was synchronized to the CPU." << std::endl; + return EXIT_FAILURE; + } + CheckError("Attenuated forward", RelativeL2Error(cpuForward->GetOutput(), cudaForward->GetOutput()), 5.e-4); + + cpuBack->SetInput(1, cpuForward->GetOutput()); + cpuBack->SetInput(2, cpuAttenuation->GetOutput()); + cpuBack->Modified(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cpuBack->Update()); + cudaBack->SetInput(1, cudaForward->GetOutput()); + cudaBack->SetInput(2, cudaAttenuation->GetOutput()); + cudaBack->Modified(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(cudaBack->Update()); + if (!cudaBack->GetOutput()->GetCudaDataManager()->IsCPUBufferDirty()) + { + std::cerr << "Attenuated CUDA Zeng backprojection output was synchronized to the CPU." << std::endl; + return EXIT_FAILURE; + } + CheckError("Attenuated backprojection", RelativeL2Error(cpuBack->GetOutput(), cudaBack->GetOutput()), 2.e-3); + + std::cout << "CUDA Zeng CPU comparison test PASSED." << std::endl; + return EXIT_SUCCESS; +} diff --git a/wrapping/rtkCudaZengBackProjectionImageFilter.wrap b/wrapping/rtkCudaZengBackProjectionImageFilter.wrap new file mode 100644 index 000000000..a56c77121 --- /dev/null +++ b/wrapping/rtkCudaZengBackProjectionImageFilter.wrap @@ -0,0 +1,3 @@ +if(RTK_USE_CUDA) + itk_wrap_simple_class("rtk::CudaZengBackProjectionImageFilter" POINTER) +endif() diff --git a/wrapping/rtkCudaZengForwardProjectionImageFilter.wrap b/wrapping/rtkCudaZengForwardProjectionImageFilter.wrap new file mode 100644 index 000000000..af2e43e79 --- /dev/null +++ b/wrapping/rtkCudaZengForwardProjectionImageFilter.wrap @@ -0,0 +1,3 @@ +if(RTK_USE_CUDA) + itk_wrap_simple_class("rtk::CudaZengForwardProjectionImageFilter" POINTER) +endif() From c9137987e8a018dd4a3a55186c535251497a2dff Mon Sep 17 00:00:00 2001 From: "C. Constanza" Date: Wed, 16 Sep 2026 20:41:25 +0200 Subject: [PATCH 2/4] ENH: Address CUDA Zeng review comments --- .../rtkbackprojections/rtkbackprojections.py | 27 +++++------ .../rtkforwardprojections.py | 16 +++---- applications/rtkprojectors_group.py | 40 +++++++--------- .../rtkCudaZengBackProjectionImageFilter.h | 20 ++++---- .../rtkCudaZengForwardProjectionImageFilter.h | 20 ++++---- src/rtkCudaZengBackProjectionImageFilter.cxx | 46 +++++++----------- ...tkCudaZengForwardProjectionImageFilter.cxx | 48 +++++++------------ src/rtkCudaZengProjectionImageFilter.cu | 3 +- 8 files changed, 90 insertions(+), 130 deletions(-) diff --git a/applications/rtkbackprojections/rtkbackprojections.py b/applications/rtkbackprojections/rtkbackprojections.py index 11d1d26f0..4e907dbad 100644 --- a/applications/rtkbackprojections/rtkbackprojections.py +++ b/applications/rtkbackprojections/rtkbackprojections.py @@ -131,23 +131,20 @@ 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: - bp.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf: - bp.SetAlpha(args_info.alphapsf) - - elif args_info.bp == "CudaZeng": - if hasattr(itk, "CudaImage"): + 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() - if args_info.sigmazero: - bp.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf: - bp.SetAlpha(args_info.alphapsf) else: - print("The program has not been compiled with cuda option") - sys.exit(1) + bp = rtk.ZengBackProjectionImageFilter[ + OutputImageType, OutputImageType + ].New() + if args_info.sigmazero is not None: + bp.SetSigmaZero(args_info.sigmazero) + if args_info.alphapsf is not None: + bp.SetAlpha(args_info.alphapsf) elif args_info.bp == "CudaFDKBackProjection": if hasattr(itk, "CudaImage"): diff --git a/applications/rtkforwardprojections/rtkforwardprojections.py b/applications/rtkforwardprojections/rtkforwardprojections.py index 789de1664..e1f3bd043 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.py +++ b/applications/rtkforwardprojections/rtkforwardprojections.py @@ -158,16 +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 == "CudaZeng": - if hasattr(itk, "CudaImage"): + 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: - print("The program has not been compiled with cuda option") - sys.exit(1) + forwardProjection = rtk.ZengForwardProjectionImageFilter[ + OutputImageType, OutputImageType + ].New() elif args_info.fp == "MIP": forwardProjection = rtk.MaximumIntensityProjectionImageFilter[ OutputImageType, OutputImageType diff --git a/applications/rtkprojectors_group.py b/applications/rtkprojectors_group.py index 85ba0d887..72657f85b 100644 --- a/applications/rtkprojectors_group.py +++ b/applications/rtkprojectors_group.py @@ -38,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( @@ -90,16 +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) - if args_info.sigmazero is not None: - recon.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf is not None: - recon.SetAlphaPSF(args_info.alphapsf) - if args_info.attenuationmap is not None: - recon.SetAttenuationMap(attenuation_map) - elif args_info.bp == "CudaZeng": - recon.SetBackProjectionFilter(ReconType.BackProjectionType_BP_CUDAZENG) + 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: @@ -128,16 +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) - if args_info.sigmazero is not None: - recon.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf is not None: - recon.SetAlphaPSF(args_info.alphapsf) - if args_info.attenuationmap is not None: - recon.SetAttenuationMap(attenuation_map) - elif args_info.fp == "CudaZeng": - recon.SetForwardProjectionFilter(ReconType.ForwardProjectionType_FP_CUDAZENG) + 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: diff --git a/include/rtkCudaZengBackProjectionImageFilter.h b/include/rtkCudaZengBackProjectionImageFilter.h index 6de7519f6..b210bb2d2 100644 --- a/include/rtkCudaZengBackProjectionImageFilter.h +++ b/include/rtkCudaZengBackProjectionImageFilter.h @@ -21,7 +21,7 @@ #include "rtkConfiguration.h" #ifdef RTK_USE_CUDA -# include "rtkBackProjectionImageFilter.h" +# include "rtkZengBackProjectionImageFilter.h" # include "RTKExport.h" # include # include @@ -36,13 +36,16 @@ namespace rtk * 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, - BackProjectionImageFilter, itk::CudaImage>> + : public itk::CudaInPlaceImageFilter< + itk::CudaImage, + itk::CudaImage, + ZengBackProjectionImageFilter, itk::CudaImage>> { public: ITK_DISALLOW_COPY_AND_MOVE(CudaZengBackProjectionImageFilter); @@ -50,7 +53,7 @@ class RTK_EXPORT CudaZengBackProjectionImageFilter /** 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 = BackProjectionImageFilter; + using ProjectorType = ZengBackProjectionImageFilter; using Self = CudaZengBackProjectionImageFilter; using Superclass = itk::CudaInPlaceImageFilter; using Pointer = itk::SmartPointer; @@ -59,11 +62,6 @@ class RTK_EXPORT CudaZengBackProjectionImageFilter itkNewMacro(Self); itkOverrideGetNameOfClassMacro(CudaZengBackProjectionImageFilter); - itkGetConstMacro(SigmaZero, double); - itkSetMacro(SigmaZero, double); - itkGetConstMacro(Alpha, double); - itkSetMacro(Alpha, double); - protected: CudaZengBackProjectionImageFilter(); ~CudaZengBackProjectionImageFilter() override; @@ -71,8 +69,6 @@ class RTK_EXPORT CudaZengBackProjectionImageFilter GPUGenerateData() override; private: - double m_SigmaZero{ 1.5417233052142099 }; - double m_Alpha{ 0.016241189545787734 }; void * m_CudaWorkspace{ nullptr }; }; diff --git a/include/rtkCudaZengForwardProjectionImageFilter.h b/include/rtkCudaZengForwardProjectionImageFilter.h index 8b4d2fa15..04dd0ce60 100644 --- a/include/rtkCudaZengForwardProjectionImageFilter.h +++ b/include/rtkCudaZengForwardProjectionImageFilter.h @@ -21,7 +21,7 @@ #include "rtkConfiguration.h" #ifdef RTK_USE_CUDA -# include "rtkForwardProjectionImageFilter.h" +# include "rtkZengForwardProjectionImageFilter.h" # include "RTKExport.h" # include # include @@ -36,13 +36,16 @@ namespace rtk * The implementation reproduces the slice recursion of * ZengForwardProjectionImageFilter, 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 CudaZengForwardProjectionImageFilter - : public itk::CudaInPlaceImageFilter, - itk::CudaImage, - ForwardProjectionImageFilter, itk::CudaImage>> + : public itk::CudaInPlaceImageFilter< + itk::CudaImage, + itk::CudaImage, + ZengForwardProjectionImageFilter, itk::CudaImage>> { public: ITK_DISALLOW_COPY_AND_MOVE(CudaZengForwardProjectionImageFilter); @@ -50,7 +53,7 @@ class RTK_EXPORT CudaZengForwardProjectionImageFilter /** 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 = ForwardProjectionImageFilter; + using ProjectorType = ZengForwardProjectionImageFilter; using Self = CudaZengForwardProjectionImageFilter; using Superclass = itk::CudaInPlaceImageFilter; using Pointer = itk::SmartPointer; @@ -59,11 +62,6 @@ class RTK_EXPORT CudaZengForwardProjectionImageFilter itkNewMacro(Self); itkOverrideGetNameOfClassMacro(CudaZengForwardProjectionImageFilter); - itkGetConstMacro(SigmaZero, double); - itkSetMacro(SigmaZero, double); - itkGetConstMacro(Alpha, double); - itkSetMacro(Alpha, double); - protected: CudaZengForwardProjectionImageFilter(); ~CudaZengForwardProjectionImageFilter() override; @@ -71,8 +69,6 @@ class RTK_EXPORT CudaZengForwardProjectionImageFilter GPUGenerateData() override; private: - double m_SigmaZero{ 1.5417233052142099 }; - double m_Alpha{ 0.016241189545787734 }; void * m_CudaWorkspace{ nullptr }; }; diff --git a/src/rtkCudaZengBackProjectionImageFilter.cxx b/src/rtkCudaZengBackProjectionImageFilter.cxx index 40a7cbcf6..d95bd65f7 100644 --- a/src/rtkCudaZengBackProjectionImageFilter.cxx +++ b/src/rtkCudaZengBackProjectionImageFilter.cxx @@ -27,30 +27,9 @@ namespace rtk namespace { using MatrixType = ThreeDCircularProjectionGeometry::ThreeDHomogeneousMatrixType; - -MatrixType -TransformMatrix(const itk::CenteredEuler3DTransform * transform) -{ - MatrixType matrix; - matrix.SetIdentity(); - for (unsigned int r = 0; r < 3; ++r) - { - for (unsigned int c = 0; c < 3; ++c) - matrix[r][c] = transform->GetMatrix()[r][c]; - matrix[r][3] = transform->GetOffset()[r]; - } - return matrix; -} -void -StoreMatrix(const MatrixType & matrix, float * destination) -{ - for (unsigned int r = 0; r < 3; ++r) - for (unsigned int c = 0; c < 4; ++c) - destination[4 * r + c] = static_cast(matrix[r][c]); -} } // namespace -CudaZengBackProjectionImageFilter::CudaZengBackProjectionImageFilter() { this->InPlaceOff(); } +CudaZengBackProjectionImageFilter::CudaZengBackProjectionImageFilter() = default; CudaZengBackProjectionImageFilter::~CudaZengBackProjectionImageFilter() { @@ -133,11 +112,20 @@ CudaZengBackProjectionImageFilter::GPUGenerateData() auto inverse = TransformType::New(); if (!transform->GetInverse(inverse)) itkGenericExceptionMacro(<< "Could not invert Zeng rotation transform."); - const MatrixType matrix = GetPhysicalPointToIndexMatrix(rotatedImage.GetPointer()).GetVnlMatrix() * - TransformMatrix(inverse).GetVnlMatrix() * - GetIndexToPhysicalPointMatrix(this->GetInput(0)).GetVnlMatrix() * - volumeIndexTranslation.GetVnlMatrix(); - StoreMatrix(matrix, matrices.data() + 12 * local); + MatrixType inverseMatrix; + inverseMatrix.SetIdentity(); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 3; ++c) + inverseMatrix[r][c] = inverse->GetMatrix()[r][c]; + inverseMatrix[r][3] = inverse->GetOffset()[r]; + } + const MatrixType matrix = + GetPhysicalPointToIndexMatrix(rotatedImage.GetPointer()).GetVnlMatrix() * inverseMatrix.GetVnlMatrix() * + GetIndexToPhysicalPointMatrix(this->GetInput(0)).GetVnlMatrix() * volumeIndexTranslation.GetVnlMatrix(); + for (unsigned int r = 0; r < 3; ++r) + for (unsigned int c = 0; c < 4; ++c) + matrices[12 * local + 4 * r + c] = static_cast(matrix[r][c]); const double firstDistance = geometry->GetSourceToIsocenterDistances()[projection] + origin[2]; firstSlices[local] = std::max(0, static_cast(std::ceil(-firstDistance / spacing[2]))); @@ -163,8 +151,8 @@ CudaZengBackProjectionImageFilter::GPUGenerateData() volumeOut, projections, attenuation, - static_cast(m_SigmaZero), - static_cast(m_Alpha), + static_cast(this->GetSigmaZero()), + static_cast(this->GetAlpha()), &m_CudaWorkspace); } diff --git a/src/rtkCudaZengForwardProjectionImageFilter.cxx b/src/rtkCudaZengForwardProjectionImageFilter.cxx index 9d763f1cc..2a0cb4d2a 100644 --- a/src/rtkCudaZengForwardProjectionImageFilter.cxx +++ b/src/rtkCudaZengForwardProjectionImageFilter.cxx @@ -27,34 +27,9 @@ namespace rtk namespace { using MatrixType = ThreeDCircularProjectionGeometry::ThreeDHomogeneousMatrixType; - -MatrixType -TransformMatrix(const itk::CenteredEuler3DTransform * transform) -{ - MatrixType matrix; - matrix.SetIdentity(); - for (unsigned int r = 0; r < 3; ++r) - { - for (unsigned int c = 0; c < 3; ++c) - matrix[r][c] = transform->GetMatrix()[r][c]; - matrix[r][3] = transform->GetOffset()[r]; - } - return matrix; -} - -void -StoreTextureMatrix(const MatrixType & matrix, const itk::Index<3> & bufferedIndex, float * destination) -{ - for (unsigned int r = 0; r < 3; ++r) - for (unsigned int c = 0; c < 4; ++c) - destination[4 * r + c] = static_cast(matrix[r][c]); - destination[3] -= bufferedIndex[0]; - destination[7] -= bufferedIndex[1]; - destination[11] -= bufferedIndex[2]; -} } // namespace -CudaZengForwardProjectionImageFilter::CudaZengForwardProjectionImageFilter() { this->InPlaceOff(); } +CudaZengForwardProjectionImageFilter::CudaZengForwardProjectionImageFilter() = default; CudaZengForwardProjectionImageFilter::~CudaZengForwardProjectionImageFilter() { @@ -133,10 +108,23 @@ CudaZengForwardProjectionImageFilter::GPUGenerateData() origin[2] = rotatedCenter[2] - spacing[2] * (rotatedSize[2] - 1.0) / 2.0; rotatedImage->SetOrigin(origin); + MatrixType transformMatrix; + transformMatrix.SetIdentity(); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 3; ++c) + transformMatrix[r][c] = transform->GetMatrix()[r][c]; + transformMatrix[r][3] = transform->GetOffset()[r]; + } const MatrixType matrix = GetPhysicalPointToIndexMatrix(this->GetInput(1)).GetVnlMatrix() * - TransformMatrix(transform).GetVnlMatrix() * + transformMatrix.GetVnlMatrix() * GetIndexToPhysicalPointMatrix(rotatedImage.GetPointer()).GetVnlMatrix(); - StoreTextureMatrix(matrix, volumeRegion.GetIndex(), matrices.data() + 12 * local); + for (unsigned int r = 0; r < 3; ++r) + { + for (unsigned int c = 0; c < 4; ++c) + matrices[12 * local + 4 * r + c] = static_cast(matrix[r][c]); + matrices[12 * local + 4 * r + 3] -= volumeRegion.GetIndex(r); + } farDistances[local] = static_cast(geometry->GetSourceToIsocenterDistances()[projection] + origin[2] + spacing[2] * (rotatedSize[2] - 1)); } @@ -161,8 +149,8 @@ CudaZengForwardProjectionImageFilter::GPUGenerateData() projectionOut, volume, attenuation, - static_cast(m_SigmaZero), - static_cast(m_Alpha), + static_cast(this->GetSigmaZero()), + static_cast(this->GetAlpha()), &m_CudaWorkspace); } diff --git a/src/rtkCudaZengProjectionImageFilter.cu b/src/rtkCudaZengProjectionImageFilter.cu index a09cdc8a6..1d70508aa 100644 --- a/src/rtkCudaZengProjectionImageFilter.cu +++ b/src/rtkCudaZengProjectionImageFilter.cu @@ -625,7 +625,8 @@ CUDA_zeng_back_project(const int projectionSize[3], float * current = workspace.current; float * blurred = workspace.blurred; float * gaussianScratch = workspace.gaussianScratch; - cudaMemcpy(devVolumeOut, devVolumeIn, volumeBytes, cudaMemcpyDeviceToDevice); + if (devVolumeOut != devVolumeIn) + cudaMemcpy(devVolumeOut, devVolumeIn, volumeBytes, cudaMemcpyDeviceToDevice); cudaArray * attenuationArray = nullptr; cudaTextureObject_t attenuationTexture = 0; From 66ea1a65d70b4f5d018ec003784b63ca6d80dcff Mon Sep 17 00:00:00 2001 From: "C. Constanza" Date: Thu, 17 Sep 2026 10:50:15 +0200 Subject: [PATCH 3/4] BUG: Match CPU max-slice calculation, border interpolation, and first-slice attenuation --- src/rtkCudaZengBackProjectionImageFilter.cxx | 2 +- ...tkCudaZengForwardProjectionImageFilter.cxx | 2 +- src/rtkCudaZengProjectionImageFilter.cu | 33 ++++++++++++++++--- 3 files changed, 31 insertions(+), 6 deletions(-) diff --git a/src/rtkCudaZengBackProjectionImageFilter.cxx b/src/rtkCudaZengBackProjectionImageFilter.cxx index d95bd65f7..4fdc65ceb 100644 --- a/src/rtkCudaZengBackProjectionImageFilter.cxx +++ b/src/rtkCudaZengBackProjectionImageFilter.cxx @@ -55,7 +55,7 @@ CudaZengBackProjectionImageFilter::GPUGenerateData() static_cast(volumeRegion.GetSize(2)) }; int rotatedSize[3] = { projectionSize[0], projectionSize[1], - static_cast(std::ceil(volumeRegion.GetSize(2) * std::sqrt(2.0))) }; + static_cast(volumeRegion.GetSize(2) * std::sqrt(2.0)) }; float rotatedSpacing[3] = { static_cast(this->GetInput(1)->GetSpacing()[0]), static_cast(this->GetInput(1)->GetSpacing()[1]), static_cast(this->GetInput(0)->GetSpacing()[2]) }; diff --git a/src/rtkCudaZengForwardProjectionImageFilter.cxx b/src/rtkCudaZengForwardProjectionImageFilter.cxx index 2a0cb4d2a..8e7de3cfb 100644 --- a/src/rtkCudaZengForwardProjectionImageFilter.cxx +++ b/src/rtkCudaZengForwardProjectionImageFilter.cxx @@ -57,7 +57,7 @@ CudaZengForwardProjectionImageFilter::GPUGenerateData() static_cast(volumeRegion.GetSize(2)) }; int rotatedSize[3] = { projectionSize[0], projectionSize[1], - static_cast(std::ceil(volumeRegion.GetSize(2) * std::sqrt(2.0))) }; + static_cast(volumeRegion.GetSize(2) * std::sqrt(2.0)) }; float rotatedSpacing[3] = { static_cast(this->GetInput(0)->GetSpacing()[0]), static_cast(this->GetInput(0)->GetSpacing()[1]), static_cast(this->GetInput(1)->GetSpacing()[2]) }; diff --git a/src/rtkCudaZengProjectionImageFilter.cu b/src/rtkCudaZengProjectionImageFilter.cu index 1d70508aa..57d3a714b 100644 --- a/src/rtkCudaZengProjectionImageFilter.cu +++ b/src/rtkCudaZengProjectionImageFilter.cu @@ -38,6 +38,13 @@ applyMatrix(const float * m, float x, float y, float z) __device__ float trilinearZero(const float * image, int3 size, float3 p) { + // Match ITK's linear interpolator at the half-voxel image boundary. + if (!(p.x >= -0.5f && p.x < size.x - 0.5f && p.y >= -0.5f && p.y < size.y - 0.5f && p.z >= -0.5f && + p.z < size.z - 0.5f)) + return 0.f; + p.x = fminf(fmaxf(p.x, 0.f), static_cast(size.x - 1)); + p.y = fminf(fmaxf(p.y, 0.f), static_cast(size.y - 1)); + p.z = fminf(fmaxf(p.z, 0.f), static_cast(size.z - 1)); const int x0 = static_cast(floorf(p.x)); const int y0 = static_cast(floorf(p.y)); const int z0 = static_cast(floorf(p.z)); @@ -339,6 +346,7 @@ copyProjection(const float * projections, float * slice, int count, int projecti __global__ void attenuateSlice(float * slice, cudaTextureObject_t attenuation, + int3 volumeSize, const float * rotatedToVolume, int width, int height, @@ -350,7 +358,9 @@ attenuateSlice(float * slice, if (x >= width || y >= height) return; const float3 position = applyMatrix(rotatedToVolume, x, y, z); - slice[y * width + x] *= expf(-attenuationStep * tex3D(attenuation, position.x, position.y, position.z)); + if (position.x >= 0.f && position.x < volumeSize.x && position.y >= 0.f && position.y < volumeSize.y && + position.z >= 0.f && position.z < volumeSize.z) + slice[y * width + x] *= expf(-attenuationStep * tex3D(attenuation, position.x, position.y, position.z)); } __global__ void @@ -637,7 +647,7 @@ CUDA_zeng_back_project(const int projectionSize[3], attenuationTexture, false, true, - cudaAddressModeBorder); + cudaAddressModeClamp); const dim3 block2(16, 16); const dim3 grid2(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16)); @@ -653,6 +663,15 @@ CUDA_zeng_back_project(const int projectionSize[3], const float * inverseMatrix = deviceMetadata.inverseMatrices + 12 * projection; cudaMemset(rotated, 0, rotatedBytes); copyProjection<<>>(devProjections, current, slicePixels, projection); + if (attenuationTexture) + attenuateSlice<<>>(current, + attenuationTexture, + cudaVolumeSize, + inverseMatrix, + rotatedSize[0], + rotatedSize[1], + firstSlices[projection], + rotatedSpacing[2]); float distance = nearDistances[projection]; gaussian2D(current, blurred, @@ -674,8 +693,14 @@ CUDA_zeng_back_project(const int projectionSize[3], if (z + 1 == rotatedSize[2]) break; if (attenuationTexture) - attenuateSlice<<>>( - current, attenuationTexture, inverseMatrix, rotatedSize[0], rotatedSize[1], z + 1, rotatedSpacing[2]); + attenuateSlice<<>>(current, + attenuationTexture, + cudaVolumeSize, + inverseMatrix, + rotatedSize[0], + rotatedSize[1], + z + 1, + rotatedSpacing[2]); distance += rotatedSpacing[2]; gaussian2D(current, blurred, From ee2729c752d452982c05036cd01bad30dcf38aed Mon Sep 17 00:00:00 2001 From: "C. Constanza" Date: Thu, 17 Sep 2026 10:57:47 +0200 Subject: [PATCH 4/4] ENH: Add configurable depth-slice batching to CUDA Zeng projectors --- .../rtkCudaZengBackProjectionImageFilter.h | 7 +- .../rtkCudaZengForwardProjectionImageFilter.h | 7 +- include/rtkCudaZengProjectionImageFilter.hcu | 2 + src/rtkCudaZengBackProjectionImageFilter.cxx | 1 + ...tkCudaZengForwardProjectionImageFilter.cxx | 1 + src/rtkCudaZengProjectionImageFilter.cu | 1019 ++++++++++------- 6 files changed, 620 insertions(+), 417 deletions(-) diff --git a/include/rtkCudaZengBackProjectionImageFilter.h b/include/rtkCudaZengBackProjectionImageFilter.h index b210bb2d2..2f3ee4ddf 100644 --- a/include/rtkCudaZengBackProjectionImageFilter.h +++ b/include/rtkCudaZengBackProjectionImageFilter.h @@ -62,6 +62,10 @@ class RTK_EXPORT CudaZengBackProjectionImageFilter 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; @@ -69,7 +73,8 @@ class RTK_EXPORT CudaZengBackProjectionImageFilter GPUGenerateData() override; private: - void * m_CudaWorkspace{ nullptr }; + void * m_CudaWorkspace{ nullptr }; + unsigned int m_BatchSize{ 0 }; }; } // namespace rtk diff --git a/include/rtkCudaZengForwardProjectionImageFilter.h b/include/rtkCudaZengForwardProjectionImageFilter.h index 04dd0ce60..d32c3a532 100644 --- a/include/rtkCudaZengForwardProjectionImageFilter.h +++ b/include/rtkCudaZengForwardProjectionImageFilter.h @@ -62,6 +62,10 @@ class RTK_EXPORT CudaZengForwardProjectionImageFilter itkNewMacro(Self); itkOverrideGetNameOfClassMacro(CudaZengForwardProjectionImageFilter); + /** Maximum number of projections processed together. Zero selects the automatic size (default). */ + itkSetMacro(BatchSize, unsigned int); + itkGetConstMacro(BatchSize, unsigned int); + protected: CudaZengForwardProjectionImageFilter(); ~CudaZengForwardProjectionImageFilter() override; @@ -69,7 +73,8 @@ class RTK_EXPORT CudaZengForwardProjectionImageFilter GPUGenerateData() override; private: - void * m_CudaWorkspace{ nullptr }; + void * m_CudaWorkspace{ nullptr }; + unsigned int m_BatchSize{ 0 }; }; } // namespace rtk diff --git a/include/rtkCudaZengProjectionImageFilter.hcu b/include/rtkCudaZengProjectionImageFilter.hcu index 6148f1dd3..66867b77f 100644 --- a/include/rtkCudaZengProjectionImageFilter.hcu +++ b/include/rtkCudaZengProjectionImageFilter.hcu @@ -33,6 +33,7 @@ CUDA_zeng_forward_project(const int projectionSize[3], const float *devAttenuation, float sigmaZero, float alpha, + unsigned int batchSize, void ** workspace); void RTK_EXPORT @@ -49,6 +50,7 @@ CUDA_zeng_back_project(const int projectionSize[3], const float *devAttenuation, float sigmaZero, float alpha, + unsigned int batchSize, void ** workspace); void RTK_EXPORT diff --git a/src/rtkCudaZengBackProjectionImageFilter.cxx b/src/rtkCudaZengBackProjectionImageFilter.cxx index 4fdc65ceb..868483519 100644 --- a/src/rtkCudaZengBackProjectionImageFilter.cxx +++ b/src/rtkCudaZengBackProjectionImageFilter.cxx @@ -153,6 +153,7 @@ CudaZengBackProjectionImageFilter::GPUGenerateData() attenuation, static_cast(this->GetSigmaZero()), static_cast(this->GetAlpha()), + m_BatchSize, &m_CudaWorkspace); } diff --git a/src/rtkCudaZengForwardProjectionImageFilter.cxx b/src/rtkCudaZengForwardProjectionImageFilter.cxx index 8e7de3cfb..a9f79c2b5 100644 --- a/src/rtkCudaZengForwardProjectionImageFilter.cxx +++ b/src/rtkCudaZengForwardProjectionImageFilter.cxx @@ -151,6 +151,7 @@ CudaZengForwardProjectionImageFilter::GPUGenerateData() attenuation, static_cast(this->GetSigmaZero()), static_cast(this->GetAlpha()), + m_BatchSize, &m_CudaWorkspace); } diff --git a/src/rtkCudaZengProjectionImageFilter.cu b/src/rtkCudaZengProjectionImageFilter.cu index 57d3a714b..92634ed41 100644 --- a/src/rtkCudaZengProjectionImageFilter.cu +++ b/src/rtkCudaZengProjectionImageFilter.cu @@ -27,6 +27,9 @@ namespace { +constexpr int CoefficientStride = 33; +constexpr int MaximumRadius = 32; + inline __device__ float3 applyMatrix(const float * m, float x, float y, float z) { @@ -56,7 +59,9 @@ trilinearZero(const float * image, int3 size, float3 p) for (int dy = 0; dy <= 1; ++dy) for (int dx = 0; dx <= 1; ++dx) { - const int x = x0 + dx, y = y0 + dy, z = z0 + dz; + const int x = x0 + dx; + const int y = y0 + dy; + const int z = z0 + dz; if (x >= 0 && x < size.x && y >= 0 && y < size.y && z >= 0 && z < size.z) { const float wx = dx ? fx : 1.f - fx; @@ -68,104 +73,21 @@ trilinearZero(const float * image, int3 size, float3 p) return result; } -__global__ void -sampleForwardSlice(float * slice, - const float * previous, - const float * volume, - const float * attenuation, - int3 volumeSize, - const float * matrix, - int width, - int height, - int z, - float attenuationStep, - bool addPrevious) -{ - const int x = blockIdx.x * blockDim.x + threadIdx.x; - const int y = blockIdx.y * blockDim.y + threadIdx.y; - if (x >= width || y >= height) - return; - const float3 position = applyMatrix(matrix, x, y, z); - float value = trilinearZero(volume, volumeSize, position); - if (addPrevious) - value += previous[y * width + x]; - if (attenuation) - value *= expf(-attenuationStep * trilinearZero(attenuation, volumeSize, position)); - slice[y * width + x] = value; -} - -__global__ void -gaussianXShared(const float * input, float * output, int width, int height, const float * coefficients, int radius) -{ - extern __shared__ float tile[]; - const int x = blockIdx.x * blockDim.x + threadIdx.x; - const int y = blockIdx.y * blockDim.y + threadIdx.y; - const int tileWidth = blockDim.x + 2 * radius; - const int tileElements = tileWidth * blockDim.y; - const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; - const int threadCount = blockDim.x * blockDim.y; - for (int index = threadIndex; index < tileElements; index += threadCount) - { - const int localX = index % tileWidth; - const int localY = index / tileWidth; - const int globalX = blockIdx.x * blockDim.x + localX - radius; - const int globalY = blockIdx.y * blockDim.y + localY; - tile[index] = globalX >= 0 && globalX < width && globalY < height ? input[globalY * width + globalX] : 0.f; - } - __syncthreads(); - if (x >= width || y >= height) - return; - float sum = 0.f; - for (int offset = -radius; offset <= radius; ++offset) - sum += coefficients[abs(offset)] * tile[threadIdx.y * tileWidth + threadIdx.x + radius + offset]; - output[y * width + x] = sum; -} - -__global__ void -gaussianYShared(const float * input, float * output, int width, int height, const float * coefficients, int radius) -{ - extern __shared__ float tile[]; - const int x = blockIdx.x * blockDim.x + threadIdx.x; - const int y = blockIdx.y * blockDim.y + threadIdx.y; - const int tileHeight = blockDim.y + 2 * radius; - const int tileElements = blockDim.x * tileHeight; - const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; - const int threadCount = blockDim.x * blockDim.y; - for (int index = threadIndex; index < tileElements; index += threadCount) - { - const int localX = index % blockDim.x; - const int localY = index / blockDim.x; - const int globalX = blockIdx.x * blockDim.x + localX; - const int globalY = blockIdx.y * blockDim.y + localY - radius; - tile[index] = globalX < width && globalY >= 0 && globalY < height ? input[globalY * width + globalX] : 0.f; - } - __syncthreads(); - if (x >= width || y >= height) - return; - float sum = 0.f; - for (int offset = -radius; offset <= radius; ++offset) - sum += coefficients[abs(offset)] * tile[(threadIdx.y + radius + offset) * blockDim.x + threadIdx.x]; - output[y * width + x] = sum; -} - std::vector discreteGaussianCoefficients(double variance) { if (variance <= 1.e-12) return { 1.f }; - constexpr double maximumError = 1.e-5; - constexpr size_t maximumKernelWidth = 32; const double exponential = std::exp(-variance); - std::vector coefficients; - coefficients.push_back(exponential * std::cyl_bessel_i(0., variance)); - double sum = coefficients[0]; + std::vector coefficients{ exponential * std::cyl_bessel_i(0., variance) }; + double sum = coefficients[0]; for (int order = 1; sum < 1. - maximumError; ++order) { const double coefficient = exponential * std::cyl_bessel_i(static_cast(order), variance); coefficients.push_back(coefficient); sum += 2. * coefficient; - if (coefficient <= 0. || coefficients.size() > maximumKernelWidth) + if (coefficient <= 0. || coefficients.size() > MaximumRadius) break; } std::vector normalized(coefficients.size()); @@ -175,9 +97,7 @@ discreteGaussianCoefficients(double variance) return normalized; } -constexpr size_t GaussianCoefficientStride = 33; - -struct GaussianMetadata +struct HostGaussianMetadata { std::vector coefficientsX; std::vector coefficientsY; @@ -198,17 +118,68 @@ struct MetadataKey } }; +MetadataKey +makeMetadataKey(int mode, + const int projectionSize[3], + const int volumeSize[3], + const int rotatedSize[3], + const float rotatedSpacing[3], + float sigmaZero, + float alpha, + const float * matrices, + const float * distances, + const int * firstSlices) +{ + MetadataKey key; + key.mode = mode; + key.integers.insert(key.integers.end(), projectionSize, projectionSize + 3); + key.integers.insert(key.integers.end(), volumeSize, volumeSize + 3); + key.integers.insert(key.integers.end(), rotatedSize, rotatedSize + 3); + key.integers.insert(key.integers.end(), firstSlices, firstSlices + projectionSize[2]); + key.values.insert(key.values.end(), rotatedSpacing, rotatedSpacing + 3); + key.values.push_back(sigmaZero); + key.values.push_back(alpha); + key.values.insert(key.values.end(), matrices, matrices + 12 * projectionSize[2]); + key.values.insert(key.values.end(), distances, distances + projectionSize[2]); + return key; +} + +void +appendGaussian(HostGaussianMetadata & metadata, + std::map> & cache, + float variance, + float spacingX, + float spacingY) +{ + const double variances[] = { variance / (spacingX * spacingX), variance / (spacingY * spacingY) }; + std::vector * packed[] = { &metadata.coefficientsX, &metadata.coefficientsY }; + std::vector * radii[] = { &metadata.radiiX, &metadata.radiiY }; + for (int dimension = 0; dimension < 2; ++dimension) + { + auto [it, inserted] = cache.try_emplace(variances[dimension]); + if (inserted) + it->second = discreteGaussianCoefficients(variances[dimension]); + radii[dimension]->push_back(static_cast(it->second.size()) - 1); + packed[dimension]->insert(packed[dimension]->end(), it->second.begin(), it->second.end()); + packed[dimension]->resize(packed[dimension]->size() + CoefficientStride - it->second.size(), 0.f); + } +} + struct DeviceMetadata { - float * coefficientsX{}; - float * coefficientsY{}; - float * matrices{}; - float * inverseMatrices{}; - std::vector radiiX; - std::vector radiiY; + float * coefficientsX{}; + float * coefficientsY{}; + float * matrices{}; + float * inverseMatrices{}; + int * radiiX{}; + int * radiiY{}; + int * firstSlices{}; ~DeviceMetadata() { + cudaFree(firstSlices); + cudaFree(radiiY); + cudaFree(radiiX); cudaFree(inverseMatrices); cudaFree(matrices); cudaFree(coefficientsY); @@ -216,178 +187,413 @@ struct DeviceMetadata } }; -struct CudaZengWorkspace +struct Workspace { float * current{}; float * blurred{}; - float * gaussianScratch{}; + float * scratch{}; float * rotated{}; size_t sliceCapacity{}; size_t rotatedCapacity{}; std::map metadata; - ~CudaZengWorkspace() + ~Workspace() { cudaFree(rotated); - cudaFree(gaussianScratch); + cudaFree(scratch); cudaFree(blurred); cudaFree(current); } }; void -ensureSliceBuffers(CudaZengWorkspace & workspace, size_t requiredElements) +releaseSlices(Workspace & workspace) { - if (requiredElements <= workspace.sliceCapacity) - return; - cudaFree(workspace.gaussianScratch); + cudaFree(workspace.scratch); cudaFree(workspace.blurred); cudaFree(workspace.current); - cudaMalloc(&workspace.current, requiredElements * sizeof(float)); - cudaMalloc(&workspace.blurred, requiredElements * sizeof(float)); - cudaMalloc(&workspace.gaussianScratch, requiredElements * sizeof(float)); - workspace.sliceCapacity = requiredElements; + workspace.scratch = nullptr; + workspace.blurred = nullptr; + workspace.current = nullptr; + workspace.sliceCapacity = 0; } -void -ensureRotatedBuffer(CudaZengWorkspace & workspace, size_t requiredElements) +bool +allocateBuffer(float ** buffer, size_t elements) { - if (requiredElements <= workspace.rotatedCapacity) - return; + const auto error = cudaMalloc(buffer, elements * sizeof(float)); + if (error == cudaErrorMemoryAllocation) + { + cudaGetLastError(); + return false; + } + if (error != cudaSuccess) + itkGenericExceptionMacro(<< "CUDA Zeng allocation failed: " << cudaGetErrorString(error)); + return true; +} + +bool +ensureSlices(Workspace & workspace, size_t elements) +{ + if (elements <= workspace.sliceCapacity) + return true; + releaseSlices(workspace); + if (!allocateBuffer(&workspace.current, elements) || !allocateBuffer(&workspace.blurred, elements) || + !allocateBuffer(&workspace.scratch, elements)) + { + releaseSlices(workspace); + return false; + } + workspace.sliceCapacity = elements; + return true; +} + +bool +ensureRotated(Workspace & workspace, size_t elements) +{ + if (elements <= workspace.rotatedCapacity) + return true; cudaFree(workspace.rotated); - cudaMalloc(&workspace.rotated, requiredElements * sizeof(float)); - workspace.rotatedCapacity = requiredElements; + workspace.rotated = nullptr; + workspace.rotatedCapacity = 0; + if (!allocateBuffer(&workspace.rotated, elements)) + return false; + workspace.rotatedCapacity = elements; + return true; } -MetadataKey -makeMetadataKey(int mode, - const int projectionSize[3], - const int volumeSize[3], - const int rotatedSize[3], - const float rotatedSpacing[3], - float sigmaZero, - float alpha, - const float * matrices, - const float * distances, - const int * firstSlices) +unsigned int +automaticBatchSize(const Workspace & workspace, int projections, int depth, size_t pixels, bool backward) { - MetadataKey key; - key.mode = mode; - key.integers.insert(key.integers.end(), projectionSize, projectionSize + 3); - key.integers.insert(key.integers.end(), volumeSize, volumeSize + 3); - key.integers.insert(key.integers.end(), rotatedSize, rotatedSize + 3); - if (firstSlices) - key.integers.insert(key.integers.end(), firstSlices, firstSlices + projectionSize[2]); - key.values.insert(key.values.end(), rotatedSpacing, rotatedSpacing + 3); - key.values.push_back(sigmaZero); - key.values.push_back(alpha); - key.values.insert(key.values.end(), matrices, matrices + 12 * projectionSize[2]); - key.values.insert(key.values.end(), distances, distances + projectionSize[2]); - return key; + size_t freeBytes = 0; + size_t totalBytes = 0; + cudaMemGetInfo(&freeBytes, &totalBytes); + (void)totalBytes; + const size_t slices = backward ? static_cast(depth) + 3 : 3; + const size_t bytesPerProjection = std::max(1, slices * pixels * sizeof(float)); + size_t existingCapacity = workspace.sliceCapacity / pixels; + if (backward) + existingCapacity = std::min(existingCapacity, workspace.rotatedCapacity / (static_cast(depth) * pixels)); + return static_cast(std::max( + 1, std::min(projections, std::max(freeBytes * 3 / 5 / bytesPerProjection, existingCapacity)))); } -void -appendGaussianMetadata(GaussianMetadata & metadata, - std::map> & coefficientCache, - float variance, - float spacingX, - float spacingY) +unsigned int +allocateBatch(Workspace & workspace, unsigned int batchSize, int depth, size_t pixels, bool backward) { - const double variances[] = { variance / (spacingX * spacingX), variance / (spacingY * spacingY) }; - std::vector * packed[] = { &metadata.coefficientsX, &metadata.coefficientsY }; - std::vector * radii[] = { &metadata.radiiX, &metadata.radiiY }; - for (int dimension = 0; dimension < 2; ++dimension) + for (;;) { - auto [iterator, inserted] = coefficientCache.try_emplace(variances[dimension]); - if (inserted) - iterator->second = discreteGaussianCoefficients(variances[dimension]); - const auto & coefficients = iterator->second; - radii[dimension]->push_back(static_cast(coefficients.size()) - 1); - packed[dimension]->insert(packed[dimension]->end(), coefficients.begin(), coefficients.end()); - packed[dimension]->resize(packed[dimension]->size() + GaussianCoefficientStride - coefficients.size(), 0.f); + if (ensureSlices(workspace, static_cast(batchSize) * pixels) && + (!backward || ensureRotated(workspace, static_cast(batchSize) * depth * pixels))) + return batchSize; + releaseSlices(workspace); + if (batchSize == 1) + itkGenericExceptionMacro(<< "Insufficient GPU memory for one CUDA Zeng projection."); + batchSize = std::max(1u, batchSize / 2); } } void -gaussian2D(const float * input, - float * output, - float * scratch, - const float * deviceCoefficientsX, - const float * deviceCoefficientsY, - int radiusX, - int radiusY, - int width, - int height, - dim3 grid, - dim3 block) +uploadMetadata(DeviceMetadata & device, + HostGaussianMetadata && metadata, + const float * matrices, + const std::vector & inverseMatrices, + const std::vector & firstSlices, + int projections) { - const size_t sharedX = (block.x + 2 * radiusX) * block.y * sizeof(float); - const size_t sharedY = block.x * (block.y + 2 * radiusY) * sizeof(float); - gaussianXShared<<>>(input, scratch, width, height, deviceCoefficientsX, radiusX); - gaussianYShared<<>>(scratch, output, width, height, deviceCoefficientsY, radiusY); + cudaMalloc(&device.coefficientsX, metadata.coefficientsX.size() * sizeof(float)); + cudaMalloc(&device.coefficientsY, metadata.coefficientsY.size() * sizeof(float)); + cudaMalloc(&device.radiiX, metadata.radiiX.size() * sizeof(int)); + cudaMalloc(&device.radiiY, metadata.radiiY.size() * sizeof(int)); + cudaMalloc(&device.matrices, 12 * projections * sizeof(float)); + cudaMalloc(&device.firstSlices, projections * sizeof(int)); + cudaMemcpy(device.coefficientsX, + metadata.coefficientsX.data(), + metadata.coefficientsX.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy(device.coefficientsY, + metadata.coefficientsY.data(), + metadata.coefficientsY.size() * sizeof(float), + cudaMemcpyHostToDevice); + cudaMemcpy(device.radiiX, metadata.radiiX.data(), metadata.radiiX.size() * sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(device.radiiY, metadata.radiiY.data(), metadata.radiiY.size() * sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(device.matrices, matrices, 12 * projections * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(device.firstSlices, firstSlices.data(), projections * sizeof(int), cudaMemcpyHostToDevice); + if (!inverseMatrices.empty()) + { + cudaMalloc(&device.inverseMatrices, inverseMatrices.size() * sizeof(float)); + cudaMemcpy( + device.inverseMatrices, inverseMatrices.data(), inverseMatrices.size() * sizeof(float), cudaMemcpyHostToDevice); + } + CUDA_CHECK_ERROR; } __global__ void -finishForward(const float * input, const float * zeng, float * output, int count, float thickness) +gaussianX(const float * input, + float * output, + int width, + int height, + const float * coefficients, + const int * radii, + const int * firstSlices, + int metadataStride, + int metadataZ, + int batchStart) { - const int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < count) - output[i] = input[i] + thickness * zeng[i]; + extern __shared__ float tile[]; + const int localProjection = blockIdx.z; + if (metadataZ < firstSlices[batchStart + localProjection]) + { + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x < width && y < height) + { + const size_t offset = static_cast(localProjection) * width * height + y * width + x; + output[offset] = input[offset]; + } + return; + } + const int metadataIndex = (batchStart + localProjection) * metadataStride + metadataZ; + const int radius = radii[metadataIndex]; + const int tileWidth = blockDim.x + 2 * radius; + const int tileElements = tileWidth * blockDim.y; + const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; + const int threadCount = blockDim.x * blockDim.y; + const size_t sliceOffset = static_cast(localProjection) * width * height; + for (int index = threadIndex; index < tileElements; index += threadCount) + { + const int localX = index % tileWidth; + const int localY = index / tileWidth; + const int x = blockIdx.x * blockDim.x + localX - radius; + const int y = blockIdx.y * blockDim.y + localY; + tile[index] = x >= 0 && x < width && y < height ? input[sliceOffset + y * width + x] : 0.f; + } + __syncthreads(); + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= width || y >= height) + return; + float sum = 0.f; + const float * kernel = coefficients + static_cast(metadataIndex) * CoefficientStride; + for (int offset = -radius; offset <= radius; ++offset) + sum += kernel[abs(offset)] * tile[threadIdx.y * tileWidth + threadIdx.x + radius + offset]; + output[sliceOffset + y * width + x] = sum; } __global__ void -copyProjection(const float * projections, float * slice, int count, int projection) +gaussianY(const float * input, + float * output, + int width, + int height, + const float * coefficients, + const int * radii, + const int * firstSlices, + int metadataStride, + int metadataZ, + int batchStart) +{ + extern __shared__ float tile[]; + const int localProjection = blockIdx.z; + if (metadataZ < firstSlices[batchStart + localProjection]) + { + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x < width && y < height) + { + const size_t offset = static_cast(localProjection) * width * height + y * width + x; + output[offset] = input[offset]; + } + return; + } + const int metadataIndex = (batchStart + localProjection) * metadataStride + metadataZ; + const int radius = radii[metadataIndex]; + const int tileHeight = blockDim.y + 2 * radius; + const int tileElements = blockDim.x * tileHeight; + const int threadIndex = threadIdx.y * blockDim.x + threadIdx.x; + const int threadCount = blockDim.x * blockDim.y; + const size_t sliceOffset = static_cast(localProjection) * width * height; + for (int index = threadIndex; index < tileElements; index += threadCount) + { + const int localX = index % blockDim.x; + const int localY = index / blockDim.x; + const int x = blockIdx.x * blockDim.x + localX; + const int y = blockIdx.y * blockDim.y + localY - radius; + tile[index] = x < width && y >= 0 && y < height ? input[sliceOffset + y * width + x] : 0.f; + } + __syncthreads(); + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + if (x >= width || y >= height) + return; + float sum = 0.f; + const float * kernel = coefficients + static_cast(metadataIndex) * CoefficientStride; + for (int offset = -radius; offset <= radius; ++offset) + sum += kernel[abs(offset)] * tile[(threadIdx.y + radius + offset) * blockDim.x + threadIdx.x]; + output[sliceOffset + y * width + x] = sum; +} + +void +gaussianBatch(const float * input, + float * output, + float * scratch, + const DeviceMetadata & metadata, + int width, + int height, + int batchStart, + int batchCount, + int metadataZ, + int metadataStride) { - const int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < count) - slice[i] = projections[projection * count + i]; + const dim3 block(16, 16); + const dim3 grid(iDivUp(width, 16), iDivUp(height, 16), batchCount); + const size_t sharedX = (block.x + 2 * MaximumRadius) * block.y * sizeof(float); + const size_t sharedY = block.x * (block.y + 2 * MaximumRadius) * sizeof(float); + gaussianX<<>>(input, + scratch, + width, + height, + metadata.coefficientsX, + metadata.radiiX, + metadata.firstSlices, + metadataStride, + metadataZ, + batchStart); + gaussianY<<>>(scratch, + output, + width, + height, + metadata.coefficientsY, + metadata.radiiY, + metadata.firstSlices, + metadataStride, + metadataZ, + batchStart); } __global__ void -attenuateSlice(float * slice, - cudaTextureObject_t attenuation, - int3 volumeSize, - const float * rotatedToVolume, - int width, - int height, - int z, - float attenuationStep) +sampleForward(float * current, + const float * previous, + const float * volume, + const float * attenuation, + int3 volumeSize, + const float * matrices, + const int * firstSlices, + int width, + int height, + int z, + float attenuationStep, + int batchStart, + bool addPrevious) { const int x = blockIdx.x * blockDim.x + threadIdx.x; const int y = blockIdx.y * blockDim.y + threadIdx.y; - if (x >= width || y >= height) + const int localProjection = blockIdx.z; + const int projection = batchStart + localProjection; + if (x >= width || y >= height || z < firstSlices[projection]) return; - const float3 position = applyMatrix(rotatedToVolume, x, y, z); - if (position.x >= 0.f && position.x < volumeSize.x && position.y >= 0.f && position.y < volumeSize.y && - position.z >= 0.f && position.z < volumeSize.z) - slice[y * width + x] *= expf(-attenuationStep * tex3D(attenuation, position.x, position.y, position.z)); + const int pixel = y * width + x; + const size_t offset = static_cast(localProjection) * width * height + pixel; + const float3 position = applyMatrix(matrices + 12 * projection, x, y, z); + float value = trilinearZero(volume, volumeSize, position); + if (addPrevious) + value += previous[offset]; + if (attenuation) + value *= expf(-attenuationStep * trilinearZero(attenuation, volumeSize, position)); + current[offset] = value; +} + +__global__ void +finishForward(const float * input, const float * zeng, float * output, int pixels, float thickness, int batchStart) +{ + const int pixel = blockIdx.x * blockDim.x + threadIdx.x; + const int localProjection = blockIdx.y; + if (pixel >= pixels) + return; + const int projection = batchStart + localProjection; + const size_t local = static_cast(localProjection) * pixels + pixel; + output[static_cast(projection) * pixels + pixel] = + input[static_cast(projection) * pixels + pixel] + thickness * zeng[local]; +} + +__global__ void +copyProjections(const float * projections, float * current, int pixels, int batchStart) +{ + const int pixel = blockIdx.x * blockDim.x + threadIdx.x; + const int localProjection = blockIdx.y; + if (pixel < pixels) + current[static_cast(localProjection) * pixels + pixel] = + projections[static_cast(batchStart + localProjection) * pixels + pixel]; +} + +__global__ void +storeSlices(const float * current, + float * rotated, + int pixels, + int depth, + int z, + const int * firstSlices, + int batchStart) +{ + const int pixel = blockIdx.x * blockDim.x + threadIdx.x; + const int localProjection = blockIdx.y; + const int projection = batchStart + localProjection; + if (pixel < pixels && z >= firstSlices[projection]) + rotated[(static_cast(localProjection) * depth + z) * pixels + pixel] = + current[static_cast(localProjection) * pixels + pixel]; } __global__ void -storeSlice(const float * slice, float * volume, int count, int z) +attenuate(float * current, + cudaTextureObject_t attenuation, + int3 volumeSize, + const float * inverseMatrices, + int width, + int height, + int z, + const int * firstSlices, + bool firstSlice, + float step, + int batchStart) { - const int i = blockIdx.x * blockDim.x + threadIdx.x; - if (i < count) - volume[z * count + i] = slice[i]; + const int x = blockIdx.x * blockDim.x + threadIdx.x; + const int y = blockIdx.y * blockDim.y + threadIdx.y; + const int localProjection = blockIdx.z; + const int projection = batchStart + localProjection; + if (x >= width || y >= height || (!firstSlice && z <= firstSlices[projection])) + return; + const int slice = firstSlice ? firstSlices[projection] : z; + const float3 position = applyMatrix(inverseMatrices + 12 * projection, x, y, slice); + const size_t offset = static_cast(localProjection) * width * height + y * width + x; + if (position.x >= 0.f && position.x < volumeSize.x && position.y >= 0.f && position.y < volumeSize.y && + position.z >= 0.f && position.z < volumeSize.z) + current[offset] *= expf(-step * tex3D(attenuation, position.x, position.y, position.z)); } __global__ void -addRotatedVolume(const float * rotated, - const float * input, - float * output, - int3 volumeSize, - int3 rotatedSize, - const float * volumeToRotated, - float thickness) +addRotatedBatch(const float * rotated, + float * output, + int3 volumeSize, + int3 rotatedSize, + const float * matrices, + int batchStart, + int batchCount, + float thickness) { const int x = blockIdx.x * blockDim.x + threadIdx.x; const int y = blockIdx.y * blockDim.y + threadIdx.y; const int z = blockIdx.z * blockDim.z + threadIdx.z; if (x >= volumeSize.x || y >= volumeSize.y || z >= volumeSize.z) return; - const int index = (z * volumeSize.y + y) * volumeSize.x + x; - const float3 position = applyMatrix(volumeToRotated, x, y, z); - output[index] = input[index] + thickness * trilinearZero(rotated, rotatedSize, position); + float sum = 0.f; + const size_t rotatedElements = static_cast(rotatedSize.x) * rotatedSize.y * rotatedSize.z; + for (int localProjection = 0; localProjection < batchCount; ++localProjection) + { + const int projection = batchStart + localProjection; + const float3 position = applyMatrix(matrices + 12 * projection, x, y, z); + sum += trilinearZero(rotated + localProjection * rotatedElements, rotatedSize, position); + } + const int index = (z * volumeSize.y + y) * volumeSize.x + x; + output[index] += thickness * sum; } void @@ -423,15 +629,30 @@ CUDA_zeng_forward_project(const int projectionSize[3], const float * devAttenuation, float sigmaZero, float alpha, + unsigned int requestedBatchSize, void ** workspacePointer) { - const int3 cudaVolumeSize = make_int3(volumeSize[0], volumeSize[1], volumeSize[2]); if (!*workspacePointer) - *workspacePointer = new CudaZengWorkspace; - auto & workspace = *static_cast(*workspacePointer); - const int pixelCount = rotatedSize[0] * rotatedSize[1]; - ensureSliceBuffers(workspace, pixelCount); - + *workspacePointer = new Workspace; + auto & workspace = *static_cast(*workspacePointer); + const int projections = projectionSize[2]; + const int depth = rotatedSize[2]; + const int pixels = rotatedSize[0] * rotatedSize[1]; + const int metadataStride = depth + 1; + std::vector firstSlices(projections); + std::vector nearDistances(projections); + for (int projection = 0; projection < projections; ++projection) + { + float nearDistance = farDistances[projection]; + int first = depth - 1; + while (first > 0 && nearDistance - rotatedSpacing[2] >= 0.f) + { + --first; + nearDistance -= rotatedSpacing[2]; + } + firstSlices[projection] = first; + nearDistances[projection] = nearDistance; + } const auto key = makeMetadataKey(0, projectionSize, volumeSize, @@ -441,112 +662,94 @@ CUDA_zeng_forward_project(const int projectionSize[3], alpha, rotatedToVolumeMatrices, farDistances, - nullptr); - auto [metadataIterator, inserted] = workspace.metadata.try_emplace(key); - auto & deviceMetadata = metadataIterator->second; + firstSlices.data()); + auto [iterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = iterator->second; if (inserted) { - GaussianMetadata metadata; - std::map> coefficientCache; - for (int projection = 0; projection < projectionSize[2]; ++projection) + HostGaussianMetadata metadata; + std::map> cache; + for (int projection = 0; projection < projections; ++projection) { - float distance = farDistances[projection]; - for (int z = rotatedSize[2] - 2; z >= 0 && distance - rotatedSpacing[2] >= 0.f; --z) + const float nearDistance = nearDistances[projection]; + for (int z = 0; z < depth; ++z) { - const float variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + - 2.f * rotatedSpacing[2] * alpha * sigmaZero - - alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; - appendGaussianMetadata( - metadata, coefficientCache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); - distance -= rotatedSpacing[2]; + float variance = 0.f; + if (z >= firstSlices[projection] && z + 1 < depth) + { + const float distance = nearDistance + (z + 1 - firstSlices[projection]) * rotatedSpacing[2]; + variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + 2.f * rotatedSpacing[2] * alpha * sigmaZero - + alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; + } + appendGaussian(metadata, cache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); } - const float finalVariance = (alpha * distance + sigmaZero) * (alpha * distance + sigmaZero); - appendGaussianMetadata(metadata, coefficientCache, finalVariance, rotatedSpacing[0], rotatedSpacing[1]); + const float finalVariance = (alpha * nearDistance + sigmaZero) * (alpha * nearDistance + sigmaZero); + appendGaussian(metadata, cache, finalVariance, rotatedSpacing[0], rotatedSpacing[1]); } - deviceMetadata.radiiX = std::move(metadata.radiiX); - deviceMetadata.radiiY = std::move(metadata.radiiY); - cudaMalloc(&deviceMetadata.coefficientsX, metadata.coefficientsX.size() * sizeof(float)); - cudaMalloc(&deviceMetadata.coefficientsY, metadata.coefficientsY.size() * sizeof(float)); - cudaMalloc(&deviceMetadata.matrices, 12 * projectionSize[2] * sizeof(float)); - cudaMemcpy(deviceMetadata.coefficientsX, - metadata.coefficientsX.data(), - metadata.coefficientsX.size() * sizeof(float), - cudaMemcpyHostToDevice); - cudaMemcpy(deviceMetadata.coefficientsY, - metadata.coefficientsY.data(), - metadata.coefficientsY.size() * sizeof(float), - cudaMemcpyHostToDevice); - cudaMemcpy( - deviceMetadata.matrices, rotatedToVolumeMatrices, 12 * projectionSize[2] * sizeof(float), cudaMemcpyHostToDevice); + uploadMetadata(deviceMetadata, std::move(metadata), rotatedToVolumeMatrices, {}, firstSlices, projections); } - float * current = workspace.current; - float * blurred = workspace.blurred; - float * gaussianScratch = workspace.gaussianScratch; - const dim3 block(16, 16); - const dim3 grid(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16)); - - size_t gaussianIndex = 0; - for (int projection = 0; projection < projectionSize[2]; ++projection) + const unsigned int desiredBatchSize = requestedBatchSize + ? std::min(requestedBatchSize, static_cast(projections)) + : automaticBatchSize(workspace, projections, depth, pixels, false); + const unsigned int batchSize = allocateBatch(workspace, desiredBatchSize, depth, pixels, false); + const dim3 block(16, 16); + const int3 cudaVolumeSize = make_int3(volumeSize[0], volumeSize[1], volumeSize[2]); + for (int batchStart = 0; batchStart < projections; batchStart += batchSize) { - const float * matrix = deviceMetadata.matrices + 12 * projection; - float distance = farDistances[projection]; - int z = rotatedSize[2] - 1; - sampleForwardSlice<<>>(current, - nullptr, - devVolume, - devAttenuation, - cudaVolumeSize, - matrix, - rotatedSize[0], - rotatedSize[1], - z, - rotatedSpacing[2], - false); - for (--z; z >= 0 && distance - rotatedSpacing[2] >= 0.f; --z) + const int batchCount = std::min(batchSize, projections - batchStart); + const dim3 grid(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16), batchCount); + sampleForward<<>>(workspace.current, + nullptr, + devVolume, + devAttenuation, + cudaVolumeSize, + deviceMetadata.matrices, + deviceMetadata.firstSlices, + rotatedSize[0], + rotatedSize[1], + depth - 1, + rotatedSpacing[2], + batchStart, + false); + for (int z = depth - 2; z >= 0; --z) { - gaussian2D(current, - blurred, - gaussianScratch, - deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.radiiX[gaussianIndex], - deviceMetadata.radiiY[gaussianIndex], - rotatedSize[0], - rotatedSize[1], - grid, - block); - ++gaussianIndex; - sampleForwardSlice<<>>(current, - blurred, - devVolume, - devAttenuation, - cudaVolumeSize, - matrix, - rotatedSize[0], - rotatedSize[1], - z, - rotatedSpacing[2], - true); - distance -= rotatedSpacing[2]; + gaussianBatch(workspace.current, + workspace.blurred, + workspace.scratch, + deviceMetadata, + rotatedSize[0], + rotatedSize[1], + batchStart, + batchCount, + z, + metadataStride); + sampleForward<<>>(workspace.current, + workspace.blurred, + devVolume, + devAttenuation, + cudaVolumeSize, + deviceMetadata.matrices, + deviceMetadata.firstSlices, + rotatedSize[0], + rotatedSize[1], + z, + rotatedSpacing[2], + batchStart, + true); } - gaussian2D(current, - blurred, - gaussianScratch, - deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.radiiX[gaussianIndex], - deviceMetadata.radiiY[gaussianIndex], - rotatedSize[0], - rotatedSize[1], - grid, - block); - ++gaussianIndex; - finishForward<<>>(devProjectionIn + projection * pixelCount, - blurred, - devProjectionOut + projection * pixelCount, - pixelCount, - rotatedSpacing[2]); + gaussianBatch(workspace.current, + workspace.blurred, + workspace.scratch, + deviceMetadata, + rotatedSize[0], + rotatedSize[1], + batchStart, + batchCount, + depth, + metadataStride); + finishForward<<>>( + devProjectionIn, workspace.blurred, devProjectionOut, pixels, rotatedSpacing[2], batchStart); } CUDA_CHECK_ERROR; } @@ -565,19 +768,18 @@ CUDA_zeng_back_project(const int projectionSize[3], const float * devAttenuation, float sigmaZero, float alpha, + unsigned int requestedBatchSize, void ** workspacePointer) { - const size_t volumeBytes = static_cast(volumeSize[0]) * volumeSize[1] * volumeSize[2] * sizeof(float); if (!*workspacePointer) - *workspacePointer = new CudaZengWorkspace; - auto & workspace = *static_cast(*workspacePointer); - const int slicePixels = rotatedSize[0] * rotatedSize[1]; - const size_t rotatedElements = static_cast(slicePixels) * rotatedSize[2]; - const size_t rotatedBytes = rotatedElements * sizeof(float); - ensureSliceBuffers(workspace, slicePixels); - ensureRotatedBuffer(workspace, rotatedElements); - - const auto key = makeMetadataKey(1, + *workspacePointer = new Workspace; + auto & workspace = *static_cast(*workspacePointer); + const int projections = projectionSize[2]; + const int depth = rotatedSize[2]; + const int pixels = rotatedSize[0] * rotatedSize[1]; + const int metadataStride = depth + 1; + std::vector first(firstSlices, firstSlices + projections); + const auto key = makeMetadataKey(1, projectionSize, volumeSize, rotatedSize, @@ -587,57 +789,41 @@ CUDA_zeng_back_project(const int projectionSize[3], volumeToRotatedMatrices, nearDistances, firstSlices); - auto [metadataIterator, inserted] = workspace.metadata.try_emplace(key); - auto & deviceMetadata = metadataIterator->second; + auto [iterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = iterator->second; if (inserted) { - GaussianMetadata metadata; - std::map> coefficientCache; - std::vector inverseMatrices(12 * projectionSize[2]); - for (int projection = 0; projection < projectionSize[2]; ++projection) + HostGaussianMetadata metadata; + std::map> cache; + std::vector inverseMatrices(12 * projections); + for (int projection = 0; projection < projections; ++projection) { - float distance = nearDistances[projection]; - float variance = (alpha * distance + sigmaZero) * (alpha * distance + sigmaZero); - appendGaussianMetadata(metadata, coefficientCache, variance, rotatedSpacing[0], rotatedSpacing[1]); - for (int z = firstSlices[projection]; z + 1 < rotatedSize[2]; ++z) + for (int z = 0; z < depth; ++z) { - distance += rotatedSpacing[2]; - variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + 2.f * rotatedSpacing[2] * alpha * sigmaZero - - alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; - appendGaussianMetadata( - metadata, coefficientCache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); + float variance = 0.f; + if (z >= firstSlices[projection] && z + 1 < depth) + { + const float distance = nearDistances[projection] + (z + 1 - firstSlices[projection]) * rotatedSpacing[2]; + variance = distance * 2.f * rotatedSpacing[2] * alpha * alpha + 2.f * rotatedSpacing[2] * alpha * sigmaZero - + alpha * alpha * rotatedSpacing[2] * rotatedSpacing[2]; + } + appendGaussian(metadata, cache, std::max(0.f, variance), rotatedSpacing[0], rotatedSpacing[1]); } + const float initialVariance = + (alpha * nearDistances[projection] + sigmaZero) * (alpha * nearDistances[projection] + sigmaZero); + appendGaussian(metadata, cache, initialVariance, rotatedSpacing[0], rotatedSpacing[1]); invertAffine(volumeToRotatedMatrices + 12 * projection, inverseMatrices.data() + 12 * projection); } - deviceMetadata.radiiX = std::move(metadata.radiiX); - deviceMetadata.radiiY = std::move(metadata.radiiY); - cudaMalloc(&deviceMetadata.coefficientsX, metadata.coefficientsX.size() * sizeof(float)); - cudaMalloc(&deviceMetadata.coefficientsY, metadata.coefficientsY.size() * sizeof(float)); - cudaMalloc(&deviceMetadata.matrices, 12 * projectionSize[2] * sizeof(float)); - cudaMalloc(&deviceMetadata.inverseMatrices, 12 * projectionSize[2] * sizeof(float)); - cudaMemcpy(deviceMetadata.coefficientsX, - metadata.coefficientsX.data(), - metadata.coefficientsX.size() * sizeof(float), - cudaMemcpyHostToDevice); - cudaMemcpy(deviceMetadata.coefficientsY, - metadata.coefficientsY.data(), - metadata.coefficientsY.size() * sizeof(float), - cudaMemcpyHostToDevice); - cudaMemcpy( - deviceMetadata.matrices, volumeToRotatedMatrices, 12 * projectionSize[2] * sizeof(float), cudaMemcpyHostToDevice); - cudaMemcpy(deviceMetadata.inverseMatrices, - inverseMatrices.data(), - inverseMatrices.size() * sizeof(float), - cudaMemcpyHostToDevice); + uploadMetadata(deviceMetadata, std::move(metadata), volumeToRotatedMatrices, inverseMatrices, first, projections); } - float * rotated = workspace.rotated; - float * current = workspace.current; - float * blurred = workspace.blurred; - float * gaussianScratch = workspace.gaussianScratch; + const unsigned int desiredBatchSize = requestedBatchSize + ? std::min(requestedBatchSize, static_cast(projections)) + : automaticBatchSize(workspace, projections, depth, pixels, true); + const unsigned int batchSize = allocateBatch(workspace, desiredBatchSize, depth, pixels, true); + const size_t volumeBytes = static_cast(volumeSize[0]) * volumeSize[1] * volumeSize[2] * sizeof(float); if (devVolumeOut != devVolumeIn) cudaMemcpy(devVolumeOut, devVolumeIn, volumeBytes, cudaMemcpyDeviceToDevice); - cudaArray * attenuationArray = nullptr; cudaTextureObject_t attenuationTexture = 0; if (devAttenuation) @@ -648,78 +834,81 @@ CUDA_zeng_back_project(const int projectionSize[3], false, true, cudaAddressModeClamp); - const dim3 block2(16, 16); - const dim3 grid2(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16)); const dim3 block3(8, 8, 4); const dim3 grid3(iDivUp(volumeSize[0], 8), iDivUp(volumeSize[1], 8), iDivUp(volumeSize[2], 4)); const int3 cudaVolumeSize = make_int3(volumeSize[0], volumeSize[1], volumeSize[2]); const int3 cudaRotatedSize = make_int3(rotatedSize[0], rotatedSize[1], rotatedSize[2]); - - size_t gaussianIndex = 0; - for (int projection = 0; projection < projectionSize[2]; ++projection) + for (int batchStart = 0; batchStart < projections; batchStart += batchSize) { - const float * matrix = deviceMetadata.matrices + 12 * projection; - const float * inverseMatrix = deviceMetadata.inverseMatrices + 12 * projection; - cudaMemset(rotated, 0, rotatedBytes); - copyProjection<<>>(devProjections, current, slicePixels, projection); + const int batchCount = std::min(batchSize, projections - batchStart); + float * current = workspace.current; + float * blurred = workspace.blurred; + const dim3 grid2(iDivUp(rotatedSize[0], 16), iDivUp(rotatedSize[1], 16), batchCount); + cudaMemset(workspace.rotated, 0, static_cast(batchCount) * depth * pixels * sizeof(float)); + copyProjections<<>>(devProjections, current, pixels, batchStart); if (attenuationTexture) - attenuateSlice<<>>(current, - attenuationTexture, - cudaVolumeSize, - inverseMatrix, - rotatedSize[0], - rotatedSize[1], - firstSlices[projection], - rotatedSpacing[2]); - float distance = nearDistances[projection]; - gaussian2D(current, - blurred, - gaussianScratch, - deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.radiiX[gaussianIndex], - deviceMetadata.radiiY[gaussianIndex], - rotatedSize[0], - rotatedSize[1], - grid2, - block2); - ++gaussianIndex; + attenuate<<>>(current, + attenuationTexture, + cudaVolumeSize, + deviceMetadata.inverseMatrices, + rotatedSize[0], + rotatedSize[1], + 0, + deviceMetadata.firstSlices, + true, + rotatedSpacing[2], + batchStart); + gaussianBatch(current, + blurred, + workspace.scratch, + deviceMetadata, + rotatedSize[0], + rotatedSize[1], + batchStart, + batchCount, + depth, + metadataStride); std::swap(current, blurred); - - for (int z = firstSlices[projection]; z < rotatedSize[2]; ++z) + for (int z = 0; z < depth; ++z) { - storeSlice<<>>(current, rotated, slicePixels, z); - if (z + 1 == rotatedSize[2]) + storeSlices<<>>( + current, workspace.rotated, pixels, depth, z, deviceMetadata.firstSlices, batchStart); + if (z + 1 == depth) break; if (attenuationTexture) - attenuateSlice<<>>(current, - attenuationTexture, - cudaVolumeSize, - inverseMatrix, - rotatedSize[0], - rotatedSize[1], - z + 1, - rotatedSpacing[2]); - distance += rotatedSpacing[2]; - gaussian2D(current, - blurred, - gaussianScratch, - deviceMetadata.coefficientsX + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.coefficientsY + GaussianCoefficientStride * gaussianIndex, - deviceMetadata.radiiX[gaussianIndex], - deviceMetadata.radiiY[gaussianIndex], - rotatedSize[0], - rotatedSize[1], - grid2, - block2); - ++gaussianIndex; + attenuate<<>>(current, + attenuationTexture, + cudaVolumeSize, + deviceMetadata.inverseMatrices, + rotatedSize[0], + rotatedSize[1], + z + 1, + deviceMetadata.firstSlices, + false, + rotatedSpacing[2], + batchStart); + gaussianBatch(current, + blurred, + workspace.scratch, + deviceMetadata, + rotatedSize[0], + rotatedSize[1], + batchStart, + batchCount, + z, + metadataStride); std::swap(current, blurred); } - addRotatedVolume<<>>( - rotated, devVolumeOut, devVolumeOut, cudaVolumeSize, cudaRotatedSize, matrix, rotatedSpacing[2]); + addRotatedBatch<<>>(workspace.rotated, + devVolumeOut, + cudaVolumeSize, + cudaRotatedSize, + deviceMetadata.matrices, + batchStart, + batchCount, + rotatedSpacing[2]); } - CUDA_CHECK_ERROR; if (attenuationArray) { cudaDestroyTextureObject(attenuationTexture); @@ -731,5 +920,5 @@ CUDA_zeng_back_project(const int projectionSize[3], void CUDA_zeng_release_workspace(void * workspace) { - delete static_cast(workspace); + delete static_cast(workspace); }