API Reference

This section provides documentation for the differentiable CT operators and analytical reconstruction helpers in diffct. Each autograd class is a PyTorch torch.autograd.Function that runs a CUDA kernel under the hood, enabling gradient propagation through the full projection / backprojection pipeline. The analytical helpers (ramp_filter_1d, cone_cosine_weights, cone_weighted_backproject and friends) are intentionally kept as plain functions so they can be composed with the autograd operators or used directly in one-shot FBP / FDK pipelines.

Overview

The library exposes three families of operators, grouped by geometry:

• Parallel beam (2D): parallel-beam geometry for 2D reconstruction.

• Fan beam (2D): divergent fan-beam geometry.

• Cone beam (3D): full 3D cone-beam geometry for volumetric reconstruction with a 2D flat-panel detector.

Each geometry ships a differentiable forward projector and a differentiable backprojector that are adjoints of each other (the backprojector is the gradient of the projector). For cone-beam the analytical FDK path is routed through a dedicated voxel-driven gather kernel; see the cone_weighted_backproject notes below.

Parallel Beam Operators

The parallel beam geometry assumes parallel X-ray beams, commonly used in synchrotron CT and some medical CT scanners.

class diffct.differentiable.ParallelProjectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 2D parallel beam forward projection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for parallel beam CT geometry. The forward pass computes the sinogram from a 2D image using parallel beam geometry. The backward pass computes gradients using the adjoint backprojection operation. Requires CUDA-capable hardware and a properly configured CUDA environment; all input tensors must reside on the same CUDA device.

Examples

>>> import torch
>>> from diffct.differentiable import ParallelProjectorFunction
>>>
>>> # Create a 2D image with gradient tracking
>>> image = torch.randn(128, 128, device='cuda', requires_grad=True)
>>> # Define projection parameters
>>> angles = torch.linspace(0, torch.pi, 180, device='cuda')
>>> num_detectors = 128
>>> detector_spacing = 1.0
>>> # Compute forward projection
>>> projector = ParallelProjectorFunction.apply
>>> sinogram = projector(image, angles, num_detectors, detector_spacing)
>>> # Compute loss and gradients
>>> loss = sinogram.sum()
>>> loss.backward()
>>> print(f"Gradient shape: {image.grad.shape}")  # (128, 128)

static backward(ctx, grad_sinogram)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, image, angles, num_detectors, detector_spacing=1.0, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0)

Compute the 2D parallel beam forward projection (Radon transform) of an image using CUDA acceleration.

Parameters:

• image (torch.Tensor) – 2D input image tensor of shape (H, W), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_angles,), must be on the same CUDA device as image.

• num_detectors (int) – Number of detector elements in the sinogram (columns).

• detector_spacing (float, optional) – Physical spacing between detector elements (default: 1.0).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as detector_spacing, default: 1.0).

Returns:

sinogram – 2D tensor of shape (num_angles, num_detectors) containing the forward projection (sinogram) on the same device as image.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Uses the Siddon method with interpolation for accurate ray tracing and bilinear interpolation.

Examples

>>> image = torch.randn(128, 128, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, torch.pi, 180, device='cuda')
>>> sinogram = ParallelProjectorFunction.apply(
...     image, angles, 128, 1.0
... )

class diffct.differentiable.ParallelBackprojectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 2D parallel beam backprojection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for parallel beam backprojection. The forward pass computes a 2D reconstruction from sinogram data using parallel beam backprojection, and the backward pass computes gradients via forward projection as the adjoint operation. Requires CUDA-capable hardware and consistent device placements.

Examples

>>> import torch
>>> from diffct.differentiable import ParallelBackprojectorFunction
>>>
>>> sinogram = torch.randn(180, 128, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, torch.pi, 180, device='cuda')
>>> recon = ParallelBackprojectorFunction.apply(sinogram, angles, 1.0, 128, 128)
>>> loss = recon.sum()
>>> loss.backward()
>>> print(sinogram.grad.shape)  # (180, 128)

static backward(ctx, grad_output)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, sinogram, angles, detector_spacing=1.0, H=128, W=128, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0)

