diffct: Differentiable CT Operators from Circular Orbits to Arbitrary Trajectories

1. Introduction

Computed tomography (CT) reconstruction still relies on three closely related operator types that must agree with each other numerically: a forward model that maps an attenuation field to projection data, a pure adjoint backprojection that propagates information back to image space, and reconstruction operators that approximate inversion under additional geometry assumptions [1]. Modern learning-based pipelines add a differentiability requirement so that image-domain variables, regularization parameters, and surrounding network parameters can be optimized end to end. diffct addresses this intersection by combining CUDA kernels, PyTorch autograd wrappers, and analytical reconstruction helpers in one small library [2,3].

The current repository has two public branch lineages. The main branch is the stable circular-orbit version published on PyPI. The dev branch keeps the same three geometry families and the same Siddon projector/backprojector class names, while replacing closed-form circular geometry parameters by explicit per-view trajectory tensors. In practice this means that circular, spiral, sinusoidal, saddle, random, and user-defined trajectories all flow through the same projector kernels without changing the CUDA implementation. The trajectory-array analytical helpers on dev track the overhaul merged into main 1.2.10/1.2.11; the current public main branch has since advanced to the 1.3.4 release line.

The report focuses on four contributions. First, it gives a single notation for the main circular-orbit parameterization and the dev arbitrary-trajectory generalization. Second, it records the differentiable projector/backprojector pair implemented in the source, following the discrete ray-tracing lineage of [4]. Third, it derives the analytical reconstruction helpers actually used in the examples and tests, including cosine pre-weighting, Parker weighting, ramp filtering, angular integration weights, and weighted voxel-driven backprojection [5,6]. Fourth, it summarizes the repository-level validation logic encoded in adjoint, gradcheck, reconstruction-accuracy, offset, smoke, and benchmark tests, positioning diffct relative to differentiable CT libraries such as PYRO-NN and recent geometry-learning work [7,8].

Code and availability.

The source repository is hosted at https://github.com/sypsyp97/diffct, and stable releases are distributed through PyPI at https://pypi.org/project/diffct/ (pip install diffct).

2. Related Work

diffct sits at the intersection of four established CT lines of work. The first is classical analytical reconstruction. Parallel-beam FBP, fan-beam short-scan weighting, and cone-beam FDK remain the baseline reference formulas for the circular-orbit setting [1,5,6]. The second is discrete projector design, where Siddon’s plane-intersection traversal remains one of the most widely used ray-driven baselines [4]. The third is footprint-based projector design, especially the separable-footprint family in [9], which showed how voxel footprints can be approximated efficiently by separable trapezoidal or rectangular detector-domain kernels. The fourth is differentiable CT software, where libraries such as PYRO-NN and LEAP established the value of pairing GPU projectors with modern deep-learning workflows [7,10,11].

The dev branch also benefits from a broader arbitrary-trajectory literature. Exact cone-beam inversion is governed by completeness conditions and data sufficiency constraints, classically associated with Tuy and Smith [12,13]. Recent work such as DRACO revisits these ideas in a differentiable setting and makes arbitrary-orbit reconstruction a first-class optimization target [14]. diffct occupies a narrower but useful engineering position in this landscape: it exposes explicit per-view geometry tensors, exact discrete adjoints for Siddon-style projectors, and a practical analytical reconstruction path that shares geometry conventions with PyTorch-based reconstruction workflows [3,8].

3. Branch Lineage and Scope

3.1. Why the Report Is Centered on dev

The dev branch is the most complete mathematical target for a full report because it captures the generalized Siddon-based interface for circular and non-circular trajectories. The main branch remains the reference line for the circular-orbit SF family. The public API on dev re-exports the same high-level Siddon operator families:

• ParallelProjectorFunction / ParallelBackprojectorFunction,
• FanProjectorFunction / FanBackprojectorFunction,
• ConeProjectorFunction / ConeBackprojectorFunction,
• geometry generators for parallel, fan, and cone beam scanning,
• analytical helpers for FBP/FDK pipelines adapted to trajectory arrays.

The report therefore uses the dev Siddon/trajectory-array interface as the primary mathematical model and treats the main branch as the circular-orbit specialization and SF reference line.

3.2. Main Versus Dev

Table 1 summarizes the operational distinction between the two branches.

The unresolved gap is equally important. The main branch contains separable-footprint backends that depend on closed-form circular geometry and therefore remain branch-local. In 2D fan beam, backend="sf" denotes a separable trapezoidal detector footprint. In 3D cone beam, backend="sf_tr" uses a transaxial trapezoid with an axial rectangle, while backend="sf_tt" uses trapezoids in both detector directions, following the separable-footprint taxonomy in [9]. The dev branch keeps the same autograd structure and trajectory-array analytical helpers but omits these SF kernels because their footprint construction is still specialized to circular geometry.

The branch split also creates explicit compatibility costs. Existing code written against scalar angles/sid/sdd/backend signatures must now construct per-view geometry tensors before calling the projector classes. The separable-footprint backends available on main are absent on dev. The detector-grid convention on the dev analytical helpers uses $(k-N/2)\Delta u$, whereas the current main circular closed-form helpers still center detector bins with $(k-(N-1)/2)\Delta u$, so copying parameters across branches without adjustment introduces a half-cell shift for even detector sizes. The report therefore uses the dev equations as the global interface and treats the main SF family as an important circular-orbit specialization that still needs to be documented explicitly.

4. Unified Geometry Model

4.1. Discrete Detector Coordinates

The trajectory-array analytical helpers on dev use the detector-cell convention

$$u_k=\left(k-{\displaystyle \frac{N_u}{2}}\right)\Delta u+u_0,v_{\ell }=\left(\ell -{\displaystyle \frac{N_v}{2}}\right)\Delta v+v_0,$$

(1)

where $N_u$ and $N_v$ denote detector sizes, $\Delta u$ and $\Delta v$ denote detector spacing, and $u_0,v_0$ are optional offsets. This section follows the dev analytical-helper convention with the factor $N_u/2$ rather than $(N_u-1)/2$; the current main circular closed-form helpers still use $(N_u-1)/2$.

4.2. Parallel-Beam Parameterization

For a view index n, the dev branch represents 2D parallel geometry by a unit ray direction ${\mathit{d}}_n\in {\mathbb{R}}^2$, a detector origin ${\mathit{o}}_n\in {\mathbb{R}}^2$, and a detector axis ${\mathit{e}}_{u,n}\in {\mathbb{R}}^2$ perpendicular to ${\mathit{d}}_n$. Detector sample k is anchored at

$${\mathit{r}}_{n,k}\left(0\right)={\mathit{o}}_n+u_k{\mathit{e}}_{u,n},{\mathit{r}}_{n,k}\left(t\right)={\mathit{r}}_{n,k}\left(0\right)+t{\mathit{d}}_n,{\parallel {\mathit{d}}_n\parallel }_2=1.$$

(2)

The circular-orbit helper used as the main specialization is

$$\mathit{d}\left(\theta \right)=\left[\begin{array}{c}cos\theta \\ sin\theta \end{array}\right],\mathit{o}\left(\theta \right)=\left[\begin{array}{c}-d_{\det }sin\theta \\ \phantom{-}{d}_{\det }cos\theta \end{array}\right],{\mathit{e}}_u\left(\theta \right)=\left[\begin{array}{c}-sin\theta \\ \phantom{-}cos\theta \end{array}\right].$$

(3)

