diffct
======

.. raw:: html

   <div class="landing-search">
     <form action="search.html" method="get" role="search">
       <input type="search" name="q" placeholder="Search the documentation" aria-label="Search the documentation" autocomplete="off" autofocus>
     </form>
     <div class="srch-mount"></div>
   </div>

diffct provides differentiable projectors and matched adjoints for computed
tomography (CT) in PyTorch. It supports 2D parallel and fan beams and 3D cone
and parallel beams. Each view can follow any trajectory, and autograd gives
gradients for the volume and for trajectory tensors. Detector pixels can follow
parameterized surfaces, with first-order gradients for their parameters. Views
can be split across several GPUs.

Native kernels run on NVIDIA CUDA GPUs. The portable PyTorch backend runs on
CPU, PyTorch ROCm and Apple Silicon (MPS) devices.

Code written for diffct 1.x must be updated; see :doc:`migration`.

Install
-------

diffct needs Python 3.10 or newer, and a Python version for which PyTorch
publishes wheels. Install PyTorch for your device first. Then print the CUDA
version of your PyTorch build:

.. code-block:: bash

   python -c "import torch; print(torch.version.cuda)"

Install diffct with the extra that matches this version. The extra installs the
CUDA compiler libraries that Numba CUDA needs.

.. tab-set::

   .. tab-item:: 13.x

      .. code-block:: bash

         pip install "diffct[cu13]"

   .. tab-item:: 12.x

      .. code-block:: bash

         pip install "diffct[cu12]"

   .. tab-item:: None (CPU, ROCm or MPS)

      .. code-block:: bash

         pip install diffct

      Use ``backend="torch"``. Importing diffct and using this backend require
      no CUDA runtime.

   .. tab-item:: From source

      .. code-block:: bash

         git clone https://github.com/sypsyp97/diffct.git
         cd diffct
         pip install -e ".[cu13]"    # or ".[cu12]", or "." without CUDA
         python examples/quickstart.py

      Most scripts in ``examples/`` need an NVIDIA GPU; ``plot_trajectory.py``
      and ``benchmark_trajectory_admission.py --backend torch`` run on CPU.

Check the installation:

.. code-block:: bash

   python -c "import torch, diffct; print(diffct.__version__); print(torch.cuda.is_available())"


Backend, precision and solver contracts are in :doc:`execution`.

Quickstart
----------

Create the acquisition once. Then call ``project()`` and ``backproject()``.

.. tab-set::

   .. tab-item:: CUDA

      .. code-block:: python

         import torch
         from diffct import Projector, spiral_trajectory_3d

         trajectory = spiral_trajectory_3d(360, sid=320.0, sdd=512.0, z_range=40.0, n_turns=1.0,
                                           device="cpu")
         A = Projector(trajectory, volume_shape=(128, 128, 128), detector_shape=(384, 256),
                       detector_spacing=0.8)

         volume = torch.rand(A.volume_shape, device="cuda")
         sinogram = A.project(volume)        # (360, 384, 256) line integrals
         adjoint = A.backproject(sinogram)   # matched adjoint A^T

         x = torch.zeros_like(volume, requires_grad=True)
         loss = 0.5 * (A.project(x) - sinogram).square().sum()
         loss.backward()                     # x.grad = A^T (A x - y)

   .. tab-item:: CPU (PyTorch backend)

      .. code-block:: python

         import torch
         from diffct import Projector, spiral_trajectory_3d

         trajectory = spiral_trajectory_3d(60, sid=120.0, sdd=192.0, z_range=15.0, n_turns=1.0,
                                           device="cpu")
         A = Projector(trajectory, volume_shape=(48, 48, 48), detector_shape=(144, 96),
                       detector_spacing=0.8, backend="torch")

         volume = torch.rand(A.volume_shape)
         sinogram = A.project(volume)        # (60, 144, 96) line integrals
         adjoint = A.backproject(sinogram)   # matched adjoint A^T

         x = torch.zeros_like(volume, requires_grad=True)
         loss = 0.5 * (A.project(x) - sinogram).square().sum()
         loss.backward()                     # x.grad = A^T (A x - y)

   .. tab-item:: Apple Silicon (MPS)

      .. code-block:: python

         import torch
         from diffct import Projector, spiral_trajectory_3d

         trajectory = spiral_trajectory_3d(60, sid=120.0, sdd=192.0, z_range=15.0, n_turns=1.0,
                                           device="cpu")
         A = Projector(trajectory, volume_shape=(48, 48, 48), detector_shape=(144, 96),
                       detector_spacing=0.8, backend="torch")

         volume = torch.rand(A.volume_shape, device="mps")
         sinogram = A.project(volume)        # (60, 144, 96) line integrals
         adjoint = A.backproject(sinogram)   # matched adjoint A^T

         x = torch.zeros_like(volume, requires_grad=True)
         loss = 0.5 * (A.project(x) - sinogram).square().sum()
         loss.backward()                     # x.grad = A^T (A x - y)

See :doc:`quickstart` for shapes, devices and troubleshooting, :doc:`examples`
for complete reconstructions, and :doc:`performance` for multi-GPU timings.

Citation and license
--------------------

For software and technical-report citations, see the
`README <https://github.com/sypsyp97/diffct/blob/main/README.md#citation>`_.
The project uses the
`Apache 2.0 license <https://github.com/sypsyp97/diffct/blob/main/LICENSE>`_.
This site documents the source checkout, including
`unreleased APIs <https://github.com/sypsyp97/diffct/blob/main/CHANGELOG.md#unreleased>`_.

.. toctree::
   :hidden:

   getting_started
   user_guide
   examples
   api
   changelog


Getting Started
===============

.. grid:: 1 2 2 2
   :gutter: 3
   :class-container: feature-grid

   .. grid-item-card:: Installation
      :link: installation
      :link-type: doc

      Choose the extra that matches the CUDA version of your PyTorch build, or
      use the PyTorch backend on CPU.

   .. grid-item-card:: Quickstart
      :link: quickstart
      :link-type: doc

      Project an image, apply the matched adjoint and compute a gradient.

.. toctree::
   :hidden:

   installation
   quickstart


Installation
============

Native execution uses an NVIDIA CUDA GPU. The portable PyTorch backend also
supports CPU, PyTorch ROCm and Apple Silicon (MPS) devices; see
:doc:`execution` for precision and hardware-validation limits.

diffct needs Python 3.10 or newer, and a Python version for which PyTorch
publishes wheels. Install PyTorch for your device first. Then print the CUDA
version of your PyTorch build:

.. code-block:: bash

   python -c "import torch; print(torch.version.cuda)"

Install diffct with the extra that matches this version. The extra installs the
CUDA compiler libraries that Numba CUDA needs.

.. tab-set::

   .. tab-item:: 13.x

      .. code-block:: bash

         pip install "diffct[cu13]"

   .. tab-item:: 12.x

      .. code-block:: bash

         pip install "diffct[cu12]"

   .. tab-item:: None (CPU, ROCm or MPS)

      .. code-block:: bash

         pip install diffct

      Use ``backend="torch"``. Importing diffct and using this backend require
      no CUDA runtime.

   .. tab-item:: From source

      .. code-block:: bash

         git clone https://github.com/sypsyp97/diffct.git
         cd diffct
         pip install -e ".[cu13]"    # or ".[cu12]", or "." without CUDA
         python examples/quickstart.py

      Most scripts in ``examples/`` need an NVIDIA GPU; ``plot_trajectory.py``
      and ``benchmark_trajectory_admission.py --backend torch`` run on CPU.

Check the installation:

.. code-block:: bash

   python -c "import torch, diffct; print(diffct.__version__); print(torch.cuda.is_available())"


Numba CUDA initializes on native use. The NVVM and NVJitLink packages must match
the CUDA libraries that PyTorch loads, so the extra must match
``torch.version.cuda``. A plain ``pip install diffct`` also works for native
execution if a matching CUDA Toolkit is already installed.

Troubleshooting
---------------

- **CUDA is unavailable:** check that the active Python environment has a CUDA
  build of PyTorch and that the GPU is visible to the process.
- **Kernel compilation fails:** use Numba CUDA, NVVM and NVJitLink versions that
  match each other. Do not upgrade only one of them.
- **nvJitLink error:** ``nvJitLinkError`` with "may need newer version of
  nvJitLink library" means that the diffct extra does not match the CUDA
  version of PyTorch. Print
  ``torch.version.cuda`` and reinstall with ``diffct[cu12]`` for 12.x or
  ``diffct[cu13]`` for 13.x.


Quickstart
==========

This page projects an image, applies the matched adjoint and computes a gradient
with autograd. Install diffct first; see :doc:`installation`.

.. tab-set::

   .. tab-item:: CUDA

      .. code-block:: python

         import torch
         from diffct import Projector, circular_trajectory_2d_parallel

         trajectory = circular_trajectory_2d_parallel(180, device="cpu")
         A = Projector(trajectory, (64, 64), 96, beam="parallel",
                       detector_spacing=1.0, voxel_spacing=1.0)
         image = torch.zeros(A.volume_shape, device="cuda", dtype=torch.float32)
         image[20:44, 20:44] = 1.0
         measurements = A.project(image)        # (180, 96)
         adjoint = A.backproject(measurements)  # (64, 64), not an inverse

         estimate = torch.zeros_like(image, requires_grad=True)
         loss = 0.5 * (A.project(estimate) - measurements).square().sum()
         loss.backward()                       # A^T (A estimate - measurements)
         print(measurements.shape, estimate.grad.shape)

   .. tab-item:: CPU (PyTorch backend)

      .. code-block:: python

         import torch
         from diffct import Projector, circular_trajectory_2d_parallel

         trajectory = circular_trajectory_2d_parallel(180, device="cpu")
         A = Projector(trajectory, (64, 64), 96, beam="parallel",
                       detector_spacing=1.0, voxel_spacing=1.0, backend="torch")
         image = torch.zeros(A.volume_shape, dtype=torch.float32)
         image[20:44, 20:44] = 1.0
         measurements = A.project(image)        # (180, 96)
         adjoint = A.backproject(measurements)  # (64, 64), not an inverse

         estimate = torch.zeros_like(image, requires_grad=True)
         loss = 0.5 * (A.project(estimate) - measurements).square().sum()
         loss.backward()                       # A^T (A estimate - measurements)
         print(measurements.shape, estimate.grad.shape)

``Projector`` defaults to ``beam="cone"``. Set ``beam="parallel"`` or
``beam="fan"`` for 2D data. Inputs must be floating-point CPU or CUDA tensors with
the configured volume or sinogram shape in their trailing dimensions. Leading
batch/channel dimensions are preserved. Outputs use the selected precision on
the input device. Native CPU data streams through CUDA tiles; the PyTorch
backend computes on its input device. See :doc:`execution` and :doc:`chunking`.
Trajectory tensors can stay on CPU.

Next steps
----------

- A custom or non-circular scan: see :doc:`trajectories`.
- A curved detector: see :ref:`detector-surfaces` and run
  ``python examples/curved_detector.py`` from a source checkout.
- An iterative reconstruction: run ``python examples/iterative_reconstruction.py
  --trajectory helical``. See :doc:`examples`.
- Several GPUs or nodes: see :doc:`multi_gpu`.
- Code written for diffct 1.x: see :doc:`migration`.
- Gradient-based reconstruction with autograd and an optimizer: see
  :doc:`gradient_reconstruction`.

Troubleshooting
---------------

- **Shape error:** use ``(D, H, W)`` for cone volumes and ``(views, U, V)`` for
  cone sinograms. In distributed mode, use ``A.projection_shape`` and
  ``A.view_slice`` rather than the global view count.


User Guide
==========

These pages describe the contracts behind ``Projector`` and the reconstruction
APIs. Start with :doc:`getting_started` if you have not installed diffct yet.

.. grid:: 1 2 2 3
   :gutter: 3
   :class-container: feature-grid

   .. grid-item-card:: Trajectories and geometry
      :link: trajectories
      :link-type: doc

      Per-view trajectory tensors, built-in generators, custom scans, detector
      surfaces and geometry gradients.

   .. grid-item-card:: Execution and backends
      :link: execution
      :link-type: doc

      Backends, precision modes and solver contracts.

   .. grid-item-card:: Memory and chunking
      :link: chunking
      :link-type: doc

      Memory limits, spatial, view and pixel tiles, block stores and
      disk-backed volumes.

   .. grid-item-card:: Multi-GPU and multi-node
      :link: multi_gpu
      :link-type: doc

      Several GPUs in one process, or ``torchrun`` across nodes, with view,
      space or hybrid partitions.

   .. grid-item-card:: Performance
      :link: performance
      :link-type: doc

      Forward, adjoint and CGLS timings for a 512³ volume on 1, 4 and 8 GPUs.

   .. grid-item-card:: Differentiable FBP and FDK
      :link: analytical_learning
      :link-type: doc

      Analytical reconstruction with autograd, and learnable filters and
      weights.

   .. grid-item-card:: Migrating from 1.x
      :link: migration
      :link-type: doc

      Move code from the circular-orbit API of diffct 1.x to per-view
      trajectories.

.. toctree::
   :hidden:
   :maxdepth: 2

   trajectories
   execution
   chunking
   multi_gpu
   performance
   analytical_learning
   migration


Trajectories and Geometry
=========================

Trajectory format
-----------------

``Projector`` takes one row per view in each trajectory tensor. The beam type
selects the tuple layout. The helpers in ``diffct`` create these tensors. You can
also pass measured or calibrated poses in the same layout.