Compute the 2D parallel beam backprojection (adjoint Radon transform) of a sinogram using CUDA acceleration.

Parameters:

• sinogram (torch.Tensor) – 2D input sinogram tensor of shape (num_angles, num_detectors), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_angles,), must be on the same CUDA device as sinogram.

• detector_spacing (float, optional) – Physical spacing between detector elements (default: 1.0).

• H (int, optional) – Height of the output reconstruction image (default: 128).

• W (int, optional) – Width of the output reconstruction image (default: 128).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as detector_spacing, default: 1.0).

Returns:

reco – 2D tensor of shape (H, W) containing the reconstructed image on the same device as sinogram.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Uses the Siddon method with interpolation for accurate ray tracing and bilinear interpolation.

Examples

>>> sinogram = torch.randn(180, 128, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, torch.pi, 180, device='cuda')
>>> reco = ParallelBackprojectorFunction.apply(
...     sinogram, angles, 1.0, 128, 128
... )

Fan Beam Operators

Fan beam geometry uses a point X-ray source with a fan-shaped beam, typical in medical CT scanners.

class diffct.differentiable.FanProjectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 2D fan beam forward projection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for fan beam geometry, where rays diverge from a point X-ray source to a linear detector array. The forward pass computes sinograms using divergent beam geometry. The backward pass returns the pure adjoint P^T (Siddon scatter, no distance weighting) so that it is the correct gradient of y = P(x) with respect to x. Analytical fan-beam FBP distance weighting is handled in fan_weighted_backproject, not here.

Examples

>>> import torch
>>> from diffct.differentiable import FanProjectorFunction
>>>
>>> image = torch.randn(256, 256, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2 * torch.pi, 360, device='cuda')
>>> sinogram = FanProjectorFunction.apply(image, angles, 512, 1.0, 1500.0, 1000.0)
>>> loss = sinogram.sum()
>>> loss.backward()
>>> print(image.grad.shape)  # (256, 256)

static backward(ctx, grad_sinogram)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, image, angles, num_detectors, detector_spacing, sdd, sid, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0, backend='siddon')

Compute the 2D fan beam forward projection of an image using CUDA acceleration.

Parameters:

• backend (str, optional) – Forward projector backend. "siddon" (default) is the existing ray-driven Siddon kernel with bilinear interpolation at segment midpoints. "sf" is a prototype separable-footprint (SF-TR) voxel-driven projector that projects each voxel’s footprint as a trapezoid on the detector and integrates it closed-form over each detector cell; it is mass-conserving and closer to the physical finite-width-cell integral, at the cost of being forward-only (no autograd backward yet).

• image (torch.Tensor) – 2D input image tensor of shape (H, W), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_angles,), must be on the same CUDA device as image.

• num_detectors (int) – Number of detector elements in the sinogram (columns).

• detector_spacing (float) – Physical spacing between detector elements.

• sdd (float) – Source-to-Detector Distance (SDD). The total distance from the X-ray source to the detector, passing through the isocenter.

• sid (float) – Source-to-Isocenter Distance (SID). The distance from the X-ray source to the center of rotation (isocenter).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as detector_spacing, sdd, sid, default: 1.0).

Returns:

sinogram – 2D tensor of shape (num_angles, num_detectors) containing the fan beam sinogram on the same device as image.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Fan beam geometry uses divergent rays from a point source to the detector.

• Uses the Siddon method with interpolation for accurate ray tracing and bilinear interpolation.

Examples

>>> image = torch.randn(256, 256, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2 * torch.pi, 360, device='cuda')
>>> sinogram = FanProjectorFunction.apply(
...     image, angles, 512, 1.0, 1500.0, 1000.0
... )

class diffct.differentiable.FanBackprojectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 2D fan beam backprojection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for fan beam backprojection. The forward pass runs the pure adjoint P^T of the fan forward projector (Siddon scatter with bilinear interpolation, no distance weighting), and the backward pass computes gradients via forward projection. The forward/backward pair is therefore a self-consistent autograd adjoint pair. Analytical FBP distance weighting lives in fan_weighted_backproject and is not applied here.

Examples

