CUDA Zeng forward and back projectors in RTK by CoConstan · Pull Request #977 · RTKConsortium/RTK

CoConstan · 2026-07-29T10:45:15Z

Add CUDA Zeng forward and back projectors

Summary

This PR adds CUDA implementations of RTK's rotation-based Zeng forward and
back projectors for SPECT reconstruction. The implementation reproduces the
CPU Zeng slice recurrence, including the depth-dependent Gaussian point spread
function (PSF) and optional attenuation correction.

The new projectors operate on itk::CudaImage<float, 3> and keep their main
inputs and outputs on the GPU. They are available from C++, the command-line
applications, iterative reconstruction filters, and Python wrapping.

Main changes

Add rtk::CudaZengForwardProjectionImageFilter.
Add rtk::CudaZengBackProjectionImageFilter.
Add CUDA kernels for:
- volume rotation and trilinear interpolation;
- the Zeng forward and adjoint slice recurrences;
- separable discrete Gaussian PSF convolution;
- optional exponential attenuation.
Match the Gaussian coefficients generated by ITK's discrete Gaussian
implementation, including Bessel-based coefficients, normalization, and
constant-zero boundary behavior.
Keep volumes, projection stacks, attenuation maps, intermediate images, and
outputs GPU-resident between CUDA filters.
Reuse persistent CUDA workspaces and cache geometry/PSF metadata for repeated
iterative reconstruction subsets.
Use shared-memory tiles for the separable Gaussian convolution kernels.
Add CudaZeng to the forward/back projector selectors used by:
- rtkforwardprojections;
- rtkbackprojections;
- iterative reconstruction filters;
- Python applications.
Add Python wrapping for both CUDA filters.

Automated tests

The existing CUDA Zeng tests now cover:

functional forward projection;
output GPU residency;
forward/back adjointness with an active PSF;
adjointness with attenuation;
GPU residency with and without attenuation.

This PR also adds rtkZengProjectionComparisonCudaTest, which performs direct
CPU/CUDA comparisons using separate itk::Image and itk::CudaImage inputs.
It checks:

forward projection with SigmaZero=1.5 and Alpha=0.016;
backprojection with the same PSF;
attenuated forward projection with a uniform attenuation coefficient of
0.01;
attenuated backprojection;
finite output values;
GPU residency before CPU comparison.

The comparison tolerances are:

Comparison	Relative L2 tolerance
Forward projection	`5e-5`
Backprojection	`5e-5`
Attenuated forward projection	`5e-4`
Attenuated backprojection	`2e-3`

Measured values on the development system were:

Comparison	Relative L2 error
Forward projection	`1.31e-5`
Backprojection	`1.30e-5`
Attenuated forward projection	`2.31e-4`
Attenuated backprojection	`1.23e-3`

Compute Sanitizer was also run on the forward projector, backprojector, and a
short OSEM reconstruction. It reported zero memory errors and zero leaked
bytes.

Performance

Application-level measurements were made on an NVIDIA RTX 3500 Ada Laptop GPU
using a 64 x 64 x 64 volume, 60 projections of 64 x 64 pixels, 4 mm
spacing, SigmaZero=1.5, and Alpha=0.016.

Operator	CPU median	CUDA median	Speedup
Zeng forward projection	2.25 s	0.38 s	5.92x
Zeng backprojection	2.01 s	0.62 s	3.24x

An OSEM benchmark with five iterations and ten projections per subset produced
a median runtime of approximately 5.23 seconds. Persistent workspaces and
shared-memory convolution preserve bit-for-bit identical OSEM output compared
with the previous CUDA implementation.

Validation limitation of the local environment

The complete RTKCudaTestDriver cannot currently be linked in the development
environment because the installed older CudaCommon version fails while
compiling the unrelated 4D ROOSTER CUDA test. The new Zeng comparison and
adjointness test translation units compile successfully. A standalone driver
was used to execute the new comparison test, and all four comparisons passed.

This limitation is unrelated to the CUDA Zeng implementation and should not be
present with a current compatible ITK/CudaCommon build.

axel-grc · 2026-08-24T09:32:41Z

Hi @CoConstan , thanks for this nice contribution! I tested your branch and made three small changes locally:

Squashed your first four commits into a single ENH: Add CUDA Zeng projectors commit.
Replaced the static_cast workaround in rtkGgoFunctions.h. The -Wswitch warning comes from rtkiterativefdk, which defines its own fp option (without including rtkprojectors_section.ggo), so its generated enum lacked CudaZeng.
Documentation: added CudaZeng to the CUDA forward/back projector lists in documentation/docs/Projectors.md.

SimonRit

Very impressive work. Can you or @axel-grc (to be discussed) go over the suggestions? I would also indicate in the commit message how was this coded (which AI tool).

