A pseudospectral method for solution of the radiative transport equation

The radiative transport equation accurately describes light transport in participating media such as biological tissues, though analytic solutions are known only for simple geometries. We present a pseudospectral technique to efficiently compute numerical solutions to the time-dependent transport equation in arbitrary three-dimensional heterogeneous media with anisotropic scattering. Our GPU-accelerated implementation of the technique is validated by comparison with a Monte-Carlo simulation, demonstrating excellent agreement.


INTRODUCTION
The radiative transport equation (RTE) provides an accurate description of light transport in biological tissues [1,2]. The RTE is a mono-energetic form of the Boltzmann transport equation, which models the propagation, attenuation and scattering of particles within participating media [3], and has been applied in numerous fields, including neutron transport, atmospheric physics.
The RTE, which can be derived directly from Maxwell's equations [4], describes the distribution of specific intensity, or radiance within a domain Ω, 1 c ∂ ∂t +ŝ · ∇ + µ t (r) φ(r,ŝ, t) = µ s (r) S N −1 p(ŝ ·ŝ )φ(r,ŝ , t)dŝ + q(r,ŝ, t) (1) where the change in radiance φ(r,ŝ, t) at a point r ∈ Ω in directionŝ at time t, is given as a balance of energy lost from an attenuation term µ t = µ a + µ s accounting for absorption (µ a ) and out-scattering (µ s ), inwards scattering from s to s, and any sources q(r,ŝ, t). In the biomedical imaging modalities of interest here, the Henyey-Greenstein phase function [5] is typically employed to describe the anisotropic scattering process, p(ŝ ·ŝ ) = 1 4π Boundary conditions for the RTE specify that φ(m,ŝ, t) = 0 forŝ ·n < 0, wheren is the outward normal to the boundary of the domain at m ∈ ∂Ω: that is to say that in the absence of sources, there is no component of the radiance inwards across the boundary. The high dimensionality of the field variable φ(r,ŝ, t) ∈ R N × S N −1 , and the presence of non-local spatial and angular operators, is such that analytically solutions to the RTE have been found only for homogeneous infinite and semi-infinite geometries [6], and layered media [7].
Numerical solutions to the RTE in arbitrary threedimensional geometries are typically sought by stochastic Monte-Carlo (MC) estimates. This approach converges to the solution in the limit of an infinite number of trials, and is inherently parallelisable. Despite the availability of a number of high performance codes, it remains a challenge to generate high quality estimates at distances far from a source, where statistical noise dominates. Furthermore, practical implementations require the use of variance reduction techniques which produce complex geometry dependent noise statistics which prevents modelling of the noise statistics as part of, e.g., an image reconstruction procedure.
Deterministic methods can be categorised by how they treat the angular discretisation.
• Discrete ordinate (S N ) methods discretise the angular basis over a set of vectors equally spaced over the unit-circle or sphere. S N methods suffer from ray effects and may require stabilisation which affects both the performance and accuracy of their solutions.
• Spherical harmonic discretisations (P N methods) expand the solution in the natural orthogonal basis over the sphere, and truncate the series to allow numerical solution. P N methods can suffer from wave-like distortion owing to the truncation of the angular basis, but perform well in highly scattering media such as biological tissue.
Discretisation in space can be achieved with the finiteelement method [8][9][10], though the large size of the resultant system matrices has limited the practical application of this method. Alternatively, finite-difference [11] and finite-volume [12] methods have been proposed, but the inherent derivative operations require fine discretisation and the use of higher order derivative stencils, leading to limitations in their application in large domains. The lack of accurate and computationally efficient deterministic methods often leads to the application of the arXiv:1801.06364v2 [physics.comp-ph] 8 Mar 2018 diffusion approximation (DA) to the RTE. The DA fails in regions of low scattering, and where the gradient of the energy density is large. Furthermore, the DA is parabolic (non-causal) in contrast with the hyperbolic RTE from which it is derived. Whilst these differences have until recently been tolerated in the context of diffuse optics, new developments in high-density, time-domain, and hybrid diffuse optical methods require the use of the RTE to achieve their full potential.