>>> import torch
>>> from diffct.differentiable import FanBackprojectorFunction
>>>
>>> sinogram = torch.randn(360, 512, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2 * torch.pi, 360, device='cuda')
>>> recon = FanBackprojectorFunction.apply(sinogram, angles, 1.0, 256, 256, 1500.0, 1000.0)
>>> loss = recon.sum()
>>> loss.backward()
>>> print(sinogram.grad.shape)  # (360, 512)

static backward(ctx, grad_output)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, sinogram, angles, detector_spacing, H, W, sdd, sid, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0, backend='siddon')

Compute the 2D fan beam backprojection of a sinogram using CUDA acceleration.

Parameters:

• backend (str, optional) – Adjoint backend, matches FanProjectorFunction. "siddon" (default) uses the ray-driven Siddon scatter kernel, "sf" uses the voxel-driven gather adjoint of the SF-TR projector.

• sinogram (torch.Tensor) – 2D input fan beam sinogram tensor of shape (num_angles, num_detectors), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_angles,), must be on the same CUDA device as sinogram.

• detector_spacing (float) – Physical spacing between detector elements.

• H (int) – Height of the output reconstruction image.

• W (int) – Width of the output reconstruction image.

• sdd (float) – Source-to-Detector Distance (SDD). The total distance from the X-ray source to the detector, passing through the isocenter.

• sid (float) – Source-to-Isocenter Distance (SID). The distance from the X-ray source to the center of rotation (isocenter).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as detector_spacing, sdd, sid, default: 1.0).

Returns:

reco – 2D tensor of shape (H, W) containing the reconstructed image on the same device as sinogram.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Fan beam geometry uses divergent rays from a point source to the detector.

• Uses the Siddon method with interpolation for accurate ray tracing and bilinear interpolation.

Examples

>>> sinogram = torch.randn(360, 512, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2*torch.pi, 360, device='cuda')
>>> reco = FanBackprojectorFunction.apply(
...     sinogram, angles, 1.0, 256, 256, 1000.0, 500.0
... )

Cone Beam Operators

Cone beam geometry extends fan beam to 3D with a cone-shaped X-ray beam and a 2D detector. ConeProjectorFunction and ConeBackprojectorFunction are exact adjoints of each other (Siddon traversal with trilinear interpolation) and drive the iterative cone example. For analytical FDK, use cone_weighted_backproject — it dispatches to a dedicated voxel-driven gather kernel with the classical (sid/U)^2 distance weighting.

Projector backends

Both the fan and cone Projector / Backprojector Function pairs accept a backend keyword that selects the underlying CUDA kernel family:

"siddon" (default)

Ray-driven Siddon traversal with bilinear (2D) / trilinear (3D) interpolation. Fastest, byte-for-byte unchanged from 1.2.11 and earlier, and is the right default for the overwhelming majority of use cases.

"sf" (fan only)

Voxel-driven separable-footprint projector (Long-Fessler-Balter, IEEE TMI 2010). Each voxel’s footprint is a trapezoid built from the four projected corners and is closed-form integrated over each detector cell. Matched forward/adjoint kernels so autograd and the standalone Backprojector both work; use for iterative reconstruction under a mass-conserving cell-integral forward model.

"sf_tr" (cone only)

3D voxel-driven SF with a trapezoidal transaxial footprint and a rectangular axial footprint evaluated at voxel-centre magnification.

"sf_tt" (cone only)

Same transaxial trapezoid as SF-TR but the axial footprint is also a trapezoid built from the four (U_near, U_far) × (z_bot, z_top) projections. Captures axial magnification variation inside a single voxel at large cone angles. Strictly more expressive than SF-TR but 2–4× slower.

All three SF backends go through matched scatter/gather kernel pairs with byte-accurate adjoints (verified by test_adjoint_inner_product and test_gradcheck).

Which backend to pick. At the thin-ray sampling level, measured against the analytical ellipse / ellipsoid Radon transform, Siddon is slightly more accurate than SF (~3 % better RMSE, ~15 % better max error) because diffct’s Siddon uses smooth interpolation rather than nearest-neighbour sampling. However, in a complete analytical reconstruction pipeline (forward projection -> ramp filter -> fan_weighted_backproject / cone_weighted_backproject) SF measurably improves the reconstruction MSE: on the standard 256² examples/fbp_fan.py phantom the raw MSE drops from 0.00216 (Siddon) to 0.00178 (SF), a ~17 % improvement; on the 128³ examples/fdk_cone.py phantom it drops from 0.00325 (Siddon) to 0.00271 (SF-TR), again ~17 %. The mechanism is that SF’s mass-conserving cell-integrated sinogram pairs more naturally with the voxel-driven FBP / FDK gather backprojector than Siddon’s thin-ray sampling, which feeds extra high-frequency ringing into the ramp filter.

So:

• "siddon" — default for backward compatibility and raw forward speed. Best choice when you only need a forward projection and not a reconstruction.

• "sf" / "sf_tr" — recommended for analytical FBP / FDK reconstruction when reconstruction quality matters more than forward-projection runtime. ~2–3× slower forward, ~17 % lower reconstruction MSE on standard CT phantoms.

• "sf_tt" — strictly more expressive than SF-TR but its axial-trapezoid refinement improves RMSE by only ~0.001 at the cost of another ~40 % forward runtime; use for extreme cone angles or algorithm research only.

For iterative reconstruction with autograd, any backend is valid; SF has the theoretical advantage of exact voxel-scale mass conservation on both forward and adjoint, but in practice Siddon’s smoothing helps optimizer convergence.

class diffct.differentiable.ConeProjectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 3D cone beam forward projection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for 3D cone beam geometry. Rays emanate from a point X-ray source to a 2D detector array capturing volumetric projection data. The forward pass computes 3D projections. The backward pass returns the pure adjoint P^T of the forward projector - a Siddon ray-driven scatter with trilinear interpolation weights and no distance weighting - so that P^T is the mathematically correct gradient of y = P(x) with respect to x. Analytical FDK distance weighting is handled separately in cone_weighted_backproject, not here.

Examples

>>> import torch
>>> from diffct.differentiable import ConeProjectorFunction
>>>
>>> volume = torch.randn(128, 128, 128, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2 * torch.pi, 360, device='cuda')
>>> projections = ConeProjectorFunction.apply(volume, angles, 256, 256, 1.0, 1.0, 1500.0, 1000.0)
>>> loss = projections.sum()
>>> loss.backward()
>>> print(volume.grad.shape)  # (128, 128, 128)

static backward(ctx, grad_sinogram)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, volume, angles, det_u, det_v, du, dv, sdd, sid, voxel_spacing=1.0, detector_offset_u=0.0, detector_offset_v=0.0, center_offset_x=0.0, center_offset_y=0.0, center_offset_z=0.0, backend='siddon')

Compute the 3D cone beam forward projection of a volume using CUDA acceleration.

Parameters:

• backend (str, optional) – Forward projector backend. "siddon" (default) is the existing ray-driven trilinear Siddon kernel. "sf_tr" and "sf_tt" are separable-footprint voxel-driven projectors (Long et al., IEEE TMI 2010). SF-TR uses a trapezoidal transaxial footprint and a rectangular axial footprint evaluated at voxel-centre magnification; SF-TT uses trapezoids in BOTH directions, with the axial trapezoid built from U_near and U_far across the voxel to capture axial magnification variation at large cone angles. Both SF backends have matched voxel-driven gather adjoint kernels and support autograd.

• volume (torch.Tensor) – 3D input volume tensor of shape (D, H, W), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_views,), must be on the same CUDA device as volume.

• det_u (int) – Number of detector elements along the u-axis (width).

• det_v (int) – Number of detector elements along the v-axis (height).

• du (float) – Physical spacing between detector elements along the u-axis.

• dv (float) – Physical spacing between detector elements along the v-axis.

• sdd (float) – Source-to-Detector Distance (SDD). The total distance from the X-ray source to the detector, passing through the isocenter.

• sid (float) – Source-to-Isocenter Distance (SID). The distance from the X-ray source to the center of rotation (isocenter).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as du, dv, sdd, sid, default: 1.0).

Returns:

sino – 3D tensor of shape (num_views, det_u, det_v) containing the cone beam projections on the same device as volume.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Cone beam geometry uses a point source and a 2D detector array.

• Uses the Siddon method with interpolation for accurate 3D ray tracing and trilinear interpolation.

Examples

>>> volume = torch.randn(128, 128, 128, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2*torch.pi, 360, device='cuda')
>>> sino = ConeProjectorFunction.apply(
...     volume, angles, 256, 256, 1.0, 1.0, 1500.0, 1000.0
... )

class diffct.differentiable.ConeBackprojectorFunction(*args, **kwargs)

Bases: Function

Summary

PyTorch autograd function for differentiable 3D cone beam backprojection.

Notes

Provides a differentiable interface to the CUDA-accelerated Siddon ray-tracing method with interpolation for 3D cone beam backprojection. The forward pass runs the pure adjoint P^T of the cone forward projector: a Siddon ray-driven scatter with trilinear interpolation weights and no distance weighting. The backward pass computes gradients via 3D cone beam forward projection, which is exactly the adjoint of this pure P^T - so the forward/backward pair is self-consistent for autograd. Analytical FDK distance weighting belongs in cone_weighted_backproject and is not applied here. Requires CUDA-capable hardware and consistent device placements.

This operation may be memory- and computationally-intensive due to 3D geometry. Consider using gradient checkpointing, smaller volumes, or distributed computing for large-scale applications, and ensure sufficient GPU memory is available.

Examples

>>> import torch
>>> from diffct.differentiable import ConeBackprojectorFunction
>>>
>>> projections = torch.randn(360, 256, 256, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2 * torch.pi, 360, device='cuda')
>>> D, H, W = 128, 128, 128
>>> du, dv = 1.0, 1.0
>>> sdd, sid = 1500.0, 1000.0
>>> backprojector = ConeBackprojectorFunction.apply
>>> volume = backprojector(projections, angles, D, H, W, du, dv, sdd, sid)
>>> loss = volume.sum()
>>> loss.backward()
>>> print(f"Projection gradient shape: {projections.grad.shape}")  # (360, 256, 256)

static backward(ctx, grad_output)

Define a formula for differentiating the operation with backward mode automatic differentiation.

This function is to be overridden by all subclasses. (Defining this function is equivalent to defining the vjp function.)

It must accept a context ctx as the first argument, followed by as many outputs as the forward() returned (None will be passed in for non tensor outputs of the forward function), and it should return as many tensors, as there were inputs to forward(). Each argument is the gradient w.r.t the given output, and each returned value should be the gradient w.r.t. the corresponding input. If an input is not a Tensor or is a Tensor not requiring grads, you can just pass None as a gradient for that input.

The context can be used to retrieve tensors saved during the forward pass. It also has an attribute ctx.needs_input_grad as a tuple of booleans representing whether each input needs gradient. E.g., backward() will have ctx.needs_input_grad[0] = True if the first input to forward() needs gradient computed w.r.t. the output.

static forward(ctx, sinogram, angles, D, H, W, du, dv, sdd, sid, voxel_spacing=1.0, detector_offset_u=0.0, detector_offset_v=0.0, center_offset_x=0.0, center_offset_y=0.0, center_offset_z=0.0, backend='siddon')

Compute the 3D cone beam backprojection of a projection sinogram using CUDA acceleration.

Parameters:

• backend (str, optional) – Adjoint backend, matches ConeProjectorFunction. "siddon" (default), "sf_tr", or "sf_tt". See ConeProjectorFunction for details.

• sinogram (torch.Tensor) – 3D input cone beam projection tensor of shape (num_views, det_u, det_v), must be on a CUDA device and of type float32.

• angles (torch.Tensor) – 1D tensor of projection angles in radians, shape (num_views,), must be on the same CUDA device as sinogram.

• D (int) – Depth (z-dimension) of the output reconstruction volume.

• H (int) – Height (y-dimension) of the output reconstruction volume.

• W (int) – Width (x-dimension) of the output reconstruction volume.

• du (float) – Physical spacing between detector elements along the u-axis.

• dv (float) – Physical spacing between detector elements along the v-axis.

• sdd (float) – Source-to-Detector Distance (SDD). The total distance from the X-ray source to the detector, passing through the isocenter.

• sid (float) – Source-to-Isocenter Distance (SID). The distance from the X-ray source to the center of rotation (isocenter).

• voxel_spacing (float, optional) – Physical size of one voxel (in same units as du, dv, sdd, sid, default: 1.0).

Returns:

vol – 3D tensor of shape (D, H, W) containing the reconstructed volume on the same device as sinogram.

Return type:

torch.Tensor

Notes

• All input tensors must be on the same CUDA device.

• The operation is fully differentiable and supports autograd.

• Cone beam geometry uses a point source and a 2D detector array.

• Uses the Siddon method with interpolation for accurate 3D ray tracing and trilinear interpolation.

Examples

>>> projections = torch.randn(360, 256, 256, device='cuda', requires_grad=True)
>>> angles = torch.linspace(0, 2*torch.pi, 360, device='cuda')
>>> vol = ConeBackprojectorFunction.apply(
...     projections, angles, 128, 128, 128, 1.0, 1.0, 1500.0, 1000.0
... )

Analytical Reconstruction Helpers

These helpers build the per-view pre-weights, angle-integration weights, filter, and backprojection pieces of an analytical FBP/FDK pipeline. They are plain functions (no autograd state) and can be freely mixed with the autograd operators above.

diffct.differentiable.detector_coordinates_1d(num_detectors, detector_spacing, detector_offset=0.0, device=None, dtype=torch.float32)

Return centered detector coordinates in physical units.

Coordinates follow the convention (i - (N-1)/2) * spacing + offset. This avoids a half-pixel center bias for even detector counts.

diffct.differentiable.angular_integration_weights(angles, redundant_full_scan=True)

Compute per-view integration weights from the provided angle samples.

Parameters:

• angles (torch.Tensor) – 1D projection angles in radians.

• redundant_full_scan (bool, optional) – If True, applies a 0.5 factor for near-2*pi scans to account for view redundancy in reconstruction formulas using full circular data.

diffct.differentiable.fan_cosine_weights(num_detectors, detector_spacing, sdd, detector_offset=0.0, device=None, dtype=torch.float32)

Return fan-beam cosine pre-weights cos(gamma) for each detector bin.

diffct.differentiable.cone_cosine_weights(det_u, det_v, du, dv, sdd, detector_offset_u=0.0, detector_offset_v=0.0, device=None, dtype=torch.float32)

Return FDK cosine pre-weights D/sqrt(D^2 + u^2 + v^2) on a 2D detector.

diffct.differentiable.parker_weights(angles, num_detectors, detector_spacing, sdd, detector_offset=0.0)

Return Parker redundancy weights for fan/cone short-scan geometries.

For full scans (near 2*pi), this returns all ones.

diffct.differentiable.ramp_filter_1d(sinogram_tensor, dim=-1, sample_spacing=1.0, pad_factor=1, window=None, use_rfft=True)

Apply a 1D ramp filter along dim using FFT.

Parameters:

• sinogram_tensor (torch.Tensor) – Real-valued sinogram tensor.

• dim (int, optional) – Axis along which to filter. Default -1.

• sample_spacing (float, optional) – Physical spacing between detector samples along dim. The continuous ramp filter has units of 1/length, so the FFT output is rescaled by 1/sample_spacing to respect physical units. Default 1.0 (pure sample-unit behavior, matching the historical call signature).

• pad_factor (int, optional) – Zero-pad the signal to pad_factor * N along dim before the FFT to suppress circular-convolution wrap-around artifacts. Default 1 (no padding).

• window (str or None, optional) – Frequency-domain apodization: None/"ram-lak" (unwindowed), "hann", "hamming", "cosine", or "shepp-logan". Default None.

• use_rfft (bool, optional) – Use the real-valued FFT path when True. Default True.

Notes

Pre-existing callers that pass only (sinogram_tensor, dim) keep the same filter shape because the defaults preserve the original behavior up to the rfft path (which is numerically equivalent for real inputs) and the 1/sample_spacing scale (which is exactly 1 when sample_spacing=1.0).