Here $\theta \in [0,\pi )$ is the non-redundant ray-direction angle used by the implementation, and $d_{\det }$ is the detector-to-isocenter distance. $\theta$ differs from the classical Radon-transform normal angle by a fixed $\pi /2$ reparameterization.

4.3. Fan-Beam Parameterization

For fan beam CT the implementation stores a source position ${\mathit{S}}_n\in {\mathbb{R}}^2$, a detector center ${\mathit{C}}_n\in {\mathbb{R}}^2$, and a detector axis ${\mathit{e}}_{u,n}\in {\mathbb{R}}^2$. The detector point associated with cell k is

$${\mathit{D}}_{n,k}={\mathit{C}}_n+u_k{\mathit{e}}_{u,n},{\widehat{\mathit{d}}}_{n,k}={\displaystyle \frac{{\mathit{D}}_{n,k}-{\mathit{S}}_n}{\parallel {\mathit{D}}_{n,k}-{\mathit{S}}_n{\parallel }_2}}.$$

(4)

The circular-orbit specialization used on the stable branch is

$$\mathit{S}\left(\beta \right)=\left[\begin{array}{c}-\mathrm{sid}sin\beta \\ \phantom{-}\mathrm{sid}cos\beta \end{array}\right],\mathit{C}\left(\beta \right)=\left[\begin{array}{c}(\mathrm{sdd}-\mathrm{sid})sin\beta \\ -(\mathrm{sdd}-\mathrm{sid})cos\beta \end{array}\right],{\mathit{e}}_u\left(\beta \right)=\left[\begin{array}{c}cos\beta \\ sin\beta \end{array}\right].$$

(5)

Here $\beta \in [0,2\pi )$ is the view angle, $\mathrm{sid}$ denotes the source-to-isocenter distance, and $\mathrm{sdd}$ denotes the source-to-detector distance. Per-view quantities reuse the subscript n (e.g. ${\mathrm{sid}}_n$ below).

4.4. Cone-Beam Parameterization

For cone beam CT each view uses a source ${\mathit{S}}_n\in {\mathbb{R}}^3$, detector center ${\mathit{C}}_n\in {\mathbb{R}}^3$, and orthonormal detector basis vectors ${\mathit{e}}_{u,n},{\mathit{e}}_{v,n}\in {\mathbb{R}}^3$. Detector coordinate $(k,\ell )$ maps to

$${\mathit{D}}_{n,k,\ell }={\mathit{C}}_n+u_k{\mathit{e}}_{u,n}+v_{\ell }{\mathit{e}}_{v,n},{\widehat{\mathit{d}}}_{n,k,\ell }={\displaystyle \frac{{\mathit{D}}_{n,k,\ell }-{\mathit{S}}_n}{\parallel {\mathit{D}}_{n,k,\ell }-{\mathit{S}}_n{\parallel }_2}}.$$

(6)

The circular-orbit 3D helper is

$$\mathit{S}\left(\beta \right)=\left[\begin{array}{c}-\mathrm{sid}sin\beta \\ \phantom{-}\mathrm{sid}cos\beta \\ 0\end{array}\right],\mathit{C}\left(\beta \right)=\left[\begin{array}{c}(\mathrm{sdd}-\mathrm{sid})sin\beta \\ -(\mathrm{sdd}-\mathrm{sid})cos\beta \\ 0\end{array}\right],{\mathit{e}}_u\left(\beta \right)=\left[\begin{array}{c}cos\beta \\ sin\beta \\ 0\end{array}\right],{\mathit{e}}_v=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right].$$

(7)

4.5. Non-Circular Trajectories

The arbitrary-trajectory generators in diffct.geometry instantiate several useful families. For example, the spiral cone-beam trajectory uses

$${\beta }_n={\beta }_0+{\displaystyle \frac{2\pi N_{\mathrm{turns}}}{N_{\mathrm{views}}}}n,z_n\in \left[-{\displaystyle \frac{z_{max}}{2}},{\displaystyle \frac{z_{max}}{2}}\right],$$

(8)

where $N_{\mathrm{views}}$ is the number of views in the scan, $N_{\mathrm{turns}}$ is the number of full helical turns, ${\beta }_0$ is the start angle, and $z_n$ is the axial source offset of view n bounded by $z_{max}$.

$${\mathit{S}}_n^{\mathrm{spiral}}=\left[\begin{array}{c}-\mathrm{sid}sin{\beta }_n\\ \phantom{-}\mathrm{sid}cos{\beta }_n\\ z_n\end{array}\right],{\mathit{C}}_n^{\mathrm{spiral}}=\left[\begin{array}{c}(\mathrm{sdd}-\mathrm{sid})sin{\beta }_n\\ -(\mathrm{sdd}-\mathrm{sid})cos{\beta }_n\\ z_n\end{array}\right].$$

(9)

The sinusoidal and saddle variants modulate the source-to-isocenter distance and axial position,

$${\mathrm{sid}}_n=\mathrm{sid}+Asin\left(\omega {\beta }_n\right),$$

(10)

$$z_n=A_{z}cos\left(2{\beta }_n\right),{\mathrm{sid}}_n=\mathrm{sid}+A_{r}sin\left(2{\beta }_n\right),$$

(11)

where $A,A_r$ are radial modulation amplitudes, $A_z$ is the axial modulation amplitude, and $\omega \in {\mathbb{Z}}^{+}$ is the angular frequency of the sinusoidal wobble. and the custom generators accept user-defined source paths or ray fields. The mathematical consequence is simple: every downstream operator only sees the per-view geometry tensors, so the projector definition is trajectory agnostic.

This abstraction is stronger than the analytical theory that supports it. The forward and adjoint Siddon operators only require valid ray definitions, so they transfer immediately from circular to non-circular orbits. Exact cone-beam inversion, by contrast, depends on data-completeness conditions of the Tuy–Smith type [12,13]. diffct therefore treats arbitrary trajectories as a first-class capability for projection, adjoint backprojection, and iterative optimization, while its analytical fan/cone reconstruction helpers should be interpreted as geometry-aware practical reconstructions whose exactness is guaranteed in the classical circular settings and whose usefulness outside that regime is empirical.

Figure 1. Unified geometry parameterization used by diffct. The parallel-beam panel shows the ray direction $\mathit{d}\left(\theta \right)$, detector origin $\mathit{o}\left(\theta \right)$, and detector axis ${\mathit{e}}_u$; the fan- and cone-beam panels highlight the source $\mathit{S}\left(\beta \right)$, detector plane, and basis vectors ${\mathit{e}}_u,{\mathit{e}}_v$; the fourth panel previews circular, sinusoidal, spiral, and saddle trajectories available through the diffct.geometry generators.

Figure 1. Unified geometry parameterization used by diffct. The parallel-beam panel shows the ray direction $\mathit{d}\left(\theta \right)$, detector origin $\mathit{o}\left(\theta \right)$, and detector axis ${\mathit{e}}_u$; the fan- and cone-beam panels highlight the source $\mathit{S}\left(\beta \right)$, detector plane, and basis vectors ${\mathit{e}}_u,{\mathit{e}}_v$; the fourth panel previews circular, sinusoidal, spiral, and saddle trajectories available through the diffct.geometry generators.

5. Differentiable Forward and Adjoint Operators

5.1. Continuous Forward Model

Let $f:{\mathbb{R}}^d\to \mathbb{R}$ denote the attenuation field, with $d=2$ for the parallel/fan-beam operators and $d=3$ for the cone-beam operator, and let $t\in \mathbb{R}$ be the ray arc-length parameter. The projector implemented in diffct computes line integrals along rays determined by the geometry tensors. For parallel beam,