THE PSEUDO-SPECTRAL METHOD
Pseudo-spectral methods are typically applied to solve partial differential equations which contain pointwise potential operators. The technique exploits the O(N log N ) performance of the Fast Fourier Transform (FFT) to take spatial derivatives in the spectral domain (k-space) before applying an inverse transform to apply the potential operators in Cartesian space: thus, all operators are applied in the space in which they are diagonalised, and a far coarser spatial discretisation can be employed, limited only by the Nyquist criterion. Such techniques thus require far less memory and computational effort than their finite-difference or finite-element counterparts.
Pseudo-spectral schemes can be constructed for the RTE using both the S N and P N methods, however the elegance of our proposed technique is to exploit the fact that the scattering operator of the RTE is also diagonalised in the P N approximation, and that directional derivative term forms a recurrence relation over the spherical harmonics which is trivial to compute. We thus operate solely in the angular-spectral domain.
To derive our method we expand the radiance and source terms in a Spherical Harmonic basis, and employ the addition theorem for spherical harmonics to restate the phase function, Equations 3 through 5 are substituted into equation 1, and the result is projected into the spherical harmonic basis by taking the inner product withȲ m (ŝ). Following some algebra, this yields a set of coupled first order PDEs [13], where the terms α m = (l + m)(l − m) and β m = (l + m)(l + m − 1) are related to the Clebsch-Gordan coefficients [8,13].
To implement our pseudo-spectral scheme we perform an N dimensional spatial Fourier transform of the radiance distribution, and k = {k x , k y , k z } is the k-space wave-vector. Inserting this result into equation 6 results in a set of coupled ODEs, where F −1 denotes the inverse spatial Fourier transform, andψ

Discretisation, truncation, transient and steady-state solution
To solve the RTE numerically we form a P N approximation by closing the recurrence under the assumption that ψ m = 0 for l > N , leading to (N + 1) 2 coupled first-order ODEs. The radiance distribution is sampled in space and the continuous Fourier transform replaced with its discrete counterpart. Finally, we sample in the temporal domain and implement the temporal derivative using, for example, a first order explicit scheme, The use of an explicit temporal discretisation provides conditional stability, requiring that a suitable choice of time-step ∆t be employed. Equation 11 can be directly iterated to develop a timedomain solution. To compute a steady-state solution, this time-domain solution can be numerically integrated. Alternatively, the time-domain form can be integrated analytically, and the resulting operator used as part of an implicit iterative solver.

Boundary conditions
The primary limitation of pseudo-spectral methods is that they imply periodic boundary conditions (the exact form of which depends upon the choice of transform). In many applications a-periodic boundary conditions are enforced by the development of a perfectly matched layer (PML), which attenuates outgoing radiation before it reaches the boundary of the computational domain. Since the transport operator in the RTE is advective (that is to say that the action of transport on functions in the domain is to translate only along the directionsŝ), a PML can be implemented by the simple addition of an additional isotropic loss term outside of the real domain. Further, to adhere to the specified boundary conditions, which state that there is no radiation inwards across the boundary of the domain, we must set the scattering coefficient µ s = 0 in those regions which do not correspond the real domain interest. A full analysis of such boundary conditions was recently analysed in [14]. We thus solve, where µ b is the PML attenuation term, Ω c is the full computational domain, and Ω is the real domain of interest. The exact spatial distribution of µ b should be chosen to ensure that the solution is attenuated to a sufficient degree on the boundary of Ω c . Choice of µ b is straightforward since within Ω c the radiance decays at a rate rate exp(−µ b |r b |), where r b is the distance from the boundary in a given direction. It is also desirable for numerical reasons to choose a distribution which leads to a smooth decay of the radiance, in order to limit the introduction of higher spatial frequencies in the spectral representation.

RESULTS
The pseudo-spectral algorithm was implemented in the Julia programming language [15]. Optional GPU acceleration was implemented using the nVidia cuFFT library, and a custom kernel to apply the operator directly on the GPU. All simulations were executed on an nVidia Tesla K40 GPU driven by a 32-core Xeon workstation.
In the numerical demonstrations which follow, we in each case compare the results of our pseudo-spectral method with a Radiance Monte-Carlo technique we developed previously [16,17].