diffct.differentiable.parallel_weighted_backproject(sinogram, angles, detector_spacing, H, W, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0)

Parallel-beam backprojection for analytical FBP pipelines.

Dispatches to the dedicated parallel-beam FBP voxel-driven gather kernel (_parallel_2d_fbp_backproject_kernel). The shared Siddon-based parallel backward kernel continues to serve the autograd adjoint path used by ParallelProjectorFunction / ParallelBackprojectorFunction; this function is only the analytical FBP path. There is no distance weighting in parallel geometry (no source, no magnification), so the kernel simply gathers one linearly-interpolated detector sample per view per pixel and sums. The analytical 1/(2*pi) FBP constant is applied inside the wrapper so the returned image is already amplitude-calibrated - a unit-density disk reconstructs to amplitude 1.

diffct.differentiable.fan_weighted_backproject(sinogram, angles, detector_spacing, H, W, sdd, sid, voxel_spacing=1.0, detector_offset=0.0, center_offset_x=0.0, center_offset_y=0.0, backend='siddon')

Fan-beam weighted backprojection for analytical FBP pipelines.

Dispatches to one of two gather kernels based on backend:

• "siddon" (default): the bilinear voxel-driven fan FBP gather _fan_2d_fbp_backproject_kernel — fastest, bilinearly samples the filtered sinogram at each pixel’s projected u-coordinate.

• "sf": the separable-footprint fan FBP gather _fan_2d_sf_fbp_backproject_kernel — integrates the filtered sinogram over each pixel’s trapezoidal footprint. At nominal voxel size (detector_spacing * sid / sdd) this gives nearly the same MTF as the bilinear gather; at sub-nominal voxel sizes it gives measurably higher spatial resolution, and at supra- nominal it gives higher SNR (see LEAP’s SF vs VD analysis).

Both paths apply the (SID / U)^2 FBP weight per view inside the kernel and the sdd / (2 * pi * sid) analytical FBP scale in Python, so amplitude calibration is preserved across backends. The Siddon-based fan autograd adjoint (_fan_2d_backward_kernel) and the SF-matched adjoint (_fan_2d_sf_backward_kernel) stay unchanged — they continue to serve FanProjectorFunction / FanBackprojectorFunction and are not touched by this analytical wrapper.

diffct.differentiable.cone_weighted_backproject(sinogram, angles, D, H, W, du, dv, sdd, sid, voxel_spacing=1.0, detector_offset_u=0.0, detector_offset_v=0.0, center_offset_x=0.0, center_offset_y=0.0, center_offset_z=0.0, backend='siddon')

Cone-beam weighted backprojection for analytical FDK pipelines.

Dispatches to one of three gather kernels based on backend:

• "siddon" (default): the bilinear voxel-driven FDK gather _cone_3d_fdk_backproject_kernel — fastest, bilinearly samples the filtered sinogram at each voxel’s projected (u, v).

• "sf_tr": the separable-footprint FDK gather _cone_3d_sf_tr_fdk_backproject_kernel — integrates the filtered sinogram over each voxel’s transaxial trapezoidal footprint and axial rectangular footprint at the voxel-centre magnification.

• "sf_tt": the separable-footprint FDK gather _cone_3d_sf_tt_fdk_backproject_kernel — same transaxial trapezoid as "sf_tr" but the axial footprint is also a trapezoid built from four z-corner projections, which more faithfully captures the axial magnification variation inside one voxel at large cone angles.

At nominal voxel size (du * sid / sdd) the three backends give nearly identical MTFs; at sub-nominal voxels the SF variants give measurably higher spatial resolution, and at supra-nominal voxels they give higher SNR (see LEAP’s SF vs VD analysis). All three paths apply the (SID / U)^2 FDK weight per view inside the kernel and the sdd / (2 * pi * sid) analytical FDK scale in Python, so amplitude calibration is preserved across backends. The Siddon-based cone autograd adjoint (_cone_3d_backward_kernel) and the SF-matched adjoints (_cone_3d_sf_tr_backward_kernel / _cone_3d_sf_tt_backward_kernel) stay unchanged — they continue to serve ConeProjectorFunction / ConeBackprojectorFunction and are not touched by this wrapper.