$$P_{\mathrm{par}}\left[f\right](n,k)={\int }_{\mathbb{R}}f\left({\mathit{r}}_{n,k}\left(t\right)\right)dt,$$

(12)

with ${\mathit{r}}_{n,k}\left(t\right)$ from Equation (2). For fan beam,

$$P_{\mathrm{fan}}\left[f\right](n,k)={\int }_0^{\infty }f\left({\mathit{S}}_n+t{\widehat{\mathit{d}}}_{n,k}\right)dt,$$

(13)

and for cone beam,

$$P_{\mathrm{cone}}\left[f\right](n,k,\ell )={\int }_0^{\infty }f\left({\mathit{S}}_n+t{\widehat{\mathit{d}}}_{n,k,\ell }\right)dt.$$

(14)

These equations are branch-independent. The difference between main and dev lies only in how the trajectory tensors are produced.

5.2. Siddon Traversal with Cell-Constant Integration

The CUDA kernels implement a Siddon-style voxel traversal [4]. Consider a ray $\mathit{r}\left(t\right)=\mathit{o}+t\mathit{d}$ passing through an axis-aligned voxel grid. Intersections with grid planes satisfy

$$t_m^{\left(q\right)}={\displaystyle \frac{b_m^{\left(q\right)}-o_m}{d_m}},m\in \{x,y,z\},$$

(15)

where $\mathit{o}$ and $\mathit{d}$ are the ray origin and direction, $o_m$ and $d_m$ are their m-th components, q indexes the grid planes along axis m, and $b_m^{\left(q\right)}$ is the corresponding grid-plane coordinate. Sorting all valid intersection parameters yields segment boundaries $t_0<t_1<\dots <t_M$ inside the volume, where M is the number of segments the ray cuts across the grid. On each segment the attenuation field is taken to be constant at the value of the single pixel (2D) or voxel (3D) that the segment traverses, so the continuous line integral is approximated as

$$P\left[f\right]\left(\mathrm{ray}\right)\approx {\displaystyle \sum _{m=0}^{M-1}}\Delta t_m{f}_{\mathrm{cell}\left(m\right)},\Delta t_m=t_{m+1}-t_m,$$

(16)

where $\mathrm{cell}\left(m\right)$ denotes the discrete grid index of the cell traversed by segment m and $\Delta t_m$ is the exact chord length of the segment inside that cell. This matches the inner loop of the diffct projector kernels, which step from one grid-plane crossing to the next and accumulate d_image[iy, ix] * seg_len (2D) or d_vol[iz, iy, ix] * seg_len (3D) without any sub-cell interpolation. The forward/adjoint kernel pair is therefore matched by construction: the adjoint scatters the incoming ray-domain gradient into the same cell with the same segment length, so the discrete transpose identity holds up to floating-point tolerance (verified by test_adjoint_inner_product.py).

5.3. Adjoint Identity and Autograd

diffct exposes the projector and backprojector as paired torch.autograd.Function classes. The design target is the linear-adjoint relation

$$\langle P_{\Theta }x,y\rangle =\langle x,P_{\Theta }^{\top }y\rangle ,$$

(17)

where $\Theta$ collects all geometry tensors. The test suite checks Equation (17) for parallel, fan, and cone beam projectors. The backward pass of a forward projector therefore executes the pure adjoint backprojection kernel rather than an analytical FBP/FDK kernel.

This distinction is essential. The analytical reconstruction operators introduced later are not the same as the autograd adjoints because analytical FBP/FDK require geometry-dependent weighting and scaling. By contrast, the autograd backward uses the exact discrete transpose of the corresponding forward kernel. If

$$\mathit{y}=P_{\Theta }\mathit{x},\mathit{z}=P_{\Theta }^{\top }\mathit{s},$$

(18)

then the implemented backward rules are

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathcal{L}}{\partial \mathit{x}}}& =P_{\Theta }^{\top }{\displaystyle \frac{\partial \mathcal{L}}{\partial \mathit{y}}},{\displaystyle \frac{\partial \mathcal{L}}{\partial \Theta }}=0\mathrm{for}\mathrm{the}\mathrm{current}\mathrm{wrappers},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathcal{L}}{\partial \mathit{s}}}& =P_{\Theta }{\displaystyle \frac{\partial \mathcal{L}}{\partial \mathit{z}}},{\displaystyle \frac{\partial \mathcal{L}}{\partial \Theta }}=0\mathrm{for}\mathrm{the}\mathrm{current}\mathrm{wrappers}.\hfill \end{array}$$

(20)

The zero geometry gradient in is an implementation fact, not a mathematical impossibility: the source and detector tensors are saved in the autograd context and reused during backward, but the backward methods return None for geometry arguments. In other words, diffct currently differentiates with respect to image and sinogram variables while treating geometry as fixed runtime input. This is consistent with the repository’s present optimization examples, which solve for reconstructions while holding trajectories fixed.

For a least-squares inverse problem

$$\widehat{f}=arg\underset{f}{min}{\parallel P_{\Theta }f-p\parallel }_2^2+\lambda R\left(f\right),$$

(21)

the data-term gradient is

$${\nabla }_f{\displaystyle \frac{1}{2}}{\parallel P_{\Theta }f-p\parallel }_2^2=P_{\Theta }^{\top }(P_{\Theta }f-p).$$

(22)

This is exactly the quantity propagated by PyTorch through the diffct projector classes. At kernel level, the forward path accumulates line integrals ray by ray, whereas the backward path revisits the same traversed segments and scatters each incoming gradient into the single pixel or voxel that the segment traverses, weighted by the segment length $\Delta t_m$ of Equation (16). Because many rays contribute to the same image cell, the backward kernels use cuda.atomic.add to preserve the discrete transpose relation.

Algorithm 1 Autograd contract of the diffct projector wrappers

Require:   image or volume $\mathit{x}$, geometry tensors $\Theta$, output sinogram shape Ensure:   sinogram $\mathit{y}=P_{\Theta }\mathit{x}$; on Backward, image gradient ${\mathit{g}}_x=P_{\Theta }^{\top }{\mathit{g}}_y$ 1: procedureForward($\mathit{x},\Theta$) 2:     $\mathit{x}\leftarrow float32\left(\mathit{x}\right).\mathtt{contiguous}\left(\right).\mathtt{cuda}\left(\right)$ 3:     $\tilde{\mathit{x}},\tilde{\Theta }\leftarrow \mathtt{cuda}\_\mathtt{as}\_\mathtt{array}\left(\mathit{x}\right),\mathtt{cuda}\_\mathtt{as}\_\mathtt{array}(\Theta )$    ▹ zero-copy Numba views 4:     allocate $\mathit{y}\leftarrow \mathbf{0}$ 5:     KernelForward($\tilde{\mathit{x}},\tilde{\Theta },\mathit{y}$)     ▹ Siddon gather, current torch stream 6:     $\mathtt{ctx}.\mathtt{save}(\Theta )$ 7:     return $\mathit{y}$ 8: end procedure 9: procedureBackward(${\mathit{g}}_y$)     ▹${\mathit{g}}_y=\partial \mathcal{L}/\partial \mathit{y}$ 10:     $\Theta \leftarrow \mathtt{ctx}.\mathtt{load}\left(\right)$ 11:     allocate ${\mathit{g}}_x\leftarrow \mathbf{0}$ 12:     KernelAdjoint(${\mathit{g}}_y,\tilde{\Theta },{\mathit{g}}_x$)     ▹ cuda.atomic.add scatter 13:     return $({\mathit{g}}_x,\mathtt{None},\dots ,\mathtt{None})$    ▹ no gradient w.r.t. $\Theta$ 14: end procedure