SimonRit · 2026-09-10T06:53:33Z

+            recon.SetAlphaPSF(args_info.alphapsf)
+        if args_info.attenuationmap is not None:
+            recon.SetAttenuationMap(attenuation_map)



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

SimonRit · 2026-09-10T06:56:11Z

    )
    rtkprojectors_group.add_argument(
        "--sigmazero",
        help="PSF value at a distance of 0 meter of the detector (Zeng only)",


CudaZeng missing here

SimonRit · 2026-09-10T06:56:18Z

    )
    rtkprojectors_group.add_argument(
        "--attenuationmap",
        help="Attenuation map relative to the volume to perfom the attenuation correction (JosephAttenuated and Zeng)",


CudaZeng missing here

SimonRit · 2026-09-10T06:58:02Z

+        if args_info.alphapsf is not None:
+            recon.SetAlphaPSF(args_info.alphapsf)
+        if args_info.attenuationmap is not None:
+            recon.SetAttenuationMap(attenuation_map)


Idem, duplicates Zeng code and can be mutualized.

SimonRit · 2026-09-10T06:58:30Z

    )
    rtkprojectors_group.add_argument(
        "--alphapsf",
        help="Slope of the PSF against the detector distance (Zeng only)",


CudaZeng missing here

SimonRit · 2026-09-10T07:18:50Z

+    rotatedImage->SetOrigin(origin);
+
+    const MatrixType matrix = GetPhysicalPointToIndexMatrix(this->GetInput(1)).GetVnlMatrix() *
+                              TransformMatrix(transform).GetVnlMatrix() *


Called once, you can just do it here without function

SimonRit · 2026-09-10T07:19:39Z

+}
+} // namespace
+
+CudaZengForwardProjectionImageFilter::CudaZengForwardProjectionImageFilter() { this->InPlaceOff(); }


Why default it to InPlaceOff? I'd keep the regular expected default

SimonRit · 2026-09-10T07:21:02Z

+    if (!transform->GetInverse(inverse))
+      itkGenericExceptionMacro(<< "Could not invert Zeng rotation transform.");
+    const MatrixType matrix = GetPhysicalPointToIndexMatrix(rotatedImage.GetPointer()).GetVnlMatrix() *
+                              TransformMatrix(inverse).GetVnlMatrix() *


Called once, you can just do it here without function

SimonRit · 2026-09-10T07:21:39Z

+    const MatrixType matrix = GetPhysicalPointToIndexMatrix(this->GetInput(1)).GetVnlMatrix() *
+                              TransformMatrix(transform).GetVnlMatrix() *
+                              GetIndexToPhysicalPointMatrix(rotatedImage.GetPointer()).GetVnlMatrix();
+    StoreTextureMatrix(matrix, volumeRegion.GetIndex(), matrices.data() + 12 * local);


Called once, you can just do it here without function

SimonRit · 2026-09-10T07:22:39Z

+}
+} // namespace
+
+CudaZengBackProjectionImageFilter::CudaZengBackProjectionImageFilter() { this->InPlaceOff(); }


Why default it to InPlaceOff? I'd keep the regular expected default

axel-grc reviewed Jul 29, 2026

View reviewed changes

Comment thread src/rtkCudaZengBackProjectionImageFilter.cxx Outdated

CoConstan force-pushed the CudaZeng branch from 37970d6 to 0929cba Compare July 29, 2026 20:11

ENH: Add CUDA Zeng projectors

9070a96

axel-grc force-pushed the CudaZeng branch from 0929cba to 9070a96 Compare August 24, 2026 09:28

axel-grc requested a review from SimonRit August 24, 2026 13:17

SimonRit self-assigned this Aug 26, 2026

SimonRit added this to the RTK 3.0 milestone Aug 26, 2026

SimonRit approved these changes Sep 10, 2026

View reviewed changes

SimonRit requested a review from arobert01 September 10, 2026 08:21

Navigation Menu

Uh oh!

CUDA Zeng forward and back projectors in RTK - #977

CUDA Zeng forward and back projectors in RTK#977
CoConstan wants to merge 1 commit into
RTKConsortium:mainfrom
CoConstan:CudaZeng

CoConstan commented Jul 29, 2026

Uh oh!

Uh oh!

axel-grc commented Aug 24, 2026

Uh oh!

SimonRit left a comment

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

SimonRit Sep 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Uh oh!

Conversation

CoConstan commented Jul 29, 2026

Add CUDA Zeng forward and back projectors

Summary

Main changes

Automated tests

Performance

Validation limitation of the local environment

Uh oh!

Uh oh!

axel-grc commented Aug 24, 2026

Uh oh!

SimonRit left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants