Numerical Backward Ray-Tracing

• Numerical backward ray-tracing is a computational method that reconstructs field values and transfer functions by tracing measurement rays backward through a domain.
• It employs rigorous operator frameworks, specialized discretization with basis functions, and adjoint formulations to enhance the accuracy of applications like tomographic reconstruction and radiative transfer.
• Innovative algorithmic strategies leverage GPU parallelism, adaptive geometric predicates, and memory-efficient loops to ensure high fidelity and scalable performance.

Numerical backward ray-tracing refers to the suite of computational algorithms designed to reconstruct field values, distributions, or transfer functions along paths traced from a measurement or observation point back through a physical or mathematical domain. The approach is most prominent in applications such as tomographic image reconstruction, radiative transfer, wave-optics rendering, and geodesic or optical path computation in complex media and spacetimes. Unlike forward ray tracing, where trajectories are synthesized from sources, backward (adjoint) ray-tracing quantifies, projects, or accumulates measurements back onto model parameters, optimizing computational efficiency and alignment with inverse problems. The methodology is characterized by rigorous geometric, algebraic, or variational formulations, specialized discretization and computational strategies, and explicit handling of accuracy, stability, and parallelization constraints.

1. Mathematical Foundations and Operator Frameworks

The mathematical structure of backward ray-tracing varies with application but universally relies on a geometric or linear-operator framework. In tomographic applications, the adjoint of the projection operator is essential. For a continuous $f:\mathbb{R}^2\to\mathbb{R}$ and x-ray transform

$$\mathcal P_\theta\{f\}(y) = \int_{-\infty}^{\infty} f(t\mathbf \theta + y\,\mathbf \theta^\perp )\,dt,$$

the adjoint acts as

$$(\mathcal P^*g)(\mathbf{x}) = \int g(\theta, \langle\mathbf{x},\mathbf{\theta}^\perp\rangle)\,\varphi(\mathbf{x}-\cdot)\,d(\theta, y),$$

where $\varphi$ denotes the generator or basis function in grid models. In numerical settings, this operator is discretized as a matrix transpose (for basis-based representations) or is represented algorithmically in phase-space or mesh-based constructs. Admissible basis functions include compact-support splines and box-splines, which permit exact computation of line integrals and corresponding projectors.

For mesh-based volumetric representations (e.g., tetrahedral meshes), the operator computes the contribution of each detector measurement to every mesh element traversed by the corresponding ray, proportional to the intersection length within each element. In radiative transfer and wave-optics, the operator formalism is extended to phase-space, sometimes involving Wigner distributions or generalized rays with phase-space Gaussian structure.

2. Numerical Implementations and Algorithmic Strategies

Backward ray-tracing algorithms are implemented through structured numerical loops, often optimized for high-throughput and memory efficiency. In basis-function tomographic settings (Haouchat et al., 26 Mar 2025), the algorithm traverses rays in geometric order and distributes measured projection values backward to overlapping basis-function grid cells. A critical element is the accurate computation of pre-weighted contributions (e.g., $\varphi_{\theta}(y)$) using closed-form expressions to avoid interpolation or quadrature errors.

For tetrahedral meshes (Biguri et al., 2019), the intersection of a ray with volumetric elements is determined using robust geometric predicates, such as minimally inflated Möller–Trumbore tests for face intersections, and lengths are computed accordingly. Algorithmic steps involve initializing ray parameters, walking through acceleration structures (such as R*-trees) for rapid entry-tetrahedron search, and executing a ray-walking loop that accumulates detector values into mesh-element weights via atomic memory operations.

In higher-level physics and computational geometry, backward ray-tracing extends to phase-space mappings (Jansen et al., 2023), where the problem becomes one of recursively mapping intervals or beams in direction space ('p') through system surfaces, using composition of propagation and reflection maps while preserving étendue.

Pseudocode structures across implementations favor minimal branch conditions, heapless recursion or loop-based iteration, and tight data packing to enable GPU parallelization or efficient multi-core scaling. The independence of rays in the backward direction allows for massive data-level parallelism.

3. Analytical Properties, Accuracy, and Error Sources

The primary analytical requirement is that the backward (adjoint) operator numerically matches the forward operation to within machine precision, ensuring no bias is introduced in inverse problem solutions such as conjugate-gradient reconstructions. For basis-function tomographic ray-tracing, this is numerically confirmed through the equality $$\langle H_\varphi \mathbf{c}, \mathbf{p} \rangle = \langle \mathbf{c}, H_\varphi^\top \mathbf{p} \rangle$$.

Error sources are limited to floating-point rounding if the basis-function projections are precomputed via closed-form expressions (Haouchat et al., 26 Mar 2025). In mesh settings, geometric degeneracies can induce missed intersections or double-counting unless robust adaptive tolerances are used—these are addressed via adaptive inflation of intersection thresholds ($\varepsilon$) and double-precision arithmetic (Biguri et al., 2019).

Accuracy in physical predictions (e.g., flux, power, or convergence profiles) is generally limited by grid or mesh resolution, the support of the basis functions, and the accuracy of geometric intersection routines. Empirical results demonstrate sub-percent or percent-level agreement with analytic or high-quality reference solutions provided these elements are carefully implemented.

4. Computational Complexity and Optimization Techniques