Figure 2. Ray-driven forward gather and atomic-scatter adjoint. The forward kernel accumulates cell-constant voxel contributions along each ray (each traversed cell value weighted by its exact chord length) and produces the sinogram entry y. The adjoint kernel revisits the same traversed cells and scatters the incoming gradient back with cuda.atomic.add, which realizes the discrete transpose relation $\langle P_{\Theta }x,y\rangle =\langle x,P_{\Theta }^{\top }y\rangle$ tested by test_adjoint_inner_product.py. The third panel sketches how the paired kernels plug into the PyTorch autograd contract for $x\to y$ and $\partial \mathcal{L}/\partial y\to \partial \mathcal{L}/\partial x$.

Figure 2. Ray-driven forward gather and atomic-scatter adjoint. The forward kernel accumulates cell-constant voxel contributions along each ray (each traversed cell value weighted by its exact chord length) and produces the sinogram entry y. The adjoint kernel revisits the same traversed cells and scatters the incoming gradient back with cuda.atomic.add, which realizes the discrete transpose relation $\langle P_{\Theta }x,y\rangle =\langle x,P_{\Theta }^{\top }y\rangle$ tested by test_adjoint_inner_product.py. The third panel sketches how the paired kernels plug into the PyTorch autograd contract for $x\to y$ and $\partial \mathcal{L}/\partial y\to \partial \mathcal{L}/\partial x$.

6. Analytical Reconstruction in diffct

6.1. Detector Coordinates and Angular Weights

The analytical module builds FBP/FDK pipelines out of small helpers. Detector coordinates reuse the dev convention in Equation (1). For open angle lists and short scans, angular integration weights are trapezoidal:

$$\Delta {\beta }_1={\displaystyle \frac{{\beta }_2-{\beta }_1}{2}},\Delta {\beta }_n={\displaystyle \frac{({\beta }_n-{\beta }_{n-1})+({\beta }_{n+1}-{\beta }_n)}{2}},\Delta {\beta }_{N_{\mathrm{views}}}={\displaystyle \frac{{\beta }_{N_{\mathrm{views}}}-{\beta }_{N_{\mathrm{views}}-1}}{2}}.$$

(23)

Here ${\left\{{\beta }_n\right\}}_{n=1}^{N_{\mathrm{views}}}$ are the per-view angles and $\Delta {\beta }_n$ is the trapezoidal angular integration weight attached to view n. For periodic full circular fan or cone scans, the implementation closes the angle loop before applying the extra redundancy factor $MM1/2$. For a short scan with Parker weighting, the open-interval trapezoidal weights in Equation (23) are kept and redundancy is handled by the Parker window itself.

6.2. Ramp Filter and Apodization

Given a weighted sinogram $p_w$, the filtered projection is computed by a 1D FFT-based ramp filter along the detector axis:

$$q={\mathcal{F}}^{-1}\left\{\left|2\pi \xi \right|W\left(\xi \right)\mathcal{F}\left[p_w\right]\right\}{\displaystyle \frac{1}{\Delta s}},$$

(24)

where $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ are the 1-D Fourier transform and its inverse along the detector axis, $\xi$ is the physical Fourier frequency obtained from the normalized torch.fft.fftfreq grid after division by the detector spacing, $W\left(\xi \right)$ is an apodization window, and $\Delta s$ is the detector spacing along the filtered axis (either $\Delta u$ or $\Delta v$ depending on the pipeline). The explicit factor $1/\Delta s$ is the implementation’s conversion from sample frequency to physical units. The apodization window $W\left(\xi \right)$ supports Ram–Lak, Hann, Hamming, cosine, and Shepp–Logan choices. Using the normalized frequency $\nu =2\left|\xi \right|\in [0,1]$, the implemented window functions are

$$\begin{array}{cc}\hfill W_{\mathrm{Ram}-\mathrm{Lak}}\left(\nu \right)& =1,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill W_{\mathrm{Hann}}\left(\nu \right)& ={\displaystyle \frac{1}{2}}\left(1+cos\left(\pi \nu \right)\right),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill W_{\mathrm{Hamming}}\left(\nu \right)& =0.54+0.46cos\left(\pi \nu \right),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill W_{\mathrm{Cosine}}\left(\nu \right)& =cos\left({\displaystyle \frac{\pi \nu }{2}}\right),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill W_{\mathrm{Shepp}-\mathrm{Logan}}\left(\nu \right)& ={\displaystyle \frac{sin(\pi \nu /2)}{\pi \nu /2}}.\hfill \end{array}$$

(29)

6.3. Parallel-Beam FBP

Parallel beam requires no source-dependent distance weighting. Let $u_n\left(\mathit{x}\right)$ denote the detector coordinate obtained by projecting voxel $\mathit{x}$ onto the detector axis:

$$u_n\left(\mathit{x}\right)=\left(\mathit{x}-{\mathit{o}}_n\right)·{\mathit{e}}_{u,n}.$$

(30)

After angular weighting and ramp filtering, diffct reconstructs via

$$B_{\mathrm{par}}\left[q\right]\left(\mathit{x}\right)={\displaystyle \frac{1}{2\pi }}{\displaystyle \sum _{n=1}^{N_{\mathrm{views}}}}{q}_n\left(u_n\left(\mathit{x}\right)\right).$$

(31)

The factor $1/\left(2\pi \right)$ is applied inside parallel_weighted_backproject, so the wrapper already returns an amplitude-calibrated reconstruction.

6.4. Fan-Beam FBP

For fan beam geometry diffct computes a per-view detector normal ${\mathit{n}}_n$ by rotating ${\mathit{e}}_{u,n}$ by ${90}^{\circ }$ and orienting it so that $({\mathit{C}}_n-{\mathit{S}}_n)·{\mathit{n}}_n>0$. It then defines

$${\mathrm{sdd}}_n=({\mathit{C}}_n-{\mathit{S}}_n)·{\mathit{n}}_n,{\mathrm{sid}}_n=(-{\mathit{S}}_n)·{\mathit{n}}_n,U_n\left(\mathit{x}\right)=(\mathit{x}-{\mathit{S}}_n)·{\mathit{n}}_n.$$

(32)

Here ${\mathrm{sdd}}_n$ and ${\mathrm{sid}}_n$ are the signed source-to-detector and source-to-isocenter distances along ${\mathit{n}}_n$, and $U_n\left(\mathit{x}\right)$ is the signed distance from the source to the plane orthogonal to ${\mathit{n}}_n$ through $\mathit{x}$; it reduces to ${\mathrm{sid}}_n$ when $\mathit{x}$ is at the isocenter. For the classical circular short-scan specialization with constant $\mathrm{sdd}$, the fan angle and cosine weight are

$$\gamma \left(u\right)=arctan\left({\displaystyle \frac{u}{\mathrm{sdd}}}\right),w_{cos}^{\mathrm{fan}}\left(u\right)={\displaystyle \frac{\mathrm{sdd}}{\sqrt{{\mathrm{sdd}}^2+u^2}}}=cos\gamma \left(u\right).$$

(33)

For a short scan the Parker window used by the code is