Ramp Filter Options

ramp_filter_1d is a generic 1D ramp filter used by every analytical reconstruction example. Its call signature is:

ramp_filter_1d(sinogram_tensor, dim=-1, sample_spacing=1.0,
               pad_factor=1, window=None, use_rfft=True)

sample_spacing

Physical detector-cell spacing along dim (e.g. du for the cone-beam case). The filter is rescaled by 1 / sample_spacing internally so the output is in physical units and the reconstruction amplitude stays calibrated when detector pitch changes. Pass 1.0 to reproduce the historical sample-unit behavior.

pad_factor

Zero-pad the signal to pad_factor * N samples along dim before the FFT. Common values:

• 1 (default): no padding, fastest but prone to circular convolution wrap-around at the detector edges.

• 2: recommended for cone-beam FDK, good trade-off.

• 4: stricter padding for objects close to the detector edge.

window

Frequency-domain apodization multiplied onto the bare Ram-Lak ramp:

• None or "ram-lak": unwindowed Ram-Lak, sharpest, highest noise.

• "hann": Hann window, strong high-frequency suppression.

• "hamming": slightly milder than Hann.

• "cosine": half-cosine rolloff.

• "shepp-logan": sinc rolloff, classical choice.

use_rfft

Use torch.fft.rfft / irfft when the input is real. Defaults to True; set to False only if you need the complex FFT path.

Analytical FBP / FDK architecture

Each of the three analytical backprojection helpers dispatches to a dedicated voxel-driven gather kernel (parallel / fan / cone), separate from the Siddon-based ray-driven scatter kernels that drive autograd. The autograd kernels are the pure adjoints P^T of the forward projectors (no distance weighting, no magnification scale), so the autograd classes form self-consistent forward/backward adjoint pairs. The analytical helpers on top apply the appropriate FBP/FDK distance weighting and the analytical scale so the returned image is already amplitude-calibrated.

Scale factors applied inside the analytical helpers:

• parallel_weighted_backproject: multiplies by 1 / (2 * pi) (Fourier-convention constant; parallel beam has no source and no magnification).

• fan_weighted_backproject: multiplies by sdd / (2 * pi * sid) (Fourier constant times physical-detector-to-isocenter-plane magnification).

• cone_weighted_backproject: multiplies by sdd / (2 * pi * sid) (same derivation as the fan case; the extra third dimension does not change the ramp-filter convention).

All three helpers use (sid / U)^2 as the FDK/FBP distance weight for the divergent-beam cases (fan and cone), where U = sid + x*sin(phi) - y*cos(phi) is the distance from the source to the pixel/voxel along the central ray direction. Parallel beam has no distance weighting.

Usage Notes

Memory Management

• All operators work on CUDA tensors for optimal performance.

• Ensure sufficient GPU memory for the chosen volume / sinogram size.

• Use torch.cuda.empty_cache() to release cached allocations when switching to a large job.

Gradient Computation

• All autograd operators support automatic differentiation.

• Gradients flow through both forward and adjoint paths.

• Set requires_grad=True on input tensors to enable gradients.

• Analytical helpers (ramp_filter_1d, angular_integration_weights, cone_weighted_backproject, …) are also differentiable in the torch sense (they are pure tensor ops / autograd.Functions), so you can use them inside a loss.

Performance Considerations

• Pass contiguous tensors with the expected dimension order ((D, H, W) for 3D volumes, (num_views, det_u, det_v) for cone sinograms). The library validates layout and will raise if a non-contiguous or transposed tensor is passed.

• Use pad_factor=2 in ramp_filter_1d for cone-beam FDK; the extra FFT cost is usually negligible compared to the backprojection.

Coordinate Systems

• Image / volume coordinates: integer indices, (0, 0, 0) at the corner of the array.

• Detector coordinates: centered at the detector array center, with u in-plane and v axial for cone beam.

• Rotation: counter-clockwise around the z-axis (right-hand rule).