Computational complexity analysis reveals that backward ray-tracing can be made highly efficient, generally scaling as $\mathcal{O}(L\,K^2)$ per ray for $L$ traversed cells and $K$ basis-support radius (Haouchat et al., 26 Mar 2025), or $\mathcal{O}(P_\text{ray})$ for mesh-based settings with $P_\text{ray}$ mean number of intersected elements (Biguri et al., 2019). In practical high-resolution 3D problems, the overhead for more expressive basis functions or adaptive geometries remains only a small multiple over standard pixel- or voxel-based adjoints.

Optimization strategies include:

• Pre-computation and explicit storage of basis-projection functions for all discrete angles.
• Use of GPU-optimized libraries (e.g., Dr.Jit), with branch-free traversal and atomic memory operations to avoid synchronization overheads.
• Sparing use of branch conditions, typically employing integer comparison and neighbor-skipping logic to further reduce redundant computations.
• Explicit memory management for large, scattered-write accumulators—preferring atomic adds or block-level buffering in GPU implementations.
• Load balancing by assigning one ray per thread or process.

The algorithms exhibit near-perfect strong scaling due to the ray-independence of the adjoint operator, facilitating efficient deployment on large-scale parallel hardware.

5. Application Domains and Generalization

Backward ray-tracing algorithms are utilized in multiple domains:

• Medical and industrial tomographic image reconstruction, especially for non-standard grids or high-order representations (Haouchat et al., 26 Mar 2025, Biguri et al., 2019).
• Optical system analysis in both geometric and phase-space frameworks (Jansen et al., 2023), enabling rigorous accounting of étendue and accurate mapping through complex, possibly curved, optical elements.
• Radiative transfer in astrophysical environments, where computational speed and adaptive accuracy controls are realized via tree-based reverse ray-tracing over spatial hierarchies (Grond et al., 2019), supporting rigorous error bounds via opening-angle and absorption-depth criteria.
• Physical and wave-optics rendering, via backward tracing of generalized rays or Gaussian beams for full phase-space propagation and physically accurate image synthesis (Steinberg et al., 2023).
• Geodesic and field-equation integration in general relativity for cosmological or strong-field astrophysical predictions, employing adaptive ODE solvers and the Sachs/Jacobi formalism (Pihajoki et al., 2016, Koksbang et al., 2015, McDonald et al., 2023).

The generality of the approach is underpinned by the use of line-integral (ray-tracing) engines in both Euclidean and non-Euclidean geometries. For arbitrary projection geometries—parallel-beam, cone-beam, helical, non-uniform angle sets—the adjoint recursion proceeds unchanged, provided the geometric primitives are correctly defined. Phase-space and field-theoretic extensions accommodate wave effects, anisotropies, and polarization.

6. Performance Evaluation and Empirical Benchmarks

Empirical performance analyses across implementations highlight the practical competitiveness of backward ray-tracing approaches:

• For basis-function tomographic backprojections, degree-2 box-spline projectors and adjoints achieve 0.5–1 dB PSNR gains over standard pixel-based adjoints for coarse grids ($N\leq150$), and fewer projection angles are needed to attain comparable reconstruction quality (Haouchat et al., 26 Mar 2025).
• On GPU, voxel-based backprojectors achieve $~120$ ms per $1024^2$ projection of a $128^3$ volume, with tetrahedral mesh backprojectors being an order of magnitude slower per element but requiring $10$–$100\times$ fewer elements overall (Biguri et al., 2019).
• In phase-space mapped optical systems, backward mapping schemes (CBRM) achieve accuracies of $5\times10^{-4}$ to $10^{-3}$ in sub-second runtime, outperforming both straight-surface CBRM and naive Monte Carlo ray-tracing by orders of magnitude in both accuracy and speed (Jansen et al., 2023).
• Scaling performance is robust, with observed data lying between the $O(\log^2 N)$ (optimal) and $O(N^{1/3})$ (worst-case) bounds for tree-based algorithms in radiative transfer applications (Grond et al., 2019).
• In relativistic ray-tracing, solver implementations match machine-precision analytic and forward codes, with strict adherence to accuracy controlled by adaptive step tolerances (Pihajoki et al., 2016, McDonald et al., 2023).

7. Limitations, Assumptions, and Future Directions

Limitations are domain-specific:

• In basis-function approaches, finite support and discretization limit spatial resolution. No inherent quadrature or interpolation error exists if closed-form projectors are used, but errors may accumulate if support radii are insufficient.
• Numerical robustness requires high-precision arithmetic and adaptive tolerances to avoid topological errors, especially in mesh-based and curved-geometry environments.
• For astrophysical radiative transfer, no scattering or emission is included in basic algorithms; such additions can be implemented, but with increased computational complexity.
• Wave-optical backward path tracers are limited by the finite spatial spread of generalized rays; perfect locality is lost, and Kirchhoff diffraction requires additional sampling or modules. High-frequency caustics and sharply featured BSDFs remain challenging for sampling-based solutions.
• In strong-field relativistic or field-theoretic extensions, physical approximations (e.g., collisionless, monochromatic, weak-coupling regimes) delimit algorithmic applicability.

Future work addresses the extension of backward ray-tracing to full field-theoretic adjoint problems, mesh adaptivity for complex systems, efficient mixed representation support (mesh + phase-space), and continued improvement of accuracy and parallel scalability in increasingly heterogeneous hardware environments.