Transient solution without PML
The basic algorithm can be validated without the use of the PML by evaluating the time-domain internal radiance distribution over a period before the radiance distribution has propagated to the boundary.
We consider a cube 50mm on a side of background parameters µ a = 0.01mm −1 , µ s = 1.0mm −1 , and an anisotropy factor of g = 0.9. The region x ≤ 15mm has a higher scattering coefficient µ s = 5.0mm −1 . The domain is discretised into a 64 × 64 × 64 grid of points, and the pseudo-spectral method is used to solve the system under a 21-degree spherical harmonic (P 21 ) approximation.
The sharp onset of the scattering perturbation, and otherwise low scattering, provide an artificially challenging numerical scenario for the pseudo-spectral solver, promoting the propagation of higher angular harmonics and generating sharp features in the resultant fluence field.
A time-step of 0.1ps was chosen by inspection, and 400 steps were taken such that the initial radiance distribution, consisting of an Gaussian of isotropic angular profile, with 1/e diameter 1.5mm, has time to propagate close the boundary of the domain. The simulation took 70s to complete.
The same parameters were also employed in a Monte-Carlo solver of the RTE, in which 5 × 10 8 packets were propagated, with a kernel run-time of 269s. Figure 1 demonstrates the evolution of the radiance field during propagation over a line through the domain. Figure 2 shows the resultant fluence through the domain after 60ps, when the simulation was terminated. The initial field spreads outwards through the region of lowscattering close to the speed of light, forming a circular disc in the plane. As the radiance encounters the region of increased scattering a characteristic region of backscattering causes a peak at the interface, and the further transport is inhibited and the field becomes diffuse.
An exact comparison is challenging since the timedomain Monte-Carlo solution is inherently averaged spatially and temporally (whereas the pseudospectral method is a collocation technique, with an exact solution at points in space and time). Further reducing the spatial and temporal discretisation of the MC solution leads to an unacceptable reduction in SNR, highlighting the value of the proposed technique. Nonetheless, excellent agreement is found between the two solutions.

Steady-state internal fluence with PML
The effect of the PML can be evaluated by considering the steady-state fluence in the domain resulting from a given source distribution.
The real domain consists of a cuboid of 35.5mm on a side with homogeneous optical parameters µ a = 0.01mm −1 , µ s = 10.0mm −1 , and anisotropy factor g = 0.95. The real domain is embedded into a computation domain 50mm on a side with a PML attenuation term following a half-cosine profile rising from 0.01mm −1 at the end of the real domain to 1.0mm −1 at the boundary. The domain is discretised into a 64 × 64 × 64 grid of points, and the pseudo-spectral method is used to solve the system under a 7-degree spherical harmonic (P 7 ) ap-

proximation.
A time step of 0.04ps was chosen by inspection, and 25 × 10 3 time steps were integrated to ensure the domain had reached steady state conditions. The source term was an isotropic Gaussian distribution, with 1/e diameter 1.5mm, offset from the centre of the domain by 8mm in the x-axis. The simulation took 223s to complete.
The same parameters were also employed in a Monte-Carlo solver of the RTE, in which only the real domain was simulated and a matched boundary condition employed. The MC simulation employed 1 × 10 8 packets, with a kernel run-time of 16.1s. Figures 3 and 4 demonstrate the resultant steady state fluence for both the pseudo-spectral and Monte-Carlo techniques. In each case the characteristic logarithmic decay of the fluence is evident away from the source and boundary, with an increased decay towards the boundary. The index-matched boundary conditions imply a non-zero fluence at the boundary the domain, evident in each simulation.
Excellent agreement is shown between the pseudospectral method and the Monte-Carlo solution.

DISCUSSION & CONCLUSIONS
We have introduced a new technique to solve the radiative transport equation in arbitrary heterogeneous media. The technique has been validated against a stochastic Monte-Carlo solver, demonstrating excellent agreement.
The use of a pseudo-spectral method allows a compact representation of the radiance field with sparser sampling than required in a finite-difference or finiteelement scheme. The method is inherently parallelisable: whilst our implementation only uses a single GPU, the cuFFT library and CUDA kernels implement the algorithm can be readily extended to run across a number of devices with a minimum of inter-device communication. Our implementation is limited by the performance of the FFT, and thus the method demonstrates demonstrates O(N log N ) performance scaling.
The pseudospectral technique is particularly suited to biomedical imaging where inherent scattering leads to (angularly) smooth solutions which can be simulated with modest angular discretisations (e.g., N l ≤ 7). As we have demonstrated, the technique can readily support propagation in low-scattering regions when a higher degree approximation is employed. The absence of noise and deterministic performance is such that the technique could be employed as part of a conventional image reconstruction method in modalities such as diffuse optical tomography, quantitative photoacoustic tomography, and ultrasound modulated optical tomography.
The simplicity and flexibility of the technique suggests numerous avenues for further research. In particular, we propose acceleration of the technique with modified (temporal) finite-difference schemes, hybrid MCpseudospectral methods to treat regions of low scattering without excessive degrees of P N approximation, the implementation of reflecting boundary conditions, and custom iterative Krylov solvers to accelerate the computation of steady-state solutions.