$$w_P(\beta ,\gamma )=\left\{\begin{array}{cc}{sin}^2\left({\displaystyle \frac{\pi }{4}}{\displaystyle \frac{\beta }{{\gamma }_{max}-\gamma }}\right),\hfill & 0\le \beta <2({\gamma }_{max}-\gamma ),\hfill \\ 1,\hfill & 2({\gamma }_{max}-\gamma )\le \beta \le \pi -2\gamma ,\hfill \\ {sin}^2\left({\displaystyle \frac{\pi }{4}}{\displaystyle \frac{\pi +2{\gamma }_{max}-\beta }{{\gamma }_{max}+\gamma }}\right),\hfill & \pi -2\gamma <\beta \le \pi +2{\gamma }_{max},\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(34)

with ${\gamma }_{max}={max}_u\left|\gamma \left(u\right)\right|$. This Parker window is the circular short-scan helper used by the code; the general trajectory-array fan pipeline on dev reuses without elevating Parker weighting to a generic arbitrary-trajectory claim. A voxel $\mathit{x}$ is projected onto the detector plane through the source:

$${\mathit{H}}_n\left(\mathit{x}\right)={\mathit{S}}_n+{\displaystyle \frac{{\mathrm{sdd}}_n}{U_n\left(\mathit{x}\right)}}\left(\mathit{x}-{\mathit{S}}_n\right),u_n\left(\mathit{x}\right)=\left({\mathit{H}}_n\left(\mathit{x}\right)-{\mathit{C}}_n\right)·{\mathit{e}}_{u,n}.$$

(35)

diffct then applies the weighted FBP gather

$$B_{\mathrm{fan}}\left[q\right]\left(\mathit{x}\right)={\displaystyle \frac{\overline{{\mathrm{sdd}}_n}}{2\pi \overline{{\mathrm{sid}}_n}}}{\displaystyle \sum _{n=1}^{N_{\mathrm{views}}}}{\left({\displaystyle \frac{{\mathrm{sid}}_n}{U_n\left(\mathit{x}\right)}}\right)}^2{q}_n\left(u_n\left(\mathit{x}\right)\right),$$

(36)

where the overline denotes the mean over views. For a circular orbit Equation (36) reduces to the textbook factor

$$U(\beta ;x,y)=\mathrm{sid}+xsin\beta -ycos\beta ,u(\beta ;x,y)=\mathrm{sdd}{\displaystyle \frac{xcos\beta +ysin\beta }{U(\beta ;x,y)}}.$$

(37)

6.5. Cone-Beam FDK

The cone-beam extension follows the same structure with detector normal

$${\mathit{n}}_n={\sigma }_n{\displaystyle \frac{{\mathit{e}}_{u,n}\times {\mathit{e}}_{v,n}}{\parallel {\mathit{e}}_{u,n}\times {\mathit{e}}_{v,n}{\parallel }_2}},$$

(38)

where ${\sigma }_n$ is chosen so that $({\mathit{C}}_n-{\mathit{S}}_n)·{\mathit{n}}_n>0$. The detector-normal distances are

$${\mathrm{sdd}}_n=({\mathit{C}}_n-{\mathit{S}}_n)·{\mathit{n}}_n,{\mathrm{sid}}_n=(-{\mathit{S}}_n)·{\mathit{n}}_n,U_n\left(\mathit{x}\right)=(\mathit{x}-{\mathit{S}}_n)·{\mathit{n}}_n.$$

(39)

The cosine pre-weight is

$$w_{cos}^{\mathrm{cone}}(u,v)={\displaystyle \frac{{\mathrm{sdd}}_n}{\sqrt{{\mathrm{sdd}}_n^2+u^2+v^2}}}.$$

(40)

The detector intersection of a voxel $\mathit{x}$ is

$${\mathit{H}}_n\left(\mathit{x}\right)={\mathit{S}}_n+{\displaystyle \frac{{\mathrm{sdd}}_n}{U_n\left(\mathit{x}\right)}}\left(\mathit{x}-{\mathit{S}}_n\right),$$

(41)

with

$$u_n\left(\mathit{x}\right)=\left({\mathit{H}}_n\left(\mathit{x}\right)-{\mathit{C}}_n\right)·{\mathit{e}}_{u,n},v_n\left(\mathit{x}\right)=\left({\mathit{H}}_n\left(\mathit{x}\right)-{\mathit{C}}_n\right)·{\mathit{e}}_{v,n}.$$

(42)

The weighted FDK gather used by the code is

$$B_{\mathrm{cone}}\left[q\right]\left(\mathit{x}\right)={\displaystyle \frac{\overline{{\mathrm{sdd}}_n}}{2\pi \overline{{\mathrm{sid}}_n}}}{\displaystyle \sum _{n=1}^{N_{\mathrm{views}}}}{\left({\displaystyle \frac{{\mathrm{sid}}_n}{U_n\left(\mathit{x}\right)}}\right)}^2{q}_n\left(u_n\left(\mathit{x}\right),v_n\left(\mathit{x}\right)\right).$$

(43)

For a circular orbit this becomes

$$u(\beta ;x,y)=\mathrm{sdd}{\displaystyle \frac{xcos\beta +ysin\beta }{U(\beta ;x,y)}},v(\beta ;x,y,z)=\mathrm{sdd}{\displaystyle \frac{z}{U(\beta ;x,y)}}.$$

(44)

diffct applies the ramp filter row-wise along the u direction, matching the standard FDK construction. For non-circular trajectories, Equation (43) is a principled implementation-level extension rather than an exact closed-form inversion.

6.6. Main-Branch Separable-Footprint Family

The main branch adds a second circular-orbit backprojection family alongside the Siddon-based weighted gather path. These SF options are exposed through backend selectors on the stable fan- and cone-beam APIs, remain restricted to circular geometry, and inherit detector-footprint logic from [9] together with the matched-adjoint implementation style used in LEAP [10,11]. For a 2D fan-beam pixel $\mathit{x}$, let the four projected pixel corners induce the ordered detector abscissae

$$u_{min}\le u_{\mathrm{lo}}\le u_{\mathrm{hi}}\le u_{max}.$$

(45)

The detector footprint is then approximated by the trapezoid

$$h_{n,\mathit{x}}^{\mathrm{fan}}\left(u\right)=\left\{\begin{array}{cc}{\displaystyle \frac{u-u_{min}}{u_{\mathrm{lo}}-u_{min}}},\hfill & u_{min}\le u<u_{\mathrm{lo}},\hfill \\ 1,\hfill & u_{\mathrm{lo}}\le u\le u_{\mathrm{hi}},\hfill \\ {\displaystyle \frac{u_{max}-u}{u_{max}-u_{\mathrm{hi}}}},\hfill & u_{\mathrm{hi}}<u\le u_{max},\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(46)

The SF fan-beam gather used in main can therefore be summarized as

$$B_{\mathrm{fan}}^{\mathrm{SF}}\left[q\right]\left(\mathit{x}\right)=c_{\mathrm{fan}}^{\mathrm{SF}}\sum _n{\gamma }_{n,\mathit{x}}^{\mathrm{fan}}{\int }_{\mathbb{R}}{h}_{n,\mathit{x}}^{\mathrm{fan}}\left(u\right)q_n\left(u\right)du,$$

(47)

where the implementation uses a LEAP-inspired chord-weighted factor

$${\gamma }_{n,\mathit{x}}^{\mathrm{fan}}={\displaystyle \frac{{\mathrm{sid}}_n}{U_n\left(\mathit{x}\right)}}{\ell }_{\varphi ,n}\left(\mathit{x}\right),$$

(48)

with ${\ell }_{\varphi ,n}\left(\mathit{x}\right)$ denoting the in-plane chord through the unit pixel along the source-to-pixel direction. Relative to the Siddon analytical gather in Equation (36), the distinguishing feature is the explicit integration over detector footprint support rather than point sampling at $u_n\left(\mathit{x}\right)$.

For 3D cone beam, the SF family stays separable:

$$h_{n,\mathit{x}}^{\mathrm{cone}}(u,v)=h_{n,\mathit{x}}^u\left(u\right)h_{n,\mathit{x}}^v\left(v\right).$$

(49)

The transaxial factor $h_{n,\mathit{x}}^u$ is trapezoidal in both sf_tr and sf_tt. The axial factor $h_{n,\mathit{x}}^v$ is rectangular for sf_tr and trapezoidal for sf_tt. The corresponding gather is

$$B_{\mathrm{cone}}^{\mathrm{SF}}\left[q\right]\left(\mathit{x}\right)=c_{\mathrm{cone}}^{\mathrm{SF}}\sum _n{\gamma }_{n,\mathit{x}}^{\mathrm{cone}}{\int }_{{\mathbb{R}}^2}{h}_{n,\mathit{x}}^{\mathrm{cone}}(u,v)q_n(u,v)dudv,$$

(50)

with the chord-weighted per-view factor implemented in the code base summarized by

$${\gamma }_{n,\mathit{x}}^{\mathrm{cone}}={\ell }_{\varphi ,n}\left(\mathit{x}\right)\sqrt{1+{\left({\displaystyle \frac{v_n\left(\mathit{x}\right)}{{\mathrm{sdd}}_n}}\right)}^2}.$$

(51)

This is the mathematical reason the SF backends in main require dedicated kernels and smaller launch blocks: they integrate filtered detector values against voxel footprints, whereas the Siddon analytical path samples a single detector coordinate per voxel and view.

7. Software Architecture and CUDA Realization

The library is deliberately modular. diffct.projectors contains the PyTorch-facing autograd functions, diffct.geometry generates trajectory tensors, diffct.analytical builds analytical FBP/FDK pipelines, diffct.kernels stores the Numba CUDA kernels, diffct.utils handles device and stream management, and diffct.constants exposes low-level launch parameters and numerical constants. This separation is one of the clearest architectural improvements of the dev branch over the older monolithic diffct.differentiable interface.

Algorithm 2 Shared analytical preprocessing and the branch-specific backprojection stage used in diffct

Require:   sinogram $\mathit{p}$, geometry tensors $\Theta$, short-scan flag $s\in \{0,1\}$, backprojection selector b, with $b=\mathrm{siddon}$ on dev and $b\in \{\mathrm{siddon},\mathrm{sf}\}$ on circular main fan/cone paths Ensure:   reconstruction $\widehat{\mathit{f}}$ 1: compute detector coordinates $u_k,v_{\ell }$    ▹ Equation (1) 2: ${\{\Delta {\beta }_n\}}_{n=1}^{N_{\mathrm{views}}}\leftarrow$TrapezoidalWeights($\left\{{\beta }_n\right\}$)     ▹ Equation (23) 3: if geometry $\in \{\mathrm{fan},\mathrm{cone}\}$ then 4:     $p_{n,k(,\ell )}\leftarrow w_{cos}\left(u_k(,v_{\ell })\right)·p_{n,k(,\ell )}$    ▹ Equations (33) and (40) 5: end if 6: if$s=1$then 7:     $p_{n,k}\leftarrow w_P({\beta }_n,\gamma \left(u_k\right))·p_{n,k}$    ▹ Parker window, Equation (34) 8: end if 9: $q_{n,·}\leftarrow {\mathcal{F}}^{-1}\left\{W\left(\xi \right)\left|\xi \right|\mathcal{F}\left\{p_{n,·}\right\}\right\}/\Delta u$    ▹ Equation (24) 10: $q_{n,·}\leftarrow \Delta {\beta }_n·q_{n,·}$    ▹ angular integration weight 11: $\widehat{\mathit{f}}\leftarrow \mathbf{0}$ 12: for each voxel ${\mathit{x}}_i$ do 13:     if $b=\mathrm{siddon}$ then 14:         ${\widehat{f}}_i\leftarrow {\sum }_{n=1}^{N_{\mathrm{views}}}{w}_{n,i}{\tilde{q}}_n\left(u_n\left({\mathit{x}}_i\right),v_n\left({\mathit{x}}_i\right)\right)$    ▹ point-sampled gather 15:     else 16:         ${\widehat{f}}_i\leftarrow {\sum }_{n=1}^{N_{\mathrm{views}}}{\sum }_{k,\ell }{\phi }_{n,i}(u_k,v_{\ell })q_{n,k,\ell }$    ▹ separable-footprint, Equation (51) 17:     end if 18: end for 19: $\widehat{\mathit{f}}\leftarrow C_{\mathrm{geom}}·\widehat{\mathit{f}}$    ▹ geometry-dependent global scale 20: return$\widehat{\mathit{f}}$

7.1. PyTorch and Numba Execution Model

diffct combines PyTorch’s tensor/autograd ecosystem [3] with Numba’s CUDA JIT compilation model [15]. The forward wrappers first normalize all inputs to contiguous float32 tensors on the active CUDA device. They then call TorchCUDABridge.tensor_to_cuda_array, which exposes zero-copy Numba views through cuda.as_cuda_array. This detail matters because the projector kernels operate directly on PyTorch-owned GPU memory rather than on copied buffers.

Stream semantics are equally important. diffct retrieves the active PyTorch stream and wraps it with a cached numba.cuda.external_stream handle. As a result, Numba kernel launches respect the same dependency order as surrounding PyTorch code. The implementation therefore avoids an entire class of subtle synchronization bugs in which PyTorch and Numba would otherwise enqueue work on unrelated streams.

7.2. Kernel Families and Decorators

The code base uses two explicit numerical kernel families.

• Siddon projector and pure adjoint kernels are decorated with cuda.jit(cache=True, fastmath=True) to maximize throughput for ray traversal and the cell-constant segment accumulation.
• Analytical FBP/FDK gather kernels are decorated with cuda.jit(cache=True, fastmath=False) because the Fourier-convention constants and the geometry-dependent magnification weights amplify round-off errors more strongly than the projector kernels do.

This design explains why the code distinguishes between “autograd backprojection” and “analytical backprojection” even when both appear to map projection data back into image space. The former is the exact discrete transpose used by backward; the latter is a weighted reconstruction operator intended to approximate inversion.

The default thread-block sizes are $(16,16)$ for 2D kernels and $(8,8,8)$ for 3D kernels. The circular-orbit SF cone-beam kernels on main use smaller launch blocks than the default 3D Siddon kernels because the footprint integrals are more register-intensive. Cone-beam tensors are also permuted internally from $(D,H,W)$ to $(W,H,D)$ for memory-coalesced access in the 3D gather path and then permuted back before returning to the user.

7.3. Gradient Transport and Atomic Accumulation

The gradient path follows the same hardware design as the forward path. Each projector forward method stores geometry tensors with ctx.save_for_backward and keeps output-shape metadata in the autograd context. During backward, diffct allocates a zero-filled gradient image or volume, converts the incoming PyTorch gradient tensor into a Numba CUDA view, and launches the matched adjoint kernel on the same CUDA stream.

The backward kernel mirrors the cell-constant forward in Equation (16): forward kernels accumulate the traversed cell value weighted by the chord length $\Delta t_m$ along each ray, while backward kernels scatter the incoming ray-domain gradient back into the same cell with the same $\Delta t_m$. Because multiple rays update the same pixel or voxel, the scatter step uses cuda.atomic.add. This is the implementation mechanism that makes the discrete adjoint relation in Equation (17) true to machine precision up to the expected floating-point tolerance.

8. Experimental Protocol and Repository Validation

8.1. Reference Phantoms and Example Setups

The analytical accuracy tests use Shepp–Logan phantoms in 2D and 3D. The public main branch additionally ships synthetic real-data pipelines and a walnut CBCT example. The circular phantom examples documented in this report follow three reference configurations:

• parallel-beam FBP on a $256\times 256$ phantom with 360 views and 512 detectors,
• fan-beam FBP on a $256\times 256$ phantom with 360 views, 600 detectors, $\mathrm{sid}=500$, and $\mathrm{sdd}=800$,
• cone-beam FDK on a ${128}^3$ phantom with 360 views, a $256\times 256$ detector, $\mathrm{sid}=600$, and $\mathrm{sdd}=900$.

The non-circular iterative examples then replace circular trajectories by sinusoidal, wobbling, elliptical, spiral, saddle, and custom figure-eight variants while reusing the same projector classes.

8.2. Circular Analytical Reconstructions

The trajectory-array analytical helper line on dev inherits the amplitude correction introduced on main 1.2.10/1.2.11 by routing the filtered sinogram through the dedicated weighted gather kernels. The change log records the corresponding raw reconstruction errors listed in Table 2.

8.3. Iterative and Arbitrary-Trajectory Reconstructions

The iterative examples solve Equation (21) with a learned correction variable and an AdamW optimizer. Parallel-beam experiments compare sinusoidal and wobbling trajectories. Fan-beam experiments compare sinusoidal and elliptical source paths. Cone-beam experiments compare spiral, sinusoidal, saddle, and custom figure-eight trajectories. The dev branch ships generator scripts and plotting utilities for all of these families, so reproducing the iterative experiments requires no additional modelling effort beyond swapping the trajectory factory.

8.4. Correctness Tests

The repository-level validation strategy is unusually strong for a compact research library. On the public dev branch, the default test suite contains 62 tests; the current public main branch contains 71 tests, while an additional 27 benchmark cases are kept opt-in. The test categories summarized in Table 3 describe the shared validation themes across the two lines.

One representative test threshold is the cone-beam analytical reconstruction guard: the 3D Shepp–Logan FDK test requires a raw RMSE below $0.1$, a center-slice RMSE below $0.1$, bounded axial-profile deviation, finite values everywhere, and a reconstruction maximum below $1.5$. This combination protects against the exact classes of bugs that commonly appear in analytical CT implementations: wrong Fourier constants, omitted cosine weights, omitted ${(\mathrm{sid}/U)}^2$ factors, or detector-axis scaling mistakes.

8.5. Performance Benchmarking

The benchmark suite measures three kernels for each geometry: forward projector, pure adjoint backprojector, and the full analytical pipeline. The analytical benchmarks include the complete sequence of cosine weighting, ramp filtering, angular weighting, and weighted gather backprojection because this is the user-visible reconstruction cost. In the current public benchmark scripts, the measured kernels are the Siddon forward path, the pure adjoint path, and the full analytical pipeline; the SF backends are part of the stable main API but are not benchmarked separately in the published benchmark files.

9. Discussion

diffct now occupies a clean position in the CT software landscape. It is small enough to audit mathematically, but broad enough to support analytical reconstruction, differentiable inversion, and arbitrary-trajectory experimentation inside a single code base. The current dev line is especially attractive for research because it preserves the classical circular-orbit formulas as special cases while making trajectory tensors first-class objects, even though the public main branch has already moved beyond the 1.2.11 analytical-sync point.

The main limitation is equally clear. The weighted fan and cone analytical operators implemented in dev are a geometry-aware extension of FBP/FDK based on detector-normal distances and per-view mean scaling. For circular orbits they collapse to the familiar closed forms. For non-circular orbits they provide a practical and principled reconstruction heuristic, but they do not claim an exact inversion formula beyond the classical setting. A second limitation is branch asymmetry: the main branch still owns the SF backends, so the highest-fidelity detector-footprint model is currently available only for circular trajectories. A third limitation is gradient scope. The present autograd wrappers propagate gradients through images and sinograms, while geometry tensors are runtime inputs with no returned gradients. A fourth limitation is hardware: the entire library assumes CUDA-capable execution and uses PyTorch plus Numba as its software substrate.

The branch split also reveals a natural roadmap. The main branch proves that separable-footprint projectors are useful in the circular setting. The dev branch proves that the generic per-view tensor API is the right abstraction for arbitrary trajectories. Bridging those two lines of work will require a deeper reformulation of footprint models for detector planes that are no longer aligned with a closed-form circular geometry.

10. Conclusions

This report rewrites diffct from the code upward. The stable main branch contributes the circular-orbit baseline, while the dev branch turns the Siddon-based projector/backprojector interface into a general arbitrary-trajectory CT framework. The resulting mathematical picture is compact: projector kernels implement Siddon traversal with cell-constant segment integration; autograd uses exact discrete adjoints with fixed geometry tensors; analytical reconstruction uses explicit cosine weighting, Parker short-scan weighting in the circular fan-beam specialization, FFT-based ramp filtering, angular integration weights, and dedicated weighted gather backprojectors. With this structure in place, diffct is already well positioned for deep-learning integration, gradient-based reconstruction, and arbitrary-trajectory CT research.

Detector coordinates.

$$u_k=\left(k-{\displaystyle \frac{N_u}{2}}\right)\Delta u+u_0,v_{\ell }=\left(\ell -{\displaystyle \frac{N_v}{2}}\right)\Delta v+v_0.$$

(A1)

Forward operators.

$$\begin{array}{cc}\hfill P_{\mathrm{par}}\left[f\right](n,k)& ={\int }_{\mathbb{R}}f\left({\mathit{o}}_n+u_k{\mathit{e}}_{u,n}+t{\mathit{d}}_n\right)dt,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill P_{\mathrm{fan}}\left[f\right](n,k)& ={\int }_0^{\infty }f\left({\mathit{S}}_n+t{\widehat{\mathit{d}}}_{n,k}\right)dt,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill P_{\mathrm{cone}}\left[f\right](n,k,\ell )& ={\int }_0^{\infty }f\left({\mathit{S}}_n+t{\widehat{\mathit{d}}}_{n,k,\ell }\right)dt.\hfill \end{array}$$

(55)

Siddon discretization.

$$P\left[f\right]\left(\mathrm{ray}\right)\approx {\displaystyle \sum _{m=0}^{M-1}}\Delta t_m{f}_{\mathrm{cell}\left(m\right)}.$$

(A2)

Adjoint relation and inverse problem.

$$\langle P_{\Theta }x,y\rangle =\langle x,P_{\Theta }^{\top }y\rangle ,\widehat{f}=arg\underset{f}{min}{\parallel P_{\Theta }f-p\parallel }_2^2+\lambda R\left(f\right),{\nabla }_f\mathcal{L}\left(f\right)=2P_{\Theta }^{\top }(P_{\Theta }f-p)+\lambda \nabla R\left(f\right).$$

(A3)

$${\displaystyle \frac{\partial \mathcal{L}}{\partial x}}=P_{\Theta }^{\top }{\displaystyle \frac{\partial \mathcal{L}}{\partial y}},{\displaystyle \frac{\partial \mathcal{L}}{\partial s}}=P_{\Theta }{\displaystyle \frac{\partial \mathcal{L}}{\partial z}},{\displaystyle \frac{\partial \mathcal{L}}{\partial \Theta }}=0\mathrm{in}\mathrm{the}\mathrm{current}\mathrm{implementation}.$$

(A4)

Angular weights and ramp filter.

$$\Delta {\beta }_1={\displaystyle \frac{{\beta }_2-{\beta }_1}{2}},\Delta {\beta }_n={\displaystyle \frac{({\beta }_n-{\beta }_{n-1})+({\beta }_{n+1}-{\beta }_n)}{2}},\Delta {\beta }_N={\displaystyle \frac{{\beta }_N-{\beta }_{N-1}}{2}},$$

(A5)

$$q={\mathcal{F}}^{-1}\left\{\left|2\pi \xi \right|W\left(\xi \right)\mathcal{F}\left[p_w\right]\right\}{\displaystyle \frac{1}{\Delta s}}.$$

(A6)

Cosine weights.

$$w_{cos}^{\mathrm{fan}}\left(u\right)={\displaystyle \frac{\mathrm{sdd}}{\sqrt{{\mathrm{sdd}}^2+u^2}}},w_{cos}^{\mathrm{cone}}(u,v)={\displaystyle \frac{\mathrm{sdd}}{\sqrt{{\mathrm{sdd}}^2+u^2+v^2}}}.$$

(A7)

Weighted analytical backprojection.

$$\begin{array}{cc}\hfill B_{\mathrm{par}}\left[q\right]\left(\mathit{x}\right)& ={\displaystyle \frac{1}{2\pi }}\sum _n{q}_n\left(u_n\left(\mathit{x}\right)\right),\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill B_{\mathrm{fan}}\left[q\right]\left(\mathit{x}\right)& ={\displaystyle \frac{\overline{{\mathrm{sdd}}_n}}{2\pi \overline{{\mathrm{sid}}_n}}}\sum _n{\left({\displaystyle \frac{{\mathrm{sid}}_n}{U_n\left(\mathit{x}\right)}}\right)}^2{q}_n\left(u_n\left(\mathit{x}\right)\right),\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill B_{\mathrm{cone}}\left[q\right]\left(\mathit{x}\right)& ={\displaystyle \frac{\overline{{\mathrm{sdd}}_n}}{2\pi \overline{{\mathrm{sid}}_n}}}\sum _n{\left({\displaystyle \frac{{\mathrm{sid}}_n}{U_n\left(\mathit{x}\right)}}\right)}^2{q}_n\left(u_n\left(\mathit{x}\right),v_n\left(\mathit{x}\right)\right).\hfill \end{array}$$

(64)

Separable-footprint backprojection on main.

$$h_{n,\mathit{x}}^{\mathrm{fan}}\left(u\right)=\left\{\begin{array}{cc}{\displaystyle \frac{u-u_{min}}{u_{\mathrm{lo}}-u_{min}}},\hfill & u_{min}\le u<u_{\mathrm{lo}},\hfill \\ 1,\hfill & u_{\mathrm{lo}}\le u\le u_{\mathrm{hi}},\hfill \\ {\displaystyle \frac{u_{max}-u}{u_{max}-u_{\mathrm{hi}}}},\hfill & u_{\mathrm{hi}}<u\le u_{max},\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(A8)

$$\begin{array}{c}\hfill B_{\mathrm{fan}}^{\mathrm{SF}}\left[q\right]\left(\mathit{x}\right)=c_{\mathrm{fan}}^{\mathrm{SF}}\sum _n{\gamma }_{n,\mathit{x}}^{\mathrm{fan}}\int h_{n,\mathit{x}}^{\mathrm{fan}}\left(u\right)q_n\left(u\right)du,\end{array}$$

(66)

$$\begin{array}{c}\hfill B_{\mathrm{cone}}^{\mathrm{SF}}\left[q\right]\left(\mathit{x}\right)=c_{\mathrm{cone}}^{\mathrm{SF}}\sum _n{\gamma }_{n,\mathit{x}}^{\mathrm{cone}}\int \int h_{n,\mathit{x}}^u\left(u\right)h_{n,\mathit{x}}^v\left(v\right)q_n(u,v)dudv.\end{array}$$

(67)

Table A1. Crosswalk from the report’s mathematical notation to the current public diffct API, including projector classes, trajectory generators, analytical helpers, and the main-only SF backend names.

Report object	Code object
$P_{\mathrm{par}},P_{\mathrm{fan}},P_{\mathrm{cone}}$	ParallelProjectorFunction, FanProjectorFunction, ConeProjectorFunction
$P_{\mathrm{par}}^{\top },P_{\mathrm{fan}}^{\top },P_{\mathrm{cone}}^{\top }$	ParallelBackprojectorFunction, FanBackprojectorFunction, ConeBackprojectorFunction
Circular trajectories	circular_trajectory_2d_parallel, circular_trajectory_2d_fan, circular_trajectory_3d
Non-circular trajectories	spiral_trajectory_3d, sinusoidal_trajectory_*, saddle_trajectory_3d, random_trajectory_3d, custom_trajectory_*
$u_k,v_{\ell }$ detector coordinates	detector_coordinates_1d
$\Delta {\beta }_n$ angular weights	angular_integration_weights
$w_{cos}^{\mathrm{fan}},w_{cos}^{\mathrm{cone}}$	fan_cosine_weights, cone_cosine_weights
$w_P$ Parker window	parker_weights
Ramp filter q	ramp_filter_1d
$B_{\mathrm{par}},B_{\mathrm{fan}},B_{\mathrm{cone}}$	parallel_weighted_backproject, fan_weighted_backproject, cone_weighted_backproject
Main-branch SF backend selectors	backend="sf" for fan beam, backend="sf_tr" and backend="sf_tt" for cone beam

Table A1. Crosswalk from the report’s mathematical notation to the current public diffct API, including projector classes, trajectory generators, analytical helpers, and the main-only SF backend names.

Report object	Code object
$P_{\mathrm{par}},P_{\mathrm{fan}},P_{\mathrm{cone}}$	ParallelProjectorFunction, FanProjectorFunction, ConeProjectorFunction
$P_{\mathrm{par}}^{\top },P_{\mathrm{fan}}^{\top },P_{\mathrm{cone}}^{\top }$	ParallelBackprojectorFunction, FanBackprojectorFunction, ConeBackprojectorFunction
Circular trajectories	circular_trajectory_2d_parallel, circular_trajectory_2d_fan, circular_trajectory_3d
Non-circular trajectories	spiral_trajectory_3d, sinusoidal_trajectory_*, saddle_trajectory_3d, random_trajectory_3d, custom_trajectory_*
$u_k,v_{\ell }$ detector coordinates	detector_coordinates_1d
$\Delta {\beta }_n$ angular weights	angular_integration_weights
$w_{cos}^{\mathrm{fan}},w_{cos}^{\mathrm{cone}}$	fan_cosine_weights, cone_cosine_weights
$w_P$ Parker window	parker_weights
Ramp filter q	ramp_filter_1d
$B_{\mathrm{par}},B_{\mathrm{fan}},B_{\mathrm{cone}}$	parallel_weighted_backproject, fan_weighted_backproject, cone_weighted_backproject
Main-branch SF backend selectors	backend="sf" for fan beam, backend="sf_tr" and backend="sf_tt" for cone beam