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..4e907dbad 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", @@ -130,11 +131,19 @@ def process(args_info: argparse.Namespace): bp = rtk.JosephBackAttenuatedProjectionImageFilter[ OutputImageType, OutputImageType ].New() - elif args_info.bp == "Zeng": - bp = rtk.ZengBackProjectionImageFilter[OutputImageType, OutputImageType].New() - if args_info.sigmazero: + elif args_info.bp in ("Zeng", "CudaZeng"): + if args_info.bp == "CudaZeng": + if not hasattr(itk, "CudaImage"): + print("The program has not been compiled with cuda option") + sys.exit(1) + bp = rtk.CudaZengBackProjectionImageFilter.New() + else: + bp = rtk.ZengBackProjectionImageFilter[ + OutputImageType, OutputImageType + ].New() + if args_info.sigmazero is not None: bp.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf: + if args_info.alphapsf is not None: bp.SetAlpha(args_info.alphapsf) elif args_info.bp == "CudaFDKBackProjection": @@ -156,7 +165,12 @@ def process(args_info: argparse.Namespace): print("The program has not been compiled with cuda option") sys.exit(1) - if args_info.bp in ["CudaFDKBackProjection", "CudaBackProjection", "CudaRayCast"]: + if args_info.bp in [ + "CudaFDKBackProjection", + "CudaBackProjection", + "CudaRayCast", + "CudaZeng", + ]: bp.SetInput(itk.cuda_image_from_image(constantImageSource.GetOutput())) bp.SetInput(1, itk.cuda_image_from_image(reader.GetOutput())) if attenuation_map: 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..e1f3bd043 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.py +++ b/applications/rtkforwardprojections/rtkforwardprojections.py @@ -44,6 +44,7 @@ def build_parser(): "JosephAttenuated", "CudaRayCast", "Zeng", + "CudaZeng", "MIP", "CudaWrapRayCast", ], @@ -157,10 +158,16 @@ def process(args_info): "Please provide it using --attenuationmap." ) sys.exit(1) - elif args_info.fp == "Zeng": - forwardProjection = rtk.ZengForwardProjectionImageFilter[ - OutputImageType, OutputImageType - ].New() + elif args_info.fp in ("Zeng", "CudaZeng"): + if args_info.fp == "CudaZeng": + if not hasattr(itk, "CudaImage"): + print("The program has not been compiled with cuda option") + sys.exit(1) + forwardProjection = rtk.CudaZengForwardProjectionImageFilter.New() + else: + forwardProjection = rtk.ZengForwardProjectionImageFilter[ + OutputImageType, OutputImageType + ].New() elif args_info.fp == "MIP": forwardProjection = rtk.MaximumIntensityProjectionImageFilter[ OutputImageType, OutputImageType @@ -191,7 +198,7 @@ def process(args_info): print("The program has not been compiled with cuda option") sys.exit(1) - if args_info.fp in ["CudaWrapRayCast", "CudaRayCast"]: + if args_info.fp in ["CudaWrapRayCast", "CudaRayCast", "CudaZeng"]: forwardProjection.SetInput( itk.cuda_image_from_image(constantImageSource.GetOutput()) ) @@ -218,9 +225,9 @@ def process(args_info): forwardProjection.SetSuperiorClipImage(superiorClipImage) elif args_info.fp == "MIP": forwardProjection.SetSuperiorClipImage(superiorClipImage) - if args_info.sigmazero and args_info.fp == "Zeng": + if args_info.sigmazero and args_info.fp in ["Zeng", "CudaZeng"]: forwardProjection.SetSigmaZero(args_info.sigmazero) - if args_info.alphapsf and args_info.fp == "Zeng": + if args_info.alphapsf and args_info.fp in ["Zeng", "CudaZeng"]: forwardProjection.SetAlpha(args_info.alphapsf) forwardProjection.SetGeometry(geometry) 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..72657f85b 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", ) @@ -37,16 +38,16 @@ def add_rtkprojectors_group(parser): ) rtkprojectors_group.add_argument( "--attenuationmap", - help="Attenuation map relative to the volume to perfom the attenuation correction (JosephAttenuated and Zeng)", + help="Attenuation map relative to the volume to perform the attenuation correction (JosephAttenuated, Zeng and CudaZeng)", ) rtkprojectors_group.add_argument( "--sigmazero", - help="PSF value at a distance of 0 meter of the detector (Zeng only)", + help="PSF value at a distance of 0 meter of the detector (Zeng and CudaZeng only)", type=float, ) rtkprojectors_group.add_argument( "--alphapsf", - help="Slope of the PSF against the detector distance (Zeng only)", + help="Slope of the PSF against the detector distance (Zeng and CudaZeng only)", type=float, ) rtkprojectors_group.add_argument( @@ -89,8 +90,13 @@ def SetBackProjectionFromArgParse(args_info, recon): if args_info.superiorclipimage is not None: recon.SetSuperiorClipImage(superior_clip_image) recon.SetAttenuationMap(attenuation_map) - elif args_info.bp == "Zeng": # bp_arg_RotationBased - recon.SetBackProjectionFilter(ReconType.BackProjectionType_BP_ZENG) + elif args_info.bp in ("Zeng", "CudaZeng"): + projector = ( + ReconType.BackProjectionType_BP_ZENG + if args_info.bp == "Zeng" + else ReconType.BackProjectionType_BP_CUDAZENG + ) + recon.SetBackProjectionFilter(projector) if args_info.sigmazero is not None: recon.SetSigmaZero(args_info.sigmazero) if args_info.alphapsf is not None: @@ -119,8 +125,13 @@ def SetForwardProjectionFromArgParse(args_info, recon): ReconType.ForwardProjectionType_FP_JOSEPHATTENUATED ) recon.SetAttenuationMap(attenuation_map) - elif args_info.fp == "Zeng": # fp_arg_RotationBased - recon.SetForwardProjectionFilter(ReconType.ForwardProjectionType_FP_ZENG) + elif args_info.fp in ("Zeng", "CudaZeng"): + projector = ( + ReconType.ForwardProjectionType_FP_ZENG + if args_info.fp == "Zeng" + else ReconType.ForwardProjectionType_FP_CUDAZENG + ) + recon.SetForwardProjectionFilter(projector) if args_info.sigmazero is not None: recon.SetSigmaZero(args_info.sigmazero) if args_info.alphapsf is not None: 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..2f3ee4ddf --- /dev/null +++ b/include/rtkCudaZengBackProjectionImageFilter.h @@ -0,0 +1,82 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef rtkCudaZengBackProjectionImageFilter_h +#define rtkCudaZengBackProjectionImageFilter_h + +#include "rtkConfiguration.h" +#ifdef RTK_USE_CUDA + +# include "rtkZengBackProjectionImageFilter.h" +# include "RTKExport.h" +# include +# 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). + * 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, + itk::CudaImage, + ZengBackProjectionImageFilter, 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 = ZengBackProjectionImageFilter; + using Self = CudaZengBackProjectionImageFilter; + using Superclass = itk::CudaInPlaceImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + itkNewMacro(Self); + itkOverrideGetNameOfClassMacro(CudaZengBackProjectionImageFilter); + + /** Maximum number of projections processed together. Zero selects the automatic size (default). */ + itkSetMacro(BatchSize, unsigned int); + itkGetConstMacro(BatchSize, unsigned int); + +protected: + CudaZengBackProjectionImageFilter(); + ~CudaZengBackProjectionImageFilter() override; + void + GPUGenerateData() override; + +private: + void * m_CudaWorkspace{ nullptr }; + unsigned int m_BatchSize{ 0 }; +}; + +} // namespace rtk +#endif // RTK_USE_CUDA +#endif // rtkCudaZengBackProjectionImageFilter_h diff --git a/include/rtkCudaZengForwardProjectionImageFilter.h b/include/rtkCudaZengForwardProjectionImageFilter.h new file mode 100644 index 000000000..d32c3a532 --- /dev/null +++ b/include/rtkCudaZengForwardProjectionImageFilter.h @@ -0,0 +1,82 @@ +/*========================================================================= + * + * Copyright RTK Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef rtkCudaZengForwardProjectionImageFilter_h +#define rtkCudaZengForwardProjectionImageFilter_h + +#include "rtkConfiguration.h" +#ifdef RTK_USE_CUDA + +# include "rtkZengForwardProjectionImageFilter.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). + * 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, + itk::CudaImage, + ZengForwardProjectionImageFilter, 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 = ZengForwardProjectionImageFilter; + using Self = CudaZengForwardProjectionImageFilter; + using Superclass = itk::CudaInPlaceImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + 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; + void + GPUGenerateData() override; + +private: + void * m_CudaWorkspace{ nullptr }; + unsigned int m_BatchSize{ 0 }; +}; + +} // 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..66867b77f --- /dev/null +++ b/include/rtkCudaZengProjectionImageFilter.hcu @@ -0,0 +1,59 @@ +/*========================================================================= + * + * 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, + unsigned int batchSize, + 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, + unsigned int batchSize, + 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..868483519 --- /dev/null +++ b/src/rtkCudaZengBackProjectionImageFilter.cxx @@ -0,0 +1,161 @@ +/*========================================================================= + * + * 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; +} // namespace + +CudaZengBackProjectionImageFilter::CudaZengBackProjectionImageFilter() = default; + +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(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."); + 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]))); + 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(this->GetSigmaZero()), + static_cast(this->GetAlpha()), + m_BatchSize, + &m_CudaWorkspace); +} + +} // namespace rtk +#endif diff --git a/src/rtkCudaZengForwardProjectionImageFilter.cxx b/src/rtkCudaZengForwardProjectionImageFilter.cxx new file mode 100644 index 000000000..a9f79c2b5 --- /dev/null +++ b/src/rtkCudaZengForwardProjectionImageFilter.cxx @@ -0,0 +1,159 @@ +/*========================================================================= + * + * 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; +} // namespace + +CudaZengForwardProjectionImageFilter::CudaZengForwardProjectionImageFilter() = default; + +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(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); + + 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.GetVnlMatrix() * + GetIndexToPhysicalPointMatrix(rotatedImage.GetPointer()).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]); + matrices[12 * local + 4 * r + 3] -= volumeRegion.GetIndex(r); + } + 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(this->GetSigmaZero()), + static_cast(this->GetAlpha()), + m_BatchSize, + &m_CudaWorkspace); +} + +} // namespace rtk +#endif diff --git a/src/rtkCudaZengProjectionImageFilter.cu b/src/rtkCudaZengProjectionImageFilter.cu new file mode 100644 index 000000000..92634ed41 --- /dev/null +++ b/src/rtkCudaZengProjectionImageFilter.cu @@ -0,0 +1,924 @@ +/*========================================================================= + * + * 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 +{ +constexpr int CoefficientStride = 33; +constexpr int MaximumRadius = 32; + +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) +{ + // 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)); + 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; + 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; + 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; +} + +std::vector +discreteGaussianCoefficients(double variance) +{ + if (variance <= 1.e-12) + return { 1.f }; + constexpr double maximumError = 1.e-5; + const double exponential = std::exp(-variance); + 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() > MaximumRadius) + break; + } + std::vector normalized(coefficients.size()); + std::transform(coefficients.begin(), coefficients.end(), normalized.begin(), [sum](double value) { + return static_cast(value / sum); + }); + return normalized; +} + +struct HostGaussianMetadata +{ + 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); + } +}; + +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{}; + int * radiiX{}; + int * radiiY{}; + int * firstSlices{}; + + ~DeviceMetadata() + { + cudaFree(firstSlices); + cudaFree(radiiY); + cudaFree(radiiX); + cudaFree(inverseMatrices); + cudaFree(matrices); + cudaFree(coefficientsY); + cudaFree(coefficientsX); + } +}; + +struct Workspace +{ + float * current{}; + float * blurred{}; + float * scratch{}; + float * rotated{}; + size_t sliceCapacity{}; + size_t rotatedCapacity{}; + std::map metadata; + + ~Workspace() + { + cudaFree(rotated); + cudaFree(scratch); + cudaFree(blurred); + cudaFree(current); + } +}; + +void +releaseSlices(Workspace & workspace) +{ + cudaFree(workspace.scratch); + cudaFree(workspace.blurred); + cudaFree(workspace.current); + workspace.scratch = nullptr; + workspace.blurred = nullptr; + workspace.current = nullptr; + workspace.sliceCapacity = 0; +} + +bool +allocateBuffer(float ** buffer, size_t elements) +{ + 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); + workspace.rotated = nullptr; + workspace.rotatedCapacity = 0; + if (!allocateBuffer(&workspace.rotated, elements)) + return false; + workspace.rotatedCapacity = elements; + return true; +} + +unsigned int +automaticBatchSize(const Workspace & workspace, int projections, int depth, size_t pixels, bool backward) +{ + 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)))); +} + +unsigned int +allocateBatch(Workspace & workspace, unsigned int batchSize, int depth, size_t pixels, bool backward) +{ + for (;;) + { + 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 +uploadMetadata(DeviceMetadata & device, + HostGaussianMetadata && metadata, + const float * matrices, + const std::vector & inverseMatrices, + const std::vector & firstSlices, + int projections) +{ + 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 +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) +{ + 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 +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 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 +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; + const int localProjection = blockIdx.z; + const int projection = batchStart + localProjection; + if (x >= width || y >= height || z < firstSlices[projection]) + return; + 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 +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 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 +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; + 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 +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, + unsigned int requestedBatchSize, + void ** workspacePointer) +{ + if (!*workspacePointer) + *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, + rotatedSize, + rotatedSpacing, + sigmaZero, + alpha, + rotatedToVolumeMatrices, + farDistances, + firstSlices.data()); + auto [iterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = iterator->second; + if (inserted) + { + HostGaussianMetadata metadata; + std::map> cache; + for (int projection = 0; projection < projections; ++projection) + { + const float nearDistance = nearDistances[projection]; + for (int z = 0; z < depth; ++z) + { + 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 * nearDistance + sigmaZero) * (alpha * nearDistance + sigmaZero); + appendGaussian(metadata, cache, finalVariance, rotatedSpacing[0], rotatedSpacing[1]); + } + uploadMetadata(deviceMetadata, std::move(metadata), rotatedToVolumeMatrices, {}, firstSlices, projections); + } + + 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 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) + { + 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); + } + 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; +} + +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, + unsigned int requestedBatchSize, + void ** workspacePointer) +{ + if (!*workspacePointer) + *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, + rotatedSpacing, + sigmaZero, + alpha, + volumeToRotatedMatrices, + nearDistances, + firstSlices); + auto [iterator, inserted] = workspace.metadata.try_emplace(key); + auto & deviceMetadata = iterator->second; + if (inserted) + { + HostGaussianMetadata metadata; + std::map> cache; + std::vector inverseMatrices(12 * projections); + for (int projection = 0; projection < projections; ++projection) + { + for (int z = 0; z < depth; ++z) + { + 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); + } + uploadMetadata(deviceMetadata, std::move(metadata), volumeToRotatedMatrices, inverseMatrices, first, projections); + } + + 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) + prepareScalarTextureObject(const_cast(volumeSize), + const_cast(devAttenuation), + attenuationArray, + attenuationTexture, + false, + true, + cudaAddressModeClamp); + const dim3 block2(16, 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]); + for (int batchStart = 0; batchStart < projections; batchStart += batchSize) + { + 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) + 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 = 0; z < depth; ++z) + { + storeSlices<<>>( + current, workspace.rotated, pixels, depth, z, deviceMetadata.firstSlices, batchStart); + if (z + 1 == depth) + break; + if (attenuationTexture) + 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); + } + addRotatedBatch<<>>(workspace.rotated, + devVolumeOut, + cudaVolumeSize, + cudaRotatedSize, + deviceMetadata.matrices, + batchStart, + batchCount, + rotatedSpacing[2]); + } + 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()