Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions applications/rtkbackprojections/rtkbackprojections.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
# include "rtkCudaBackProjectionImageFilter.h"
# include "rtkCudaFDKBackProjectionImageFilter.h"
# include "rtkCudaRayCastBackProjectionImageFilter.h"
# include "rtkCudaZengBackProjectionImageFilter.h"
#endif
#include "rtkCyclicDeformationImageFilter.h"

Expand Down Expand Up @@ -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();
Expand Down
2 changes: 1 addition & 1 deletion applications/rtkbackprojections/rtkbackprojections.ggo
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
19 changes: 18 additions & 1 deletion applications/rtkbackprojections/rtkbackprojections.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ def build_parser():
"Joseph",
"JosephAttenuated",
"Zeng",
"CudaZeng",
"CudaFDKBackProjection",
"CudaBackProjection",
"CudaRayCast",
Expand Down Expand Up @@ -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)

Comment on lines +141 to +151

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can be mutualized with Zeng I believe

elif args_info.bp == "CudaFDKBackProjection":
if hasattr(itk, "CudaImage"):
bp = rtk.CudaFDKBackProjectionImageFilter.New()
Expand All @@ -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:
Expand Down
16 changes: 16 additions & 0 deletions applications/rtkforwardprojections/rtkforwardprojections.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "rtkZengForwardProjectionImageFilter.h"
#ifdef RTK_USE_CUDA
# include "rtkCudaForwardProjectionImageFilter.h"
# include "rtkCudaZengForwardProjectionImageFilter.h"
#endif

#include <itkImageFileReader.h>
Expand Down Expand Up @@ -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;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
13 changes: 10 additions & 3 deletions applications/rtkforwardprojections/rtkforwardprojections.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ def build_parser():
"JosephAttenuated",
"CudaRayCast",
"Zeng",
"CudaZeng",
"MIP",
"CudaWrapRayCast",
],
Expand Down Expand Up @@ -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)
Comment on lines +168 to +170

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd remove for consistency with the code below

elif args_info.fp == "MIP":
forwardProjection = rtk.MaximumIntensityProjectionImageFilter[
OutputImageType, OutputImageType
Expand Down Expand Up @@ -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())
)
Expand All @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion applications/rtkiterativefdk/rtkiterativefdk.ggo
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
19 changes: 18 additions & 1 deletion applications/rtkprojectors_group.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -26,6 +26,7 @@ def add_rtkprojectors_group(parser):
"CudaRayCast",
"JosephAttenuated",
"Zeng",
"CudaZeng",
],
default="VoxelBasedBackProjection",
)
Expand Down Expand Up @@ -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)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code duplicates "Zeng", I guess you can just add an or in the Zeng condition above


# Mimicks SetForwardProjectionFromGgo
Expand Down Expand Up @@ -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)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Idem, duplicates Zeng code and can be mutualized.

10 changes: 5 additions & 5 deletions applications/rtkprojectors_section.ggo
Original file line number Diff line number Diff line change
@@ -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
2 changes: 2 additions & 0 deletions documentation/docs/Projectors.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)
81 changes: 81 additions & 0 deletions include/rtkCudaZengBackProjectionImageFilter.h
Original file line number Diff line number Diff line change
@@ -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 <itkCudaImage.h>
# include <itkCudaInPlaceImageFilter.h>
# include <itkImageBase.h>

namespace rtk
{

/** \class CudaZengBackProjectionImageFilter
* \brief CUDA implementation of the rotation-based Zeng backprojector.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add [Zeng et al, IEEE TMI, [Zeng et al, 1999, IEEE TMI, 10.1109/42.796285]

*
* 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<float, 3>,
itk::CudaImage<float, 3>,
BackProjectionImageFilter<itk::CudaImage<float, 3>, itk::CudaImage<float, 3>>>

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I'm not mistaking, a large simplification would be to inherit from ZengBackProjectionImageFilter. Then all common members don't need to be declared.

{
public:
ITK_DISALLOW_COPY_AND_MOVE(CudaZengBackProjectionImageFilter);
using ImageType = itk::CudaImage<float, 3>;
/** Metadata-only image used to build the rotated-grid transforms. Using
* ImageBase avoids the CUDA image factory replacing itk::Image::New(). */
using CPUImageType = itk::ImageBase<3>;
using ProjectorType = BackProjectionImageFilter<ImageType, ImageType>;
using Self = CudaZengBackProjectionImageFilter;
using Superclass = itk::CudaInPlaceImageFilter<ImageType, ImageType, ProjectorType>;
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;

itkNewMacro(Self);
itkOverrideGetNameOfClassMacro(CudaZengBackProjectionImageFilter);

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
Loading
Loading