.. list-table::
   :header-rows: 1

   * - Beam
     - Tuple (each component's shape)
     - Volume / sinogram
   * - ``parallel``
     - ``(ray_dir, det_origin, det_u)``, each ``(N, 2)``
     - ``(H, W)`` / ``(N, U)``
   * - ``fan``
     - ``(src_pos, det_center, det_u)``, each ``(N, 2)``
     - ``(H, W)`` / ``(N, U)``
   * - ``cone`` (default)
     - ``(src_pos, det_center, det_u, det_v)``, each ``(N, 3)``
     - ``(D, H, W)`` / ``(N, U, V)``
   * - ``parallel3d``
     - ``(ray_dir, det_center, det_u, det_v)``, each ``(N, 3)``
     - ``(D, H, W)`` / ``(N, U, V)``

Coordinates are Cartesian ``(x, y[, z])``. Volume indexing is ``[y, x]`` for 2D
and ``[z, y, x]`` for 3D. The origin is the volume centre. Use one length unit
for positions, voxel spacing and detector spacing. Voxel spacing is a single
positive scalar or a positive per-axis tuple in tensor order: ``(dy,dx)`` or
``(dz,dy,dx)``. Both backends support per-axis spacing; see :doc:`execution`.

Detector axes are unit direction vectors, not pixel-step vectors. The detector
spacing sets the pixel pitch. A scalar spacing applies to both 3D detector axes.
With the default flat detector, pixel ``(u, v)`` is at:

.. code-block:: text

   det_center + (u - (U - 1) / 2) * du * det_u
              + (v - (V - 1) / 2) * dv * det_v

A measured sinogram in ``(views, V, U)`` order must be transposed to
``(views, U, V)`` before use. Check the physical direction of each detector axis.
Reshaping does not correct axis orientation.

Custom acquisition example
--------------------------

This example modifies a generated helix with a per-view detector shift. For a
measured acquisition, replace the four tensors with your calibrated arrays in the
coordinate system above.

.. code-block:: python

   import torch
   from diffct import Projector, spiral_trajectory_3d

   trajectory = spiral_trajectory_3d(
       90, sid=80.0, sdd=128.0, z_range=10.0, n_turns=1.0, device="cpu")
   source, center, u_axis, v_axis = trajectory
   phase = torch.arange(90, dtype=torch.float32) * (2 * torch.pi / 90)
   center = center + (0.25 * phase.sin())[:, None] * u_axis
   A = Projector((source, center, u_axis, v_axis), (32, 32, 32), (96, 64),
                 detector_spacing=(0.8, 0.8))
   volume = torch.ones(A.volume_shape, device="cuda")
   sinogram = A.project(volume)  # (90, 96, 64)

All components must be finite floating-point tensors with the same nonzero view
count. ``ray_dir`` and ``det_u`` must be orthogonal for 2D parallel beams. ``det_u``
and ``det_v`` must be orthogonal for 3D beams. ``parallel3d`` permits oblique
detectors: the unit ray direction must not be tangent to the detector plane.
Degeneracy is checked by the frame's numerical rank in float64.
Its rays are infinite lines through each pixel. It uses the PyTorch backend
on CPU or CUDA; ``backend="auto"`` selects it automatically and explicit
``backend="numba"`` is unsupported. With flat fan/cone detectors,
the source must not coincide with the detector centre or view the detector
edge-on. The native FP32 backend requires the nearer of source and detector
centre to lie within 1e6 voxel spacings of the origin. Native surface detectors validate each actual
source/pixel pair: the endpoints must remain distinct in float32, and the
nearer endpoint of every ray must lie within 1e6 voxel spacings of the origin.
Both native endpoints must lie within 1e15 voxel spacings. Portable precision
modes validate finite, nondegenerate rays in their selected dtype; see :doc:`execution`.

Arbitrary trajectories do not guarantee enough angular coverage for an inverse
problem. Use iterative methods with regularization for sparse or incomplete
acquisitions. The analytical helpers assume a circular orbit (see :doc:`api`).

.. _detector-surfaces:

Parameterized detector surfaces
-------------------------------

``Projector(..., detector_surface=surface)`` places each pixel
on a native curve. The legacy two-argument callback places pixels
on an arc, cylinder or another parameterized surface. For 3D beams it receives
physical float64 ``u, v`` grids on CPU with shape ``(U, V)`` and ``ij`` indexing. It returns local
``(u, v, n)`` offsets shaped ``(U, V, 3)`` or ``(views, U, V, 3)``. Here ``views``
is the total trajectory view count, also in distributed mode. The world point is:

.. code-block:: text

   center + offset_u * det_u + offset_v * det_v
          + offset_n * cross(det_u, det_v)

For 2D beams, ``u, v`` have shape ``(U,)`` and ``v`` is zero. Return ``(U, 3)``
or ``(views, U, 3)`` with zero middle offsets; the normal is
``(det_u_y, -det_u_x)``. A cylindrical surface can be written as:

.. code-block:: python

   radius = torch.tensor(128.0, requires_grad=True)

   def cylinder(u, v):
       u, v = u.to(radius.device), v.to(radius.device)
       angle = u / radius
       return torch.stack((radius * angle.sin(), v,
                           radius * (angle.cos() - 1)), dim=-1)

   C = Projector(trajectory, (32, 32, 32), (96, 64),
                 detector_spacing=(0.8, 1.0), detector_surface=cylinder)

For this generated frame the normal points away from the source; negative
normal offsets and a radius equal to ``sdd`` centre the cylinder on the source.

The callback is sampled at construction and afresh on each operation, and its
captured parameters receive first-order gradients. Every sample must be a finite floating-point tensor of
the documented shape; fan/cone pixels cannot coincide with the source.
``project`` and ``backproject`` retain their matched Siddon model, with one ray
per pixel and no detector-area integration. They operate directly on the native
curved grid without resampling. Run ``python examples/curved_detector.py`` for a
complete circular cone example, also shown in :doc:`examples`.

Full-volume and streamed operations accept CPU or CUDA offsets. The operator
transfers selected offsets to the compute device before expanding world points
for that device's view shard or streamed batch. Native execution preserves dtype
promotion during composition before staging complete coordinates as float32.
The PyTorch backend uses the selected precision. See :doc:`chunking` for memory
requirements.

For bounded view/pixel sampling, pass explicit tensor parameters:

.. code-block:: python

   from diffct import ParameterizedSurface

   def sample(u, v, view_indices, r):
       return torch.stack((r * torch.sin(u / r), v,
                           r * (torch.cos(u / r) - 1)), dim=-1)

   surface = ParameterizedSurface(sample, parameters=(radius,))
   C = Projector(trajectory, (32, 32, 32), (96, 64),
                 detector_spacing=(0.8, 1.0), detector_surface=surface,
                 view_chunk_size=8, detector_chunk_shape=(32, 32))

``u,v`` are a rectangle of global physical pixel centres; ``view_indices``
contains global view IDs. Return ``(*pixel_shape,3)`` shared offsets or
``(len(view_indices),*pixel_shape,3)`` batched offsets. In 2D the middle
component remains zero. Streamed grids, IDs and explicit parameter snapshots
are CPU tensors; parameter dtypes and original gradient devices are preserved.
Native sampled-surface arithmetic promotes frames and offsets to one common
dtype, at least float32; float64 inputs retain float64 arithmetic. The PyTorch
backend uses the selected precision. Gradients return in each input's original
dtype. Offsets and world points are validated
on every executing batch.

The sampler is pure. Mutable geometry tensors must be explicit parameters;
other captured configuration must be immutable. Coupling across views uses
global IDs or explicit parameters, including values for omitted views. Each
forward captures the sampler and one set of frame/parameter values. Backward
recomputes local geometry from those snapshots and reduces pixel cotangents to
parameter gradients before communication. It retains first-order geometry
gradients and data Hessians. Explicit parameter snapshots scale with parameter
count; allocations inside user samplers remain caller-controlled. Do not cache
mutable callback values across calls.

Analytical FBP/FDK and weighting helpers continue to assume flat detectors and
do not accept the callback. Reusing a flat-detector FDK requires matching-ray
resampling of the line integrals onto a covered virtual flat grid; this adds
interpolation error and is outside the core operator.

Geometry gradients
------------------

Gradients with respect to trajectory tensors are first order only. Set
``requires_grad=True`` on each tensor before you construct ``Projector``. The
projector keeps references to all trajectory components when any component
requires gradients, and reads and validates their current values on each call.
If none require gradients, it clones the entire trajectory. For example:

.. code-block:: python

   # Continue from the custom acquisition above.
   target = sinogram.detach()
   learnable_center = center.clone().requires_grad_()
   B = Projector((source, learnable_center, u_axis, v_axis), (32, 32, 32),
                 (96, 64), detector_spacing=0.8)
   loss = (B.project(volume) - target).square().sum()
   loss.backward()
   print(learnable_center.grad.shape)  # (90, 3)

This loss uses the same target geometry, so the gradient is zero. For a real
calibration loop, see ``examples/geometry_calibration.py``.

If poses come from learnable angles or offsets, rebuild the derived tensors
inside each optimization step and pass them to the same projector:

.. code-block:: python

   from diffct import circular_trajectory_3d

   angles = torch.linspace(0, 2 * torch.pi, 91, device=volume.device)[:-1]
   angles = angles.detach().requires_grad_()
   optimizer = torch.optim.Adam([angles], lr=1e-3)
   for _ in range(10):
       optimizer.zero_grad()
       current = circular_trajectory_3d(angles, sid=80.0, sdd=128.0)
       loss = (A.project(volume, trajectory=current) - target).square().sum()
       loss.backward()
       optimizer.step()

``project``, ``backproject``, ``__call__``, ``project_into`` and
``backproject_into`` accept the keyword-only ``trajectory`` argument. It replaces
the complete tuple for that call and must keep the configured beam and view
count. Omission or ``None`` uses the constructor trajectory. Each calculation
keeps its own geometry snapshot for backward; overrides leave the default
trajectory unchanged. The ``*_into`` methods retain their numerical,
no-autograd contract. Distributed ranks must agree on whether an override is
present and supply one coherent global trajectory.

Call-time admission checks flat-detector corners and coordinate bounds, with
exact pixel checks for ambiguous views. Frames stage once per validation device.
Native curved-detector admission uses bounded batches on each compute device
after view sharding. Distributed admission errors reach all ranks before
numerical execution.

Parameterize rotations so that axes stay unit length and orthogonal. Frame
checks run at construction and on each call with live default frames or an override; sampled surface positions
are checked on every operation.

Differentiable geometry helpers
-------------------------------

``circular_trajectory_2d_parallel``, ``circular_trajectory_2d_fan``,
``circular_trajectory_3d`` and ``circular_trajectory_3d_parallel`` accept either
a view count or a nonempty floating-point angle tensor as the first argument.
Angles are in radians and may be nonuniform or unsorted. A tensor preserves its
device, dtype and autograd history. ``start_angle``, ``end_angle``, ``device``
and ``dtype`` only configure count-based generation. Parallel helpers cover
``[0, pi)`` by default; fan and cone helpers cover ``[0, 2*pi)``.
The 3D parallel helper lifts the 2D orbit into
the XY plane with ``det_v=(0,0,1)``. In both parallel helpers,
``detector_distance`` is a transverse shift along ``det_u``.

``transform_trajectory(trajectory, rotation=None, translation=None, *, beam="cone")``
applies a rigid transform about the world origin. Points become ``R @ p + t``;
directions become ``R @ d``. ``rotation`` is a proper orthogonal matrix shaped
``(d,d)`` or ``(views,d,d)``, and ``translation`` is shaped ``(d,)`` or
``(views,d)``. Use the same dtype and device as the trajectory. Both shared and
per-view transforms preserve gradients to the frames and transform parameters.
The helper validates frames and transforms; ``Projector`` checks actual
detector pixels. Surface offsets may define valid rays even when the reference
detector center coincides with the source.

.. code-block:: python

   from diffct import transform_trajectory

   shift = torch.zeros(3, device=volume.device, requires_grad=True)
   current = transform_trajectory(tuple(t.detach() for t in current), translation=shift)
   shifted = A.project(volume, trajectory=current)

Limits:

- Second derivatives with respect to geometry raise an error.
- Forward-mode derivatives (``torch.func.jvp``, ``jacfwd``) with respect to trajectory
  tensors raise an error. Use reverse mode for geometry gradients.
- Second derivatives with respect to the volume or sinogram are available,
  including when trajectory or surface parameters require gradients.
- Detector and voxel spacing are fixed physical settings, not differentiable parameters.
- Geometry derivatives are piecewise: at exact voxel edges or corners they are one-sided.
- Distributed geometry gradients are summed by the operator. See :doc:`multi_gpu`
  for loss scaling.


```{include} ../EXECUTION.md
```


Large Volumes and Chunked Execution
===================================

This page describes the native CUDA path. For portable execution and analytical
streaming, see :doc:`execution`.

Pass CPU volumes and sinograms to ``Projector`` when complete arrays exceed GPU
memory. Projection and the matched adjoint execute on CUDA; outputs and data
gradients return to the input device. Tensor methods still hold complete arrays
on that device. For data exceeding host RAM, use the caller-owned block-store
methods below and the disk-backed CGLS example.

.. code-block:: python

   import torch
   from diffct import Projector, circular_trajectory_3d

   trajectory = circular_trajectory_3d(60, sid=160, sdd=256, device="cpu")
   A = Projector(trajectory, (64, 64, 64), (192, 128), detector_spacing=0.8)
   x = torch.ones(A.volume_shape, device="cpu", requires_grad=True)
   y = A.project(x)                  # CPU tensor, real CUDA computation
   y.square().mean().backward()      # x.grad stays on CPU
   adjoint = A.backproject(y.detach())

Automatic selection and manual limits
-------------------------------------

``volume_chunk_shape=None`` is the automatic default. CPU data always streams.
For CUDA data, the original full-volume path is retained when its estimated
additional working set fits each selected card and the output device. Otherwise
the operation streams spatial tiles and view batches.

The arithmetic estimate accounts for float32 staging, cone layout/output
copies, view data, frame or sampled pixel geometry, VJP buffers and collective
scratch. It reserves 25% of available memory for runtime overhead, respects
PyTorch's per-process allocator ceiling and sizes against the smallest
participating budget. Concurrent allocations can still cause allocation errors.
Native PyTorch admission checks allocation-segment capacity separately from the
live tensor payload estimate. A free split inside a partially occupied cached
segment is not a promise that another buffer can use it. Small-pool and
large-pool segment rounding also contribute to admission; a fitting sum of
tensor sizes alone is insufficient. This check is conservative and may select
smaller automatic tiles. Other allocator backends do not inherit native segment
sizes. The segment estimate assumes standard native allocator settings and
fresh fixed-buffer inventories; arbitrary fragmentation or allocator overrides
can still cause allocation errors.
On older PyTorch versions without the public memory-fraction getter, the
estimate can observe driver/allocator memory only; use explicit limits when a
process ceiling is set there.

To override the automatic limits:

.. code-block:: python

   A = Projector(trajectory, (64, 64, 64), (192, 128), detector_spacing=0.8,
                 volume_chunk_shape=(16, 32, 32), view_chunk_size=8,
                 detector_chunk_shape=(64, 64))

Spatial limits follow tensor order: ``(H, W)`` for parallel/fan and ``(D, H, W)``
for cone beams. Limits are positive integral dimensions, clipped at the volume
edge. ``view_chunk_size`` may also be supplied alone to force view streaming.
Its automatic value is at most 32 and can shrink when geometry/data buffers
require it. ``detector_chunk_shape`` follows ``(U,)`` or ``(U,V)`` order.
Automatic execution can split a single detector image when it does not fit;
reducing the volume alone cannot solve that case. Pixel rectangles keep the
global detector origin. Explicit limits must fit the planned working set.

Tiles retain the global origin-centred voxel lattice. Projection adds disjoint
tile integrals; the adjoint accumulates each tile over view batches. Odd sizes,
partial tails and internal boundaries retain the same cell ownership. Chunking
preserves the cell-constant Siddon model and matched adjoint, with float32
rounding differences from a different summation order. Smaller tiles may cost
more transfers and kernel launches.

Resident layouts and accumulation
---------------------------------

For backprojection, each active GPU keeps one contiguous native-layout
accumulator for a spatial tile. Every view batch adds to that accumulator;
layout conversion and return to the caller happen after all views contribute.
The accumulator is cleared once per tile. Projection prepares the native volume
layout once per tile/device and reuses its ray output buffer across view
batches. Geometry VJPs also prepare each tile once across their view batches. Cone kernels
retain contiguous ``(W, H, D)`` access; tensor volumes remain ``(D, H, W)``.
Forward kernels overwrite each active ray, including missed rays and short
final batches, so they do not need unconditional output clearing.

``native_volume_layout=True`` optionally retains native full-CUDA cone output
strides ``(1, D, D * H)`` for an unbatched flat-detector backprojection. The
shape stays ``(D, H, W)``. Projection accepts that layout without conversion,
so one compute device can avoid both volume layout copies in ``A^T A``.
The default is ``False``. Several local devices, surfaces, batches, streaming
and the PyTorch backend still involve copies. Distributed ranks must agree
on this setting.

For example, a ``1024**3`` float32 volume is 4 GiB. With 720 views in batches of
32, returning a complete partial volume after each of the 23 batches would
transfer 92 GiB. Returning the completed accumulated volume transfers 4 GiB
on a single GPU, excluding projections, geometry and other solver data. This
is a copy-payload estimate; it does not imply a 23-fold reconstruction speedup.
Float32 summation order can change slightly, while the underlying Siddon cell
model and matched adjoint stay the same.

Scheduling and transfers
------------------------

``schedule="auto"`` compares a few spatial-first, view-first and bounded-window
orders, including contiguous z slabs and three-dimensional blocks. Selection
uses bounded real CUDA pilot times scaled by the full native launch count,
with estimated full-workload copy payloads as a tie-breaker. The pilots use
small residency/order prefixes, so selection is a heuristic; reduced copy
payloads alone do not establish a runtime improvement. The benchmark example
preserves the initial selection report separately from its warmed timings.
``schedule="spatial"``, ``"views"`` or ``"window"`` enables a reproducible
comparison. Backprojection always completes each retained volume accumulator
before returning it. View-first backprojection is eligible only when all its
accumulators fit together.

The transfer pipeline uses at most two reusable pinned slots per device,
separate upload, compute and download streams, and CUDA events. Events protect
input consumption, host output reads and reuse of GPU buffers. Each composite
host slot is capped at 8 MiB, with a separate 8 MiB pageable read/cast scratch
ceiling. These are memory policy ceilings. Larger resident GPU tiles use
bounded host fragments; fragmentation changes copy call counts, while the
completed volume payload is transferred once. Planned local surface geometry
also has an independent host working-set ceiling. Arbitrary allocations inside
a user sampler and OS file-cache residency remain outside these ceilings.
All additional GPU slots, layouts, geometry and collective staging contribute
to the GPU budget. Very small budgets can use a bounded serial fallback.
Single-slot plans retire prior native computation before allocating the next
tile. Compute and internal allocations use the caller stream. Repeated calls
on that stream can reuse its cache, so temporary transfer streams do not divide
the native allocator's reusable cache. A new buffer is fenced against earlier caller work before
upload. This can delay prefetch of a newly allocated volume buffer; existing
ray/geometry slots retain their bounded overlap.
Transfers hold source and destination references in the bounded completion-event
queue until they finish. Completed internal copies release these references
directly, avoiding stale transfer-stream reservations across calls.

``A.last_execution_stats`` reports the last instrumented operation's schedule,
limits, layout/kernel counts and actual upload/download/peer-copy payloads.
``estimated_gpu_bytes`` describes the tensor-buffer inventory;
``estimated_reserved_bytes`` and ``reservation_budget_bytes`` expose the
separate allocation-segment estimate and available room. Copy counts describe
transferred numerical tensor bytes, independently of segment sizes.
Distributed layout metadata and planning control messages are excluded.
Pilot timing and payloads are reported separately; candidate full-workload
costs are estimates. These counts are not physical bus measurements. The value
is ``None`` before execution or for an uninstrumented full-CUDA path. Run
``python examples/benchmark_execution.py`` to compare schedules and slab/block
shapes on your hardware.

Caller-owned blocks and disk files
------------------------------------------------------------

``TensorStore`` provides detached numerical block access to an existing tensor.
``NpyStore`` maps a float32/float64 ``.npy`` file. Creating a store reserves a new
path exclusively; existing files are not overwritten. Block reads/casts and
writes do not materialize a complete mapped array.

.. code-block:: python

   from diffct import NpyStore

   # The file has shape A.volume_shape, in tensor axis order.
   source = NpyStore("volume.npy", mode="r")
   rays = NpyStore.create("projections.npy", A.projection_shape)
   A.project_into(source, rays)
   rays.flush()
   result = NpyStore.create("adjoint.npy", A.volume_shape)
   A.backproject_into(rays, result)
   result.flush()

``project_into`` and ``backproject_into`` return the supplied output. They
accept tensors, these stores, or structural backends with ``shape``, ``dtype``,
``read(slices)``, ``write(slices,value,accumulate=False)`` and ``flush()``. Outputs must
match ``A.dtype`` and be writable. Native fast output is float32; explicit
``precision="float64"`` output is float64. Each completed output region overwrites its prior
contents; no full result is allocated internally. These numerical methods do
not build an autograd graph. Use tensor ``project``/``backproject`` for
differentiation.

Custom CUDA stores must produce returned blocks on the source device's
current stream, or make that stream wait for their producer. The upload
stream waits for this work before reading a block, including on the same GPU.

Input and output must have distinct backing data. Identical objects,
detectable overlapping tensors and the same canonical mapped file path are
rejected before execution. Custom backends must guarantee that hidden backing storage is
distinct; an opaque object's path attribute alone is not an alias guarantee.
Calling a store's ``read`` with a complete array slice can still request a
complete array; keep application-side reads and vector updates blockwise too.

Autograd, curved detectors and several GPUs
-------------------------------------------

Data gradients and data Hessians, plus first-order frame/surface gradients, use
the streamed path too. Saved CPU geometry keeps backward consistent with its
forward call if parameters change. Geometry second derivatives remain
unsupported.

For bounded geometry, use ``ParameterizedSurface(sampler, parameters=...)``.
Its sampler receives global physical pixel centres, global view IDs and
explicit parameter snapshots, and generates only the executing view/pixel
batch. Pixel cotangents are reduced through the local sampler/frame graph to
parameter gradients before communication. The captured sampler and one set of
frame/parameter snapshots preserve backward behavior after parameter mutation.
See :ref:`detector-surfaces`.

Legacy two-argument callbacks receive complete CPU detector grids and can return
CPU or CUDA offsets. The operator transfers selected offsets to each compute
device before expanding world points for its owned views or streamed batch.
Memory planning uses these local view counts and the composition dtype.
Gradient-connected offsets are snapshotted on CPU and transferred per batch.
A full single-process CUDA call instead keeps admitted CUDA snapshots on the
input GPU.
Streamed composition stages a complete admitted view/pixel batch and composes
it with one launch into reusable float32 geometry buffers. If that workspace
does not fit, a smaller fallback workspace composes one view and a bounded
detector-row segment at a time. Geometry backward replays one view and a
bounded detector-row segment at a time, without retaining complete CUDA
composition graphs.
Composition and replay buffers never coexist. Admission charges the larger
phase total and the largest allocation of either phase; the segment estimate
counts the buffers of both phases, because one phase can leave segments pinned.
Native admission starts with the owned view shard and detector rectangle,
limited by explicit hints. It reduces validation batches to fit device memory,
allocation and host-staging limits. These batches are independent of execution
chunks and their automatic 32-view cap.
Callback allocations remain caller-controlled. Callback values are not cached
across operations.

``devices=[0, 1]`` still partitions views in list order. Each card streams tiles;
the full output returns to the input device. Distributed ``partition="views"``
ranks retain local view shards and replicated volumes. Local GPU backprojection
sums and NCCL reductions occur while the tile is resident, before its final
host transfer. Gloo uses bounded CPU staging. All ranks need matching settings
and must participate in the same calls, including empty ranks. The
ownership-specific loss rules in :doc:`multi_gpu` apply. Scalar NCCL reductions
also use CUDA tensors.

With CPU streaming, the cards compute their view batches concurrently. Each
card receives every volume tile, so host transfers limit the multi-GPU speedup.

With ``partition="views"``, every rank holds the complete tensor volume and
iterative state. One process per node controlling ``devices=[...]`` avoids
extra host replicas. ``partition="space"`` instead assigns balanced contiguous
first-axis slabs to ranks. ``global_volume_shape`` retains the constructor
shape; ``volume_slice`` gives global ownership and ``volume_shape``/
``local_volume_shape`` give the required local tensor/store shape. Projection
SUMs the slab contributions into replicated ray batches; backprojection
returns only the owned slab. Physical coordinates remain on the global voxel
lattice. Empty slabs are valid. The PyTorch backend also supports
``partition="hybrid", partition_grid=(view_ranks, space_ranks)``; see
:doc:`execution` for ownership and loss divisors.

Caller-owned CUDA inputs and complete CUDA outputs still occupy their full
size. Keep arrays, trajectory tensors and learnable surface parameters on CPU
for bounded GPU residency. CUDA geometry snapshots or CUDA intermediates made
inside a callback also consume GPU memory outside the streamed tile buffers.
``fbp``/``fdk`` also support view/spatial chunks and numerical stores; see
:doc:`execution`. The helical approximation and low-level CUDA helpers retain
their original execution path.

Runnable iterative reconstruction
---------------------------------

The example reuses the existing CGLS solver with CPU arrays. It supports one
process with one or several GPUs and prints measured residuals and peak
PyTorch CUDA tensor allocation. The peak excludes driver/context memory.

.. code-block:: bash

   python examples/chunked_reconstruction.py
   python examples/chunked_reconstruction.py --chunk-shape 16 32 32 --view-chunk-size 8
   python examples/chunked_reconstruction.py --devices 0 1

The tensor CGLS helper computes FP64 norms in bounded pieces and reuses vector
update buffers. Its complete tensor states still need memory. If any rank's
CUDA state exceeds its budget, all ranks agree to keep state on CPU and
replan together, preserving collective order for unequal slabs. Disk-backed
CGLS keeps ``x,s,p`` plus ``r,q`` in stores, performs blockwise vector updates
and writes completed-iteration checkpoints without gathering the volume.
``q`` remains available until the global norm determines ``alpha`` and the
``x``/residual updates finish. Account for all three volume states together;
fitting one volume does not establish that the complete solver fits.
CUDA admission also checks each state and FP64 norm block against the feasible
individual allocation limit, together with the forward and adjoint pipeline
inventories. A small total state size cannot bypass an allocation-segment limit.
State allocations can pin complete native allocator segments. Admission reserves
these segments separately from the later pipeline or FP64 norm work, matching
the fresh-segment budget used when execution replans after state allocation.
It does not assume later calls can use free splits inside the state segments.
The state bound also includes compatible idle segments on the caller's stream:
reusing an oversized cached segment can pin more memory than a fresh allocation.
The example helpers set ``operator.last_solver_stats`` with state and scratch
capacity metadata. This report belongs to the example helpers. The library
also provides numerical CGLS/SIRT/TV solvers; see :doc:`execution`.
The disk entrypoint saves it alongside the actual residual, elapsed time and
CUDA tensor peak.

.. code-block:: bash

   python examples/disk_reconstruction.py --output disk-run
   python examples/disk_reconstruction.py --surface saddle --output saddle-run
   torchrun --standalone --nproc-per-node=2 examples/disk_reconstruction.py \
       --partition space --output spatial-run

The default disk example generates explicitly synthetic measurements
blockwise. Each rank owns its output, workspace and checkpoint files. Resume
requires matching partition, ownership, shapes, scan and measurements. Files
back the arrays; the OS may cache their pages in RAM. Local tests and bounded
models provide evidence separately from physical multi-GPU or cross-node runs.
The iteration count is the target total: resuming iteration 2 with
``--iterations 5`` performs three further iterations in a new run directory.

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/chunked_reconstruction.py
   :language: python
   :linenos:


Multiple GPUs and Nodes
=======================

By default, ``Projector`` splits views across selected GPUs and automatically streams native spatial tiles
when arrays are CPU-backed or the estimated CUDA working set does not fit.
For volumes larger than one GPU, use CPU arrays or numerical block stores;
see :doc:`chunking`.
There are two modes:

- **One process, several GPUs:** pass ``devices`` to ``Projector``.
- **One process per GPU:** launch with ``torchrun`` and pass ``distributed=True``.
  Default ``partition="views"`` holds a local view shard and a replicated volume.
  ``partition="space"`` holds an owned volume slab and replicated ray batches.
  ``partition="hybrid"`` combines view and space ownership on the PyTorch backend.
  GPU working buffers hold bounded tiles, pixels and view batches.

The CUDA examples below use NCCL. Portable library calls can use CPU with Gloo
or supported PyTorch accelerators; see :doc:`execution`.

Spatial ownership
-----------------

``partition="space"`` splits the first tensor axis into balanced half-open
slabs. ``global_volume_shape`` is the constructor shape; ``volume_slice`` gives
the rank's global ownership. ``volume_shape`` and ``local_volume_shape`` give
the local input/output shape. Global physical coordinates are preserved, and
empty slabs are valid. Projection SUMs slab contributions; backprojection
generates only the owned volume. Collective order does not depend on local
tile counts. Shared parameter gradients are reduced after the local pixel VJP.

In space mode, divide identical replicated projection losses by ``world_size``
before backward; sum losses on local backprojection slabs. CGLS SUMs owned
volume norms but counts replicated ray norms once. The library's
``total_variation`` exchanges neighboring faces; custom boundary-coupled
regularizers need their own halo exchange. The existing
TV/SIRT examples continue to use views mode. The complete ownership and loss
table is in ``docs/DISTRIBUTED.md`` in the repository.

For replicated CPU state use one process per node controlling multiple GPUs
to avoid extra host copies. For true spatial ownership use one process per
GPU, with streaming inside an owned slab when needed. Disk-backed CGLS and
rank-local checkpoints are available in ``examples/disk_reconstruction.py``;
see :doc:`chunking`. The solver capacity includes ``x,s,p`` together plus ray,
geometry and pipeline state. Hybrid mode uses
``partition_grid=(view_ranks, space_ranks)`` and parent-group rank
``view_rank * space_ranks + space_rank``. Divide projection losses by
``space_world_size`` and backprojection losses by ``view_world_size``.
The disk example accepts only views and space partitions.

One process, several GPUs
-------------------------

No process group is needed. Views are split into contiguous shards in the order
of the device list. Outputs return to the input tensor's device.

.. code-block:: python

   import torch
   from diffct import Projector, spiral_trajectory_3d

   trajectory = spiral_trajectory_3d(
       31, sid=80.0, sdd=128.0, z_range=8.0, n_turns=1.0, device="cpu")
   operator = Projector(trajectory, (32, 32, 32), (48, 40),
                        devices=["cuda:0", "cuda:1"])
   volume = torch.ones((32, 32, 32), device="cuda:0", requires_grad=True)
   sinogram = operator.project(volume)                # (31, 48, 40), cuda:0
   adjoint = operator.backproject(sinogram.detach())  # (32, 32, 32), cuda:0
   sinogram.square().sum().backward()                 # volume.grad on cuda:0

Use the ordinary full-sinogram loss in this mode. No world-size scaling applies.

One process per GPU
-------------------

Save the following as ``distributed_demo.py``. Each rank creates the same
volume, computes its local measurements, and differentiates one global data term.

.. code-block:: python

   import os

   import torch
   import torch.distributed as dist
   from diffct import Projector, spiral_trajectory_3d

   # Select this process's GPU before creating the NCCL process group.
   device = torch.device("cuda", int(os.environ["LOCAL_RANK"]))
   torch.cuda.set_device(device)
   dist.init_process_group("nccl")
   try:
       n_views = 31
       detector_shape = (48, 40)  # (U, V)
       trajectory = spiral_trajectory_3d(
           n_views, sid=80.0, sdd=128.0, z_range=8.0, n_turns=1.0, device="cpu")
       operator = Projector(trajectory, (32, 32, 32), detector_shape,
                            distributed=True)

       truth = torch.ones(operator.volume_shape, device=device)
       local_measurements = operator.project(truth)  # this rank's views only
       volume = torch.zeros_like(truth, requires_grad=True)
       local_prediction = operator.project(volume)

       # Divide the local sum by the global ray count, not by the shard size.
       total_rays = n_views * detector_shape[0] * detector_shape[1]
       local_loss = (0.5 * (local_prediction - local_measurements).square().sum()
                     / total_rays)
       local_loss.backward()  # every rank receives the full volume gradient

       # Reduce a detached copy for reporting only.
       report = local_loss.detach().clone()
       dist.all_reduce(report, op=dist.ReduceOp.SUM)
       if operator.rank == 0:
           print(f"global half-MSE: {report.item():.6g}")
   finally:
       dist.destroy_process_group()

Launch on one node with two GPUs:

.. code-block:: bash

   torchrun --standalone --nproc-per-node=2 distributed_demo.py

Launch on two nodes with two GPUs each. Run the same command on both nodes,
changing only ``--node-rank``. Use the hostname of node 0 as ``--master-addr``:

.. code-block:: bash

   # On node 0:
   torchrun --nnodes=2 --node-rank=0 --nproc-per-node=2 \
       --master-addr=<node0-hostname> --master-port=29500 distributed_demo.py

   # On node 1:
   torchrun --nnodes=2 --node-rank=1 --nproc-per-node=2 \
       --master-addr=<node0-hostname> --master-port=29500 distributed_demo.py

Both nodes need the same checkout and Python environment. Every rank must
join the same collective calls and backward passes, including ranks with zero
views.

Views mode: data on each rank
-----------------------------

- ``operator.view_slice``: the contiguous range of global views on this rank.
  To use full measurements, select ``full_measurements[operator.view_slice]``.
- ``operator.projection_shape``: the shape of ``project()`` output on this rank,
  and the required input shape of ``backproject()``.
- The full volume. All ranks must start with the same volume and trajectory.
  The operator does not broadcast them.
- ``backproject()`` sums the contributions of all ranks and returns the full
  volume on each rank. It is the matched adjoint, not an inverse.

Views mode: loss scaling
------------------------

The projector sums the volume and learnable-geometry gradients across ranks. Scale
the loss as follows:

- **Projection data terms:** form the loss from local views and call ``backward()``
  on every rank. Use a local ``.sum()`` for global least squares. Divide that sum by
  the global number of rays for a mean. A local ``.mean()`` weights unequal shards
  incorrectly and is undefined for an empty shard.
- **Losses on replicated backprojection output:** divide each rank's copy of the
  same output loss by ``operator.world_size`` before ``backward()``.
- **DDP:** do not add another DDP reduction for image or geometry gradients that
  ``Projector`` already reduces.

For regularizers, logging and the full rule set, see
`docs/DISTRIBUTED.md <https://github.com/sypsyp97/diffct/blob/main/docs/DISTRIBUTED.md>`_.

Iterative reconstruction with several GPUs
------------------------------------------

The iterative example runs on one process with several GPUs or with torchrun:

.. code-block:: bash

   python examples/iterative_reconstruction.py --devices 0 1 2 3
   torchrun --nproc-per-node=4 examples/iterative_reconstruction.py

For multi-node Slurm jobs, see ``examples/slurm/multi_node.sbatch`` in the
repository. Set the account, partition and GPU request for your cluster.

For the full multi-GPU and multi-node reference, see
`docs/DISTRIBUTED.md <https://github.com/sypsyp97/diffct/blob/main/docs/DISTRIBUTED.md>`_.


Performance
===========

Time per forward projection, adjoint and CGLS iteration for a 512³ volume on 1, 4
and 8 GPUs. Run the same workload with examples/benchmark_scaling.py.

.. image:: https://sypsyp97.github.io/diffct/_images/scaling_light.png
   :class: only-light performance-chart
   :alt: Time per forward, adjoint and CGLS iteration on 1, 4 and 8 GPUs, 512^3 volume

.. image:: https://sypsyp97.github.io/diffct/_images/scaling_dark.png
   :class: only-dark performance-chart
   :alt: Time per forward, adjoint and CGLS iteration on 1, 4 and 8 GPUs, 512^3 volume

.. list-table::
   :header-rows: 1
   :class: performance-table

   * - 512³ volume
     - 1 GPU
     - 4 GPUs, 1 node
     - 8 GPUs, 2 nodes
   * - Forward projection
     - 566.79 ms
     - 150.76 ms (3.8×)
     - 81.58 ms (6.9×)
   * - Adjoint (backprojection)
     - 1060.77 ms
     - 279.37 ms (3.8×)
     - 224.04 ms (4.7×)
   * - CGLS iteration
     - 1753.73 ms
     - 468.09 ms (3.7×)
     - 337.39 ms (5.2×)

diffct 2.4.3, FP32, A100 64 GB GPUs; circular cone beam, 360 views, 1024×768
detector, pitch 1.25. Default memory settings, without an imposed allocator cap.
Projection times are medians of seven samples after three warmups. CGLS uses
three ten-iteration solves after a two-iteration warmup; the median includes
initialization and is divided by ten. Raw samples and configuration:
`scaling.json <https://github.com/sypsyp97/diffct/blob/main/docs/assets/scaling.json>`_.


Differentiable analytical reconstruction
===========================================

``Projector.project`` and ``Projector.backproject`` implement the cell-constant
Siddon forward operator and its matched adjoint. Analytical FBP and FDK use a
different, voxel-driven interpolation operator. Its backward pass scatters the
same interpolation coefficients, including distance weights and normalization.
It does not substitute the Siddon projector for the analytical transpose.

.. figure:: https://sypsyp97.github.io/diffct/_images/learning.png
   :alt: FBP with the Ram-Lak filter and with a learned filter
   :class: example-figure

   Abdominal slice of TCIA case L213 at 256² (see :ref:`data-sources`), parallel
   beam, 360 views, 2% noise, 200 Adam steps on the filter only. These settings
   differ from ``examples/learnable_reconstruction.py`` below.

The reconstruction APIs return signed images in the selected precision on the
input device. Clamping is a caller decision. In the default native fast mode,
low-precision measurements are promoted to float32 before geometry,
weighting and FFT operations; gradients retain the original input/parameter dtype.
Sinogram, filter and explicit weight gradients propagate through the
pipeline, including higher derivatives with respect to those tensors. Acquisition
geometry and isocenter are fixed; trainable geometry is rejected explicitly.
Use the Siddon ``Projector`` for geometry calibration.

Short-scan distance and offset checks account for physical coordinate scale and
the original trajectory precision, including when float32 geometry is promoted
for float64 measurements. Each check uses only the trajectory components involved
in its calculation; detector-v precision cannot relax the detector-u offset check.
Each source/detector coordinate and detector-axis component counts with its own
rounding along the checked direction, and each allows twice the largest two-view
change that this rounding can cause. Axis components have half the smallest
subnormal as a rounding floor. Eight eps of the in-plane ray length covers
arithmetic and orbit rounding. Both checks also allow a relative tolerance of
1e-4, and the absolute tolerance is at least 1e-5. A change smaller than these
bounds cannot be told apart from rounding: low-precision or far-translated
coordinates that have a component along the checked direction widen the check.
Helical slice centers use float64 independently of
the angle dtype; low-precision angles do not round physical voxel coordinates.

``fbp`` and ``fdk`` support flat detectors, portable FP32/FP64 arithmetic,
view/spatial chunks, several local devices, distributed views/space/hybrid
ownership and numerical ``fbp_into``/``fdk_into`` stores; see :doc:`execution`.
The helical approximation and low-level CUDA helpers retain their original
execution contracts. FDK is an approximate inverse even though its implemented
discrete backward pass transposes the interpolation coefficients.

Use the complete APIs with the same trajectory tensors as ``Projector``:

.. code-block:: python

   from diffct import fbp, fdk, helical_fdk

   image = fbp(sino2d, trajectory2d, (H, W), angles,
               beam="fan", detector_spacing=du, isocenter=(0, 0))
   volume = fdk(sino3d, trajectory3d, (D, H, W), angles,
                detector_spacing=(du, dv), isocenter=(0, 0, 0))
   helical = helical_fdk(sino3d, helix, (D, H, W), unwrapped_angles,
                         detector_spacing=(du, dv), isocenter=(0, 0, 0))

Sinograms use ``(views,U)`` or ``(views,U,V)``; images use ``(H,W)`` or
``(D,H,W)``. Distances, pitches and isocenter share one physical unit.
``projection_weights`` replaces the complete pre-filter cosine weight;
``angular_weights`` replaces angular quadrature in ``fbp``/``fdk``.
``filter_response`` replaces the full real rFFT multiplier, already scaled
in physical units, with length ``pad_factor*U//2+1``. It overrides ``window``.
For circular fan/cone short scans, ``short_scan=True`` adds Parker weighting
before filtering and checks the circular orbit and required angular span.
Native curved data needs covered flat-detector resampling before these APIs.

Building and learning a pipeline
-------------------------------------

Run the examples from a source checkout:

.. code-block:: bash

   python examples/analytical_reconstruction.py --size 64
   python examples/learnable_reconstruction.py --beam parallel --size 32
   python examples/learnable_reconstruction.py --beam fan --size 32
   python examples/learnable_reconstruction.py --beam cone --size 32
   python examples/helical_fdk.py --size 32

``learnable_reconstruction.py`` shows the full pipeline in ordinary PyTorch:
detector cosine weights, a trainable frequency response, FFT filtering, angular
quadrature and the analytical gather. It optimizes the filter and projection
weights for FBP or FDK on a synthetic noisy phantom. This demonstrates
autograd and composition; the fitted parameters are specific to the example.

Helical approximation
---------------------

``helical_fdk`` filters detector rows once, then reconstructs each axial slice
from a full turn centered at the angle where the source crosses that slice.
``helical_view_weights`` exposes those per-slice quadrature weights.
``helical_weighted_backproject`` accepts replacement weights after filtering,
so users can learn or replace this stage independently.

This is a native flat-detector, slice-windowed FDK approximation for an axial
helix with radial detector normals. It is not an implementation of exact
Katsevich inversion, cone-parallel rebinning or conjugate-ray 3D weighting.
The source-z samples and explicit unwrapped angles determine the window;
nominal pitch parameters are not substituted for the actual acquisition.
Incomplete angular windows are rejected. Detector truncation remains possible
and the example uses a detector tall enough to cover its selected windows.

The full-turn window and half redundancy factor are motivated by equation 14
of `Kong et al. (2018)
<https://pmc.ncbi.nlm.nih.gov/articles/PMC6269102/>`_. That paper derives the
formula after cone-parallel rebinning; applying the window directly to the
native cone gather here is an approximation. Increasing pitch or cone angle
can increase reconstruction error.

The implementation processes one slice at a time, keeping the window stage
readable and avoiding a views-by-volume weight tensor. This trades kernel-launch
overhead for bounded temporary storage; inputs, outputs and autograd state
still reside on the CUDA device.

Example results
---------------

The 32-voxel learning examples use 180 views and 20 Adam steps on a synthetic
phantom with Gaussian projection noise. These training results are specific to
that example:

.. list-table::
   :header-rows: 1

   * - Beam
     - Initial PSNR
     - Learned filter/weight PSNR
   * - Parallel FBP
     - 31.96 dB
     - 33.81 dB
   * - Fan FBP
     - 29.90 dB
     - 32.86 dB
   * - Cone FDK
     - 30.24 dB
     - 32.09 dB

The 32-cubed, 720-view, four-turn helical example gives 32.99 dB. The integration
weights of each slice sum to pi. Manual composition and the complete API agree
within 4e-7. Float32 atomic accumulation order depends on the GPU, so the last
digit of these values can change.


```{include} ../MIGRATION.md
```


Examples
========

Each page has the commands, an explanation and one result figure. Figure
captions state the data and settings. The scripts are in the ``examples/``
directory of a source checkout; see the "From source" tab in :doc:`installation`.

Analytical reconstruction
-------------------------

- :doc:`Parallel-beam FBP <fbp_parallel_example>`: FBP with a choice of ramp-filter window.
- :doc:`Fan-beam FBP <fbp_fan_example>`: fan-beam FBP with cosine weights over a full turn.
- :doc:`Cone-beam FDK <fdk_cone_example>`: FDK on a circular orbit, with its cone-beam artifacts.
- :doc:`Learnable FBP <analytical_learning>`: learn the FBP filter through autograd.

Iterative reconstruction
------------------------

- :doc:`Gradient-based reconstruction <gradient_reconstruction>`: the volume is a parameter; project it, compare with the measurement, optimize.
- :doc:`CGLS, SIRT and TV <iterative_reco_cone_example>`: solvers on circular, helical, saddle or sinusoidal 3D scans.
- :doc:`Measured walnut <walnut_example>`: FDK, SIRT, CGLS and TV on 240 measured cone-beam views.

Clinical data and detectors
---------------------------

- :doc:`Clinical helical CT <clinical_helical_example>`: measured patient projections, flying focal spot, curved detector, 8 GPUs.
- :doc:`Curved detector <curved_detector_example>`: projection, matched adjoint and gradients on a native cylindrical detector.

All scripts
-----------

The scripts in the repository ``examples/`` directory show how to use
``diffct``. Run them from the repository root.

- ``quickstart.py``: ``Projector`` basics for parallel, fan and cone beams, with
  the adjoint check and gradients.
- ``curved_detector.py``: native cylindrical cone detector projection, matched
  backprojection, adjoint check and volume/surface gradients on a circular scan.
- ``chunked_reconstruction.py``: CPU-backed CGLS with automatic CUDA tiles,
  optional manual limits and measured peak allocation; see :doc:`chunking`.
- ``disk_reconstruction.py``: blockwise disk CGLS, optional native curves,
  rank-owned spatial slabs and completed-iteration checkpoints.
- ``benchmark_execution.py``: schedule and slab/block comparisons with actual
  copy payloads, pilot metadata, timing and CUDA tensor peaks.
- ``benchmark_trajectory_admission.py``: per-call trajectory validation without
  numerical projection. Use ``--backend torch --frame-device cpu`` on CPU.
- ``benchmark_surface_geometry.py``: forward, adjoint and geometry-gradient
  timings of both curved-detector interfaces on a synthetic helical scan.
- ``benchmark_gpu_overheads.py``: synthetic launch and transfer overheads of
  trajectories, detector surfaces, solvers and tiled FDK.
- ``analytical_reconstruction.py``: parallel-beam FBP, fan-beam FBP and cone-beam
  FDK, with a choice of ramp-filter window.
- ``learnable_reconstruction.py``: build and train a filter and detector weights
  for parallel/fan FBP or cone FDK using public helpers.
- ``helical_fdk.py``: approximate axial helical FDK with explicit slice windows
  and weight gradients.
- ``iterative_reconstruction.py``: 3D cone-beam CGLS, SIRT and TV reconstruction
  on circular, helical, saddle or sinusoidal trajectories.
- ``walnut_reconstruction.py``: FDK, SIRT, CGLS and TV reconstruction of a measured
  walnut scan.
- ``geometry_calibration.py``: recovers per-view angle errors and a detector shift
  from projections, using geometry gradients.
- ``benchmark_projector.py``: compares one GPU with several GPUs.
- ``benchmark_scaling.py``: measures the 512³ README workload on 1, 4 and 8 GPUs,
  with raw timings and numerical checks.
- ``plot_trajectory.py``: plots the trajectory generators (CPU only).

See ``examples/README.md`` in the repository for all launch modes.

.. _data-sources:

Data sources
------------

- **Patient data:** cases L213 (abdomen) and C002 (chest) of McCollough et al.,
  *Low Dose CT Image and Projection Data (LDCT-and-Projection-data)*, Version 7,
  The Cancer Imaging Archive, 2020, https://doi.org/10.7937/9npb-2637; dataset
  paper: Moen et al., Medical Physics, 2020, https://doi.org/10.1002/mp.14594.
  Chest and Liver subsets: `CC BY 4.0 <https://creativecommons.org/licenses/by/4.0/>`_.
  The data were acquired under grants EB017095 and EB017185 (Cynthia McCollough, PI)
  from the National Institute of Biomedical Imaging and Bioengineering. "Ground
  truth" panels show the vendor images, converted from HU to water units; all
  other panels are reconstructions or simulations computed with diffct.
- **Walnut:** Meaney 2022, Zenodo 6986012, CC BY 4.0; see ``examples/data/NOTICE``.
- **Shepp-Logan phantom:** computed by ``examples/_common.py``.

The figures are made by ``docs/assets/gallery/compute.py`` and
``docs/assets/gallery/render.py``; see ``docs/assets/gallery/NOTICE``.

.. toctree::
   :hidden:
   :caption: Analytical reconstruction

   fbp_parallel_example
   fbp_fan_example
   fdk_cone_example

.. toctree::
   :hidden:
   :caption: Iterative reconstruction

   gradient_reconstruction
   iterative_reco_cone_example
   walnut_example

.. toctree::
   :hidden:
   :caption: Clinical data and detectors

   clinical_helical_example
   curved_detector_example


Parallel Beam Filtered Backprojection (FBP)
===========================================

This page describes the ``examples/analytical_reconstruction.py`` pipeline. For
filter options and helper signatures, see :doc:`api`.

.. figure:: https://sypsyp97.github.io/diffct/_images/fbp_parallel.png
   :alt: Parallel-beam FBP of the Shepp-Logan phantom and of a patient slice
   :class: example-figure

   Ground truth: Shepp-Logan, and an abdominal slice of TCIA case L213 (see
   :ref:`data-sources`). 512², 720 views over 180°, 1536 detector bins at half the
   pixel size, Shepp-Logan window. PSNR in the lower corner.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`. ``--figure`` needs Matplotlib: ``pip install matplotlib``.

.. code-block:: bash

   python examples/analytical_reconstruction.py --size 64 --window shepp-logan
   python examples/analytical_reconstruction.py --size 128 --figure analytical.png

The script runs parallel, fan and cone reconstruction on one CUDA GPU and prints
PSNR against the synthetic phantom. ``--figure`` requires Matplotlib. ``--window``
accepts ``ram-lak``, ``shepp-logan``, ``cosine``, ``hamming`` or ``hann``.

Pipeline
--------

The example uses a 180-degree parallel-beam scan with 360 views. The detector has
``3 * size`` cells and a pitch of 0.5 voxels.

1. Project a phantom slice with ``Projector(beam="parallel")``.
2. Filter along the detector dimension (``dim=1``) with
   ``ramp_filter_1d(..., sample_spacing=pitch, pad_factor=2)``.
3. Multiply by ``angular_integration_weights(half_turn,
   redundant_full_scan=False)``.
4. Call ``parallel_weighted_backproject`` on the sinogram's CUDA device.

The analytical helpers are not the ``Projector.backproject()`` adjoint. Do not
use the adjoint in this pipeline, and do not add extra scale factors. The helpers
apply their own normalization.

Source
------

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/analytical_reconstruction.py
   :language: python
   :linenos:
   :caption: Parallel FBP, fan FBP and cone FDK


Fan Beam Filtered Backprojection (FBP)
======================================

This page describes the fan-beam part of ``examples/analytical_reconstruction.py``.
For filter options and helper signatures, see :doc:`api`.

.. figure:: https://sypsyp97.github.io/diffct/_images/fbp_fan.png
   :alt: Fan-beam FBP of the Shepp-Logan phantom and of a patient slice
   :class: example-figure

   Ground truth: Shepp-Logan, and a chest slice of TCIA case C002 (see
   :ref:`data-sources`), lung window. 512², 720 views over 360°, magnification 1.6,
   Shepp-Logan window. PSNR in the lower corner.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`. ``--figure`` needs Matplotlib: ``pip install matplotlib``.

.. code-block:: bash

   python examples/analytical_reconstruction.py --size 64 --window shepp-logan
   python examples/analytical_reconstruction.py --size 128 --figure analytical.png

The script runs parallel, fan and cone reconstruction on one CUDA GPU and prints
PSNR against the synthetic phantom. ``--figure`` requires Matplotlib.

Pipeline
--------

The example uses a full-turn circular fan-beam scan with 360 views,
``sid = 2.5 * size``, ``sdd = 4 * size``, ``3 * size`` detector cells and a
pitch of 0.8. Distances use the same units as the unit-spaced voxels.

1. Project a phantom slice with ``Projector(beam="fan")``.
2. Apply ``fan_cosine_weights`` for the detector pitch and source-detector distance.
3. Apply ``ramp_filter_1d`` along ``dim=1``, with ``sample_spacing=pitch``
   and ``pad_factor=2``.
4. Multiply by ``angular_integration_weights`` with ``redundant_full_scan=True``.
5. Call ``fan_weighted_backproject`` for the distance-weighted gather.

This example uses a full scan. Short-scan reconstruction needs Parker weights and
suitable angular coverage. Changing the view count does not enable it.

The analytical helpers are not the ``Projector.backproject()`` adjoint. Do not
use the adjoint in this pipeline, and do not add extra scale factors.

Source
------

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/analytical_reconstruction.py
   :language: python
   :linenos:
   :caption: Parallel FBP, fan FBP and cone FDK


Cone Beam Feldkamp-Davis-Kress (FDK)
====================================

This page describes the cone-beam part of ``examples/analytical_reconstruction.py``.
For filter options and helper signatures, see :doc:`api`.

.. figure:: https://sypsyp97.github.io/diffct/_images/fdk_cone.png
   :alt: Cone-beam projection and FDK slices of a 3D phantom
   :class: example-figure

   3D Shepp-Logan phantom, 256³, 720 views on a circular orbit. The off-centre
   slice and the coronal section show the approximation error of FDK away from the
   orbit plane.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`. ``--figure`` needs Matplotlib: ``pip install matplotlib``.

.. code-block:: bash

   python examples/analytical_reconstruction.py --size 64 --window shepp-logan
   python examples/analytical_reconstruction.py --size 128 --figure analytical.png

The script runs parallel, fan and cone reconstruction on one CUDA GPU and prints
PSNR against the synthetic phantom. ``--figure`` requires Matplotlib.

Pipeline
--------

The example uses a flat detector and a circular cone-beam scan with 360 views.
It does not reconstruct native curved-detector data directly. To reuse this
pipeline for a curved detector, first resample line integrals along matching
rays onto a covered virtual flat grid. This adds interpolation error; FDK still
retains its circular cone-beam approximation. For native curved projection and
matched backprojection, see :ref:`detector-surfaces` and :doc:`examples`.

The example projects a 3D phantom with ``Projector`` into a ``(views, U, V)``
sinogram. The helper ``examples/_common.py::fdk`` then calls :func:`diffct.reconstruction.fdk`,
which performs these steps:

1. Apply detector cosine weights computed from the actual per-view rays.
2. Filter along the detector u axis (``dim=1``) with ``sample_spacing=du``
   and ``pad_factor=2``. The last axis is v, not u.
3. Apply ``angular_integration_weights``. A full turn includes the ``1/2``
   redundancy factor.
4. Call ``cone_weighted_backproject``. The example clamps negative values to zero.

The same stages are public helpers. See :doc:`analytical_learning` to replace
or learn the weights and filter.

FDK is approximate for a circular cone-beam scan away from the central plane.
Large cone angles increase artifacts. On non-circular trajectories, FDK is a
heuristic baseline, not an exact reconstruction.

The analytical helpers are not the ``Projector.backproject()`` adjoint. Do not
use the adjoint in this pipeline, and do not add extra scale factors.

Source
------

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/analytical_reconstruction.py
   :language: python
   :linenos:
   :caption: Parallel FBP, fan FBP and cone FDK


Gradient-Based Iterative Reconstruction
=======================================

The volume is a tensor of parameters. The differentiable forward projector maps
it to projections, a loss compares them with the measurement, and a PyTorch
optimizer updates the volume. Autograd computes the gradient with the matched
adjoint ``A^T``, so no backprojection code is written by hand.

The same loop takes other losses, regularizers, constraints and network
parameters. For least squares without extra terms, the :func:`~diffct.solvers.cgls`
and :func:`~diffct.solvers.sirt` solvers converge in fewer iterations.

Objective
---------

For a projector ``A`` and measured projections ``y``, the examples minimize the
mean squared error over a nonnegative volume ``x``:

.. math::
   \hat{x} = \arg\min_{x \ge 0} \tfrac{1}{M} \lVert A x - y \rVert_2^2

``M`` is the number of measured values. After each Adam step, ``x`` is clamped
to zero or more.

2D parallel and fan beam
------------------------

.. figure:: https://sypsyp97.github.io/diffct/_images/gradient_2d.png
   :alt: Ground truth, Adam steps and data loss for parallel and fan beams
   :class: example-figure

   The loop below at 512², with 720 views, Adam learning rate 0.1, 1000 steps and
   the abdominal slice of TCIA case L213 as ground truth (see :ref:`data-sources`). PSNR uses the
   peak of the ground truth. The code below uses a 128² square so that it runs in
   seconds.

.. tab-set::

   .. tab-item:: Parallel beam

      .. code-block:: python

         import torch
         from diffct import Projector, circular_trajectory_2d_parallel

         n, n_views, n_det = 128, 360, 384
         trajectory = circular_trajectory_2d_parallel(n_views, device="cuda")
         A = Projector(trajectory, (n, n), n_det, beam="parallel", detector_spacing=0.5)

         truth = torch.zeros((n, n), device="cuda")
         truth[40:88, 40:88] = 1.0
         measured = A.project(truth).detach()

         image = torch.zeros((n, n), device="cuda", requires_grad=True)
         optimizer = torch.optim.Adam([image], lr=0.1)
         for step in range(500):
             optimizer.zero_grad()
             loss = 0.5 * (A.project(image) - measured).square().sum()
             loss.backward()
             optimizer.step()
             with torch.no_grad():
                 image.clamp_(min=0)

   .. tab-item:: Fan beam

      .. code-block:: python

         import torch
         from diffct import Projector, circular_trajectory_2d_fan

         n, n_views = 128, 360
         sid, sdd, n_det, pitch = 2.5 * n, 4.0 * n, 3 * n, 0.8
         trajectory = circular_trajectory_2d_fan(n_views, sid, sdd, device="cuda")
         A = Projector(trajectory, (n, n), n_det, beam="fan", detector_spacing=pitch)

         truth = torch.zeros((n, n), device="cuda")
         truth[40:88, 40:88] = 1.0
         measured = A.project(truth).detach()

         image = torch.zeros((n, n), device="cuda", requires_grad=True)
         optimizer = torch.optim.Adam([image], lr=0.1)
         for step in range(500):
             optimizer.zero_grad()
             loss = 0.5 * (A.project(image) - measured).square().sum()
             loss.backward()
             optimizer.step()
             with torch.no_grad():
                 image.clamp_(min=0)

The 2D loops use the summed squared error; its gradient differs from the mean
only by a constant factor, and Adam steps are nearly invariant to such a factor.

3D measured walnut
------------------

.. figure:: https://sypsyp97.github.io/diffct/_images/gradient_walnut.png
   :alt: FDK reference, Adam steps and MSE for the measured walnut
   :class: example-figure

   240 measured cone-beam views of a walnut, 256³ volume, Adam with learning
   rate 0.01. The FDK column is a reference only; the loop starts from zero.

The measured projections are scaled so that the walnut shell is about 1, as in
``examples/walnut_reconstruction.py``:

.. code-block:: python

   import torch
   from diffct import Projector

   # trajectory, detector shape and spacing from examples/walnut_reconstruction.py
   A = Projector(trajectory, (256, 256, 256), detector, detector_spacing=(du, dv), voxel_spacing=voxel)

   volume = torch.zeros(A.volume_shape, device="cuda", requires_grad=True)
   optimizer = torch.optim.Adam([volume], lr=0.01)
   for step in range(500):
       optimizer.zero_grad()
       loss = torch.nn.functional.mse_loss(A.project(volume), measured)
       loss.backward()
       optimizer.step()
       with torch.no_grad():
           volume.clamp_(min=0)

The full script is the ``gradient`` command of ``docs/assets/gallery/compute.py``.

Variations
----------

- **Regularization:** add a penalty to the loss, for example
  ``loss = loss + weight * diffct.total_variation(volume, voxel_spacing=voxel)``.
  :func:`~diffct.solvers.tv_reconstruction` packages this loop with residual
  history and checkpoints.
- **Learned components:** replace the volume tensor with the output of a network,
  or make trajectory tensors trainable to calibrate the geometry; see
  :doc:`trajectories`.
- **Large volumes and several GPUs:** the same loop works with the
  ``devices`` and ``distributed`` options of ``Projector``; see :doc:`multi_gpu`.
- **Stopping:** the learning rate and step count are starting values. On real data,
  stop with a rule that does not use a ground truth.


Cone Beam Iterative Reconstruction
==================================

``examples/iterative_reconstruction.py`` reconstructs a 3D cone-beam phantom with
``diffct.Projector``. It runs on any supported trajectory, on one GPU, on several
GPUs in one process, or on several processes and nodes.

.. figure:: https://sypsyp97.github.io/diffct/_images/iterative_cone.png
   :alt: Axial and coronal slices of a simulated helical scan
   :class: example-figure

   Simulated helical walnut scan from :doc:`walnut_example`: 720 views, 256³,
   1% noise.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`. ``--figure`` needs Matplotlib: ``pip install matplotlib``.

.. code-block:: bash

   python examples/iterative_reconstruction.py --trajectory helical
   python examples/iterative_reconstruction.py --trajectory circular --devices 0 1 2 3
   torchrun --nproc-per-node=4 examples/iterative_reconstruction.py

``--devices`` and ``torchrun`` cannot be used together. For multiple nodes, see
:doc:`multi_gpu`.

Options
-------

.. list-table::
   :header-rows: 1

   * - Option
     - Meaning
   * - ``--size``
     - Volume edge length in voxels (default 128).
   * - ``--views``
     - Number of views (default 360).
   * - ``--trajectory``
     - ``circular``, ``helical``, ``saddle`` or ``sinusoidal`` (default ``helical``).
   * - ``--algorithms``
     - One or more of ``cgls``, ``sirt`` and ``tv`` (default: all).
   * - ``--iterations``
     - Override the per-algorithm defaults (CGLS 30, SIRT 200, TV 200).
   * - ``--tv-weight``
     - TV penalty weight for ``tv`` (default 1.0).
   * - ``--noise``
     - Gaussian noise, relative to the largest projection (default 0.0).
   * - ``--phantom``
     - ``(size, size, size)`` ``.npy`` volume with values in [0, 1], used instead
       of the Shepp-Logan phantom.
   * - ``--figure``
     - Save central slices to this PNG. Requires Matplotlib.

Pipeline
--------

1. Project the phantom with ``Projector`` (cone beam, on the selected trajectory).
2. Reconstruct with CGLS and SIRT. Both use the projector and its adjoint.
3. Reconstruct with TV-regularized nonnegative least squares. The script uses
   Adam and autograd for this solver.
4. When the full sinogram is on one GPU, compute an FDK baseline. On non-circular
   trajectories, FDK is a heuristic baseline, not an exact reconstruction.

Notes
-----

- The iteration counts are demonstration settings, not a general convergence
  schedule. Choose a stopping rule for your data.
- The TV weight is an example setting. Choose it for your geometry, resolution,
  voxel spacing and data scale. The solver uses half the mean squared residual
  per ray and a mean TV term. Its voxel differences are not divided by voxel spacing.
- The script keeps complete input/output tensors on the process's CUDA device.
  Native working buffers may stream; CGLS state can fall back to CPU.
- To adapt the loss for several ranks, see :doc:`multi_gpu`.

Source
------

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/iterative_reconstruction.py
   :language: python
   :linenos:
   :caption: Iterative Reconstruction Example (3D cone beam, any trajectory)


Measured Walnut
===============

``examples/walnut_reconstruction.py`` reconstructs a measured walnut from 240
cone-beam projections on a full circular orbit. The data are from Meaney 2022
(Zenodo 6986012, CC BY 4.0; see ``examples/data/NOTICE``). The scan geometry is
passed to ``Projector`` as explicit per-view tensors, the same way as a calibrated
trajectory.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`. ``--figure`` needs Matplotlib: ``pip install matplotlib``.

.. code-block:: bash

   python examples/walnut_reconstruction.py --figure walnut.png
   python examples/walnut_reconstruction.py --size 512 --devices 0 1 2 3

The figure shows axial and coronal centre slices. The reconstructions are:

- FDK with a Hann window.
- SIRT (200 iterations).
- CGLS (20 iterations).
- TV-regularized least squares (300 Adam steps, weight 0.3).

Use ``--algorithms`` with no values to run FDK only. Use ``--window`` to change the
FDK filter.

.. image:: https://sypsyp97.github.io/diffct/_images/walnut.png
   :alt: Measured walnut: FDK, SIRT, CGLS and TV reconstructions

Simulated helical scan
----------------------

``--save-volume`` stores the FDK volume scaled to [0, 1]. Use it as the phantom for
a simulated helical scan:

.. code-block:: bash

   python examples/walnut_reconstruction.py --algorithms --save-volume walnut256.npy
   python examples/iterative_reconstruction.py --size 256 --views 720 --trajectory helical \
       --noise 0.01 --phantom walnut256.npy --figure helical.png

.. image:: _static/gallery/iterative_cone.png
   :alt: Simulated helical scan of the walnut: FDK, CGLS, SIRT and TV reconstructions


Clinical Helical CT
===================

This example reconstructs measured projections of patients from a clinical helical
CT scanner. The same script also reconstructs case C002 (chest, Siemens, flying
focal spot, 11072 views). The data are the full-dose projections of
LDCT-and-Projection-data (see :ref:`data-sources`). Each view is one detector
readout in the DICOM-CT-PD format.

.. figure:: https://sypsyp97.github.io/diffct/_images/helical_L213.png
   :alt: Axial, coronal and sagittal slices of a measured abdominal helical scan
   :class: example-figure

   Case L213 (abdomen, GE): 7140 views, 888 × 64 cylindrical detector, 1 mm voxels,
   30 CGLS iterations without regularization. Soft-tissue window.

Geometry
--------

The scanner geometry goes to ``Projector`` as per-view tensors, without
rebinning to a circular orbit or a flat detector:

- **Source:** the focal spot of each view, including the flying focal spot
  offsets in angle, radius and z.
- **Detector:** a cylindrical surface centred on the nominal focal spot, given
  as a ``detector_surface`` callback. Projection and backprojection use these
  pixel positions directly.
- **Views:** the views are split over the GPUs with ``distributed=True``;
  :func:`~diffct.solvers.cgls` runs on the distributed operator.

The detector column of the central ray is 368.56 for C002 and 442.2 for L213, as
0-based stored columns. We estimated both values beforehand with gradient-based
calibration of the column offset; that calibration is not part of this script.

.. literalinclude:: https://sypsyp97.github.io/diffct/../assets/gallery/compute.py
   :language: python
   :pyobject: helical
   :caption: docs/assets/gallery/compute.py

Run it
------

1. Download the full-dose projection and image series of a Chest or Liver case
   from The Cancer Imaging Archive.
2. Convert one case to arrays:

   .. code-block:: bash

      pip install pydicom
      python docs/assets/gallery/ldct_prep.py <projection dir> <image dir> ldct/L213

3. Reconstruct on 8 GPUs:

   .. code-block:: bash

      torchrun --nproc-per-node=8 docs/assets/gallery/compute.py helical \
          --case ldct/L213 --center 442.2 --axial-z -418 -364 -307 --out arrays


Native Curved Detector
======================

``examples/curved_detector.py`` projects a phantom onto a native cylindrical cone
detector, checks the matched backprojection and computes volume and surface
gradients on a circular scan.

.. figure:: https://sypsyp97.github.io/diffct/_images/curved_detector.png
   :alt: Projections of a patient volume on a cylindrical detector
   :class: example-figure

   Vendor images of TCIA case L213 (see :ref:`data-sources`) as a volume, projected
   onto a cylindrical detector centred on the source: 888 × 625 pixels, radius
   949 mm, source 541 mm from the isocentre. Made by the ``volumes`` command of
   ``docs/assets/gallery/compute.py``; the script below uses a smaller geometry.

Run it
------

Run the commands from the root of a source checkout; ``pip install diffct`` does not
install the ``examples/`` directory. See the "From source" tab in
:doc:`installation`.

.. code-block:: bash

   python examples/curved_detector.py

This small CUDA example uses the native curved operators directly, without
resampling or FDK. See :ref:`detector-surfaces` for the callback contract.

Source
------

.. literalinclude:: https://sypsyp97.github.io/diffct/../../examples/curved_detector.py
   :language: python
   :linenos:
   :caption: Native cylindrical detector operators


API Reference
=============

.. raw:: html

   <script>
   // The reference used to be one page. Old links such as api.html#diffct.Projector.project
   // open the generated page of the object, at the same anchor.
   (function () {
     var pages = ["diffct.Projector", "diffct.TensorStore", "diffct.NpyStore", "diffct.ParameterizedSurface", "diffct.ParallelProjectorFunction", "diffct.ParallelBackprojectorFunction", "diffct.FanProjectorFunction", "diffct.FanBackprojectorFunction", "diffct.ConeProjectorFunction", "diffct.ConeBackprojectorFunction", "diffct.geometry.circular_trajectory_3d", "diffct.geometry.circular_trajectory_3d_parallel", "diffct.geometry.random_trajectory_3d", "diffct.geometry.spiral_trajectory_3d", "diffct.geometry.sinusoidal_trajectory_3d", "diffct.geometry.saddle_trajectory_3d", "diffct.geometry.custom_trajectory_3d", "diffct.geometry.circular_trajectory_2d_parallel", "diffct.geometry.sinusoidal_trajectory_2d_parallel", "diffct.geometry.custom_trajectory_2d_parallel", "diffct.geometry.circular_trajectory_2d_fan", "diffct.geometry.sinusoidal_trajectory_2d_fan", "diffct.geometry.custom_trajectory_2d_fan", "diffct.geometry.transform_trajectory", "diffct.analytical.detector_coordinates_1d", "diffct.analytical.angular_integration_weights", "diffct.analytical.fan_cosine_weights", "diffct.analytical.cone_cosine_weights", "diffct.analytical.parker_weights", "diffct.analytical.ramp_filter_1d", "diffct.analytical.parallel_weighted_backproject", "diffct.analytical.fan_weighted_backproject", "diffct.analytical.cone_weighted_backproject", "diffct.reconstruction.fbp", "diffct.reconstruction.fdk", "diffct.reconstruction.fbp_into", "diffct.reconstruction.fdk_into", "diffct.reconstruction.helical_fdk", "diffct.reconstruction.helical_view_weights", "diffct.reconstruction.helical_weighted_backproject", "diffct.solvers.SolverResult", "diffct.solvers.cgls", "diffct.solvers.sirt", "diffct.solvers.tv_reconstruction", "diffct.solvers.cgls_into", "diffct.solvers.total_variation"];
     var target = decodeURIComponent(window.location.hash.slice(1));
     for (var name = target; name.indexOf(".") > 0; name = name.slice(0, name.lastIndexOf("."))) {
       if (pages.indexOf(name) >= 0) {
         window.location.replace("generated/" + name + ".html#" + target);
         return;
       }
     }
   })();
   </script>

High-level Projector
--------------------

.. currentmodule:: diffct

.. autosummary::
   :toctree: generated
   :template: api_class.rst
   :nosignatures:

   Projector

``Projector`` takes trajectory tensors as described in :doc:`trajectories`. For
multiple GPUs and nodes, see :doc:`multi_gpu`.
Each operation also accepts ``trajectory=...`` to replace the full trajectory
for that call, preserving its view count and the constructor default.

``detector_surface`` defines native pixel positions for both projection and
matched backprojection. See :ref:`detector-surfaces` for its coordinate and
shape contract, and :doc:`examples` for a runnable cylindrical detector.

Native CPU data uses streamed CUDA execution; ``backend="torch"`` computes
on its PyTorch device. ``precision`` selects fast FP32 or strict FP32/FP64.
Leading batch/channel axes and per-axis voxel spacing are supported; see
:doc:`execution`. ``volume_chunk_shape``,
``view_chunk_size`` and ``detector_chunk_shape`` specify optional limits;
``schedule`` allows operation-specific comparisons. ``partition="space"``
assigns rank-local volume slabs. See :doc:`chunking` for memory, storage,
geometry and autograd requirements, and :doc:`multi_gpu` for ownership.

Complete Reconstruction Pipelines
------------------------------------

.. currentmodule:: diffct.reconstruction

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   fbp
   fdk
   fbp_into
   fdk_into
   helical_fdk
   helical_view_weights
   helical_weighted_backproject

Numerical Solvers
-----------------

.. currentmodule:: diffct.solvers

.. autosummary::
   :toctree: generated
   :template: api_class.rst
   :nosignatures:

   SolverResult

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   cgls
   sirt
   tv_reconstruction
   cgls_into
   total_variation

Geometry Helpers
----------------

.. currentmodule:: diffct.geometry

**3D trajectories**

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   circular_trajectory_3d
   circular_trajectory_3d_parallel
   random_trajectory_3d
   spiral_trajectory_3d
   sinusoidal_trajectory_3d
   saddle_trajectory_3d
   custom_trajectory_3d

**2D trajectories**

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   circular_trajectory_2d_parallel
   sinusoidal_trajectory_2d_parallel
   custom_trajectory_2d_parallel
   circular_trajectory_2d_fan
   sinusoidal_trajectory_2d_fan
   custom_trajectory_2d_fan

**Rigid transforms**

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   transform_trajectory

Block stores and batched surfaces
---------------------------------

.. currentmodule:: diffct

.. autosummary::
   :toctree: generated
   :template: api_class.rst
   :nosignatures:

   TensorStore
   NpyStore
   ParameterizedSurface

``project_into`` and ``backproject_into`` use caller-owned numerical outputs.
``ParameterizedSurface`` exposes mutable geometry as explicit tensor
parameters and samples bounded global view/pixel rectangles. Legacy two-argument
surface callbacks remain available. See :ref:`detector-surfaces`.

Analytical Reconstruction Helpers
---------------------------------

.. currentmodule:: diffct.analytical

These helpers build the pre-weights, angular weights, filter and weighted
backprojection of an analytical FBP or FDK pipeline. The weighted backprojection
helpers support sinogram autograd with dedicated transpose kernels, including
higher derivatives. They use fixed geometry and reject geometry/isocenter
gradients. The returned images are already scaled. See :doc:`analytical_learning`
for complete APIs and learnable filter/weight examples.

Analytical FBP and FDK assume a flat detector and an acquisition model. These
helpers do not accept ``detector_surface``. Native curved data must be resampled
along matching rays onto a covered virtual flat detector before using these
flat-detector helpers; resampling adds interpolation error. The native
``Projector`` forward/adjoint pair needs no such conversion. A non-circular
trajectory does not make the result exact. The helpers accept the same per-view
``(src_pos, det_center, det_u[, det_v])`` arrays as ``Projector``.

Fan and cone backprojection accept a keyword-only ``isocenter`` vector in physical
units. With ``None`` (default), it is estimated from the detector normals. Give it
explicitly for non-circular or ambiguous geometry, including a single view. The
three weighted backprojection helpers need at least two detector bins along each
interpolated axis, including the cone detector's v axis.

.. autosummary::
   :toctree: generated
   :template: api_function.rst
   :nosignatures:

   detector_coordinates_1d
   angular_integration_weights
   fan_cosine_weights
   cone_cosine_weights
   parker_weights
   ramp_filter_1d
   parallel_weighted_backproject
   fan_weighted_backproject
   cone_weighted_backproject

Ramp Filter Options
-------------------

``ramp_filter_1d`` is the filter used by the analytical examples. Its 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``, for example ``du`` for cone beam.
    The output is in physical units. Pass ``1.0`` to work in sample units.

``pad_factor``
    Zero-pads the signal to ``pad_factor * N`` samples along ``dim``. Use ``2`` for
    cone-beam FDK. Use ``4`` when objects lie close to the detector edge.

``use_rfft``
    Uses the real FFT for real inputs. Set ``False`` only to force the complex FFT
    path.

The ``window`` argument selects a window applied to the Ram-Lak ramp:

.. list-table::
   :header-rows: 1

   * - ``window``
     - Effect
   * - ``None`` or ``"ram-lak"``
     - No window. Sharpest, with the most noise.
   * - ``"cosine"``
     - Half-cosine roll-off.
   * - ``"shepp-logan"``
     - Sinc roll-off. Classical choice.
   * - ``"hamming"``
     - Milder than Hann.
   * - ``"hann"`` or ``"hanning"``
     - Strongest high-frequency suppression.

Low-level Functions
-------------------

These ``autograd.Function`` wrappers take explicit trajectory tensors. Use the
``forward`` signature as the argument order for ``.apply(...)``. For new code,
prefer ``Projector``, which checks inputs and manages devices.

.. currentmodule:: diffct

.. autosummary::
   :toctree: generated
   :template: api_autograd_function.rst
   :nosignatures:

   ParallelProjectorFunction
   ParallelBackprojectorFunction
   FanProjectorFunction
   FanBackprojectorFunction
   ConeProjectorFunction
   ConeBackprojectorFunction


Projector
=========

.. currentmodule:: diffct

.. autoclass:: Projector
   :members:
   :special-members: __call__

fbp
===

.. currentmodule:: diffct.reconstruction

.. autofunction:: fbp

fdk
===

.. currentmodule:: diffct.reconstruction

.. autofunction:: fdk

fbp\_into
=========

.. currentmodule:: diffct.reconstruction

.. autofunction:: fbp_into

fdk\_into
=========

.. currentmodule:: diffct.reconstruction

.. autofunction:: fdk_into

helical\_fdk
============

.. currentmodule:: diffct.reconstruction

.. autofunction:: helical_fdk

helical\_view\_weights
======================

.. currentmodule:: diffct.reconstruction

.. autofunction:: helical_view_weights

helical\_weighted\_backproject
==============================

.. currentmodule:: diffct.reconstruction

.. autofunction:: helical_weighted_backproject

SolverResult
============

.. currentmodule:: diffct.solvers

.. autoclass:: SolverResult
   :members:
   :special-members: __call__

cgls
====

.. currentmodule:: diffct.solvers

.. autofunction:: cgls

sirt
====

.. currentmodule:: diffct.solvers

.. autofunction:: sirt

tv\_reconstruction
==================

.. currentmodule:: diffct.solvers

.. autofunction:: tv_reconstruction

cgls\_into
==========

.. currentmodule:: diffct.solvers

.. autofunction:: cgls_into

total\_variation
================

.. currentmodule:: diffct.solvers

.. autofunction:: total_variation

circular\_trajectory\_3d
========================

.. currentmodule:: diffct.geometry

.. autofunction:: circular_trajectory_3d

circular\_trajectory\_3d\_parallel
==================================

.. currentmodule:: diffct.geometry

.. autofunction:: circular_trajectory_3d_parallel

random\_trajectory\_3d
======================

.. currentmodule:: diffct.geometry

.. autofunction:: random_trajectory_3d

spiral\_trajectory\_3d
======================

.. currentmodule:: diffct.geometry

.. autofunction:: spiral_trajectory_3d

sinusoidal\_trajectory\_3d
==========================

.. currentmodule:: diffct.geometry

.. autofunction:: sinusoidal_trajectory_3d

saddle\_trajectory\_3d
======================

.. currentmodule:: diffct.geometry

.. autofunction:: saddle_trajectory_3d

custom\_trajectory\_3d
======================

.. currentmodule:: diffct.geometry

.. autofunction:: custom_trajectory_3d

circular\_trajectory\_2d\_parallel
==================================

.. currentmodule:: diffct.geometry

.. autofunction:: circular_trajectory_2d_parallel

sinusoidal\_trajectory\_2d\_parallel
====================================

.. currentmodule:: diffct.geometry

.. autofunction:: sinusoidal_trajectory_2d_parallel

custom\_trajectory\_2d\_parallel
================================

.. currentmodule:: diffct.geometry

.. autofunction:: custom_trajectory_2d_parallel

circular\_trajectory\_2d\_fan
=============================

.. currentmodule:: diffct.geometry

.. autofunction:: circular_trajectory_2d_fan

sinusoidal\_trajectory\_2d\_fan
===============================

.. currentmodule:: diffct.geometry

.. autofunction:: sinusoidal_trajectory_2d_fan

custom\_trajectory\_2d\_fan
===========================

.. currentmodule:: diffct.geometry

.. autofunction:: custom_trajectory_2d_fan

transform\_trajectory
=====================

.. currentmodule:: diffct.geometry

.. autofunction:: transform_trajectory

TensorStore
===========

.. currentmodule:: diffct

.. autoclass:: TensorStore
   :members:
   :special-members: __call__

NpyStore
========

.. currentmodule:: diffct

.. autoclass:: NpyStore
   :members:
   :special-members: __call__

ParameterizedSurface
====================

.. currentmodule:: diffct

.. autoclass:: ParameterizedSurface
   :members:
   :special-members: __call__

detector\_coordinates\_1d
=========================

.. currentmodule:: diffct.analytical

.. autofunction:: detector_coordinates_1d

angular\_integration\_weights
=============================

.. currentmodule:: diffct.analytical

.. autofunction:: angular_integration_weights

fan\_cosine\_weights
====================

.. currentmodule:: diffct.analytical

.. autofunction:: fan_cosine_weights

cone\_cosine\_weights
=====================

.. currentmodule:: diffct.analytical

.. autofunction:: cone_cosine_weights

parker\_weights
===============

.. currentmodule:: diffct.analytical

.. autofunction:: parker_weights

ramp\_filter\_1d
================

.. currentmodule:: diffct.analytical

.. autofunction:: ramp_filter_1d

parallel\_weighted\_backproject
===============================

.. currentmodule:: diffct.analytical

.. autofunction:: parallel_weighted_backproject

fan\_weighted\_backproject
==========================

.. currentmodule:: diffct.analytical

.. autofunction:: fan_weighted_backproject

cone\_weighted\_backproject
===========================

.. currentmodule:: diffct.analytical

.. autofunction:: cone_weighted_backproject

ParallelProjectorFunction
=========================

.. currentmodule:: diffct

.. autoclass:: ParallelProjectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

ParallelBackprojectorFunction
=============================

.. currentmodule:: diffct

.. autoclass:: ParallelBackprojectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

FanProjectorFunction
====================

.. currentmodule:: diffct

.. autoclass:: FanProjectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

FanBackprojectorFunction
========================

.. currentmodule:: diffct

.. autoclass:: FanBackprojectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

ConeProjectorFunction
=====================

.. currentmodule:: diffct

.. autoclass:: ConeProjectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

ConeBackprojectorFunction
=========================

.. currentmodule:: diffct

.. autoclass:: ConeBackprojectorFunction
   :class-doc-from: init
   :members:
   :undoc-members:
   :show-inheritance:

```{include} ../../CHANGELOG.md
```


:orphan:

Moved
=====

.. raw:: html

   <meta http-equiv="refresh" content="0; url=gradient_reconstruction.html">

This page moved to :doc:`gradient_reconstruction`.


:orphan:

Moved
=====

.. raw:: html

   <meta http-equiv="refresh" content="0; url=gradient_reconstruction.html">

This page moved to :doc:`gradient_reconstruction`.