Accurate and efficient Monte Carlo solutions to the radiative transport equation in the spatial frequency domain

We present an approach to solving the radiative transport equation (RTE) for layered media in the spatial frequency domain (SFD) using Monte Carlo (MC) simulations. This is done by obtaining a complex photon weight from analysis of the Fourier transform of the RTE. We also develop a modified shortcut method that enables a single MC simulation to efficiently provide RTE solutions in the SFD for any number of spatial frequencies. We provide comparisons between the modified shortcut method and conventional discrete transform methods for SFD reflectance. Further results for oblique illumination illustrate the potential diagnostic utility of the SFD phase-shifts for analysis of layered media.

Spatial frequency domain (SFD) imaging has shown great promise for noninvasive, quantitative, spectroscopic characterization of tissue composition, morphology, and physiology, which can be executed over large fields of view [1]. However, methods to extract optical (and physiological) properties from SFD images are currently limited to the diffusion approximation of the radiative transport equation (RTE) and discrete Fourier transformations of spatially resolved Monte Carlo (MC) simulations [1][2][3][4]. The availability of transport rigorous solvers in the SFD would enable accurate processing of SFD images as well as direct evaluations of approximate RTE solutions that are composed in the SFD without reliance on discrete Fourier inversion schemes [5,6].
SFD imaging illuminates the tissue surface with spatial patterns of light and measures the spatially modulated reflectance or fluorescence. While two-dimensional (2D) spatial patterns can be used, the use of one-dimensional (1D) sinusoids is most common. Current MC methods to obtain radiative transport estimates in the SFD first simulate the spatially resolved reflectance generated by a δ-function (point impulse) source and then apply a discrete Fourier transform (DFT). For 1D illumination patterns, SFD predictions can be obtained from a 1D DFT of a spatially resolved RTE solution. Importantly, use of a δfunction source at normal incidence provides radial symmetry. This symmetry allows use of the discrete Hankel transform (DHT), which is the most common approach to obtain predictions of SFD reflectance [1][2][3]. Discrete transform methods have been invaluable for the creation of lookup tables for inverse problems [2] and the study of layered tissue systems [3]. While discrete transforms provide a shortcut method for SFD solutions [7], such transforms are sensitive to the discrete sampling width and prone to inaccuracy and artifact.
Here we present the use of MC simulations to accurately and efficiently obtain transport rigorous estimates natively in the SFD without the use of discrete transforms. We accomplish this by obtaining a complex photon weight through analysis of the Fourier transform (FT) of the RTE. We then derive a modified shortcut method (MSM) that efficiently utilizes a single conventional MC simulation for a δ-function source to compute SFD estimates for multiple spatial frequencies. We verify this MSM with comparisons to discrete transform methods for a normally incident 1D sinusoidal modulated source. Moreover, we demonstrate its use for examining spatial phase-shifts resulting from SFD imaging of a thin layered medium using oblique illumination.
We begin with the source-free, time-independent RTE: (1) where L is the radiance, r the position vector, Ω the direction vector, p(Ω, Ω′) the singlescattering phase function, and μ s and μ a the local scattering and absorption coefficients, respectively. MC simulations of the RTE using continuous absorption weighting (CAW) employ photons that are attenuated continuously along each track according to Beer's Law, exp(−μ a ℓ j ), where ℓ j is the photon path length between the (j − 1) and j th collision. The photon weight following the j th track is The SFD representation of the radiance is given by the 2D FT over the transverse dimensions as , and Eq. (1) becomes (2) Besides losses due to absorption and scattering, this transformation reveals frequencydependent phase-shifts associated with radiative transport. The use of the FT over the transverse dimensions assumes constant optical properties over x and y and limits this method to layered media. From Eq. (2), a photon weight analogous to conventional CAW-MC is obtained, which includes the effects of spatial modulation. The photon weight is complex and expressed as (3) which can be formed as a triple product representing contributions from absorption and phase-shifts due to propagation along each frequency component. While Eq. (3) provides a photon deweighting scheme for the SFD, the photon weight is frequency dependent and requires updates along each photon trajectory for each spatial frequency considered.
Thus, we propose an MSM derived from Eq. (3) that exploits the fact that the frequencydependent phase accumulation is sensitive only to the net lateral displacement of the photon from the δ-function source as opposed to the detailed photon trajectory. This fact enables a single MC simulation to provide efficient RTE solutions for multiple spatial frequencies. To illustrate, we consider the f x contribution to the photon weight and express the direction cosine in terms of spatial co-ordinates: The cumulative effect of j collisions on the f x contribution to the photon weight is (4) By locating the δ-function source of the MC simulation at the origin, x 0 = y 0 = 0, and setting the initial photon weight to unity (W 0 = 1), Eq. (3) becomes (5) Equation (5) provides the complex weight of a detected photon for any spatial frequency within the MSM and can be used to bootstrap conventional CAW-MC simulations with a δfunction source. Moreover, for a normally incident 1D modulated source, the 2D complex weight reduces to exp(−2πif x x j ) exp(−2πif y y j ) = (2πf ρ ρ j ), where and is the zeroth-order Bessel function of the first kind. Thus, radial symmetry of the δ-function source reduces the calculation of the SFD photon weight to the kernel of the Hankel transform.
Functionals of the radiance can be formed from the frequency-dependent photon weights provided by the MSM. For example, the SFD reflectance for a 1D modulated light source normally incident on a homogeneous medium is (6) This formulation efficiently records for N photons the frequency-dependent contribution for any set of spatial frequencies following photon exit at location ρ n with total path length d n from the δ-function source. For a layered medium, the absorption term must be generalized to account for the photon path length and absorption coefficient within each layer according to Beer's Law.
We performed a conventional CAW-MC simulation with N = 10 7 photons for a δ-function source normally incident on a half-space with optical properties of , ℓ * = 1 mm, and g = 0.8 using the Henyey-Greenstein phase function. The exiting photon locations (x, y) and the photon weights due to absorption were saved in a database. The SFD reflectance is computed from this database using the various discrete transform and MSMs.
In Fig. 1, we examine the characteristics of the SFD reflectance provided by three equivalent MSMs for a normally incident 1D sinusoidally modulated source. We provide the percent difference of the SFD reflectance predicted using the Fourier kernels (f x , f y ) relative to the Hankel kernel (f ρ ). As the spatial frequency increases, the SFD reflectance from the Fourier kernels becomes mildly noisy relative to the Hankel kernel. This is because the Hankel kernel accounts for the phase accumulation due to photon propagation along both x and y directions, while the Fourier kernels account for the phase accumulation along only one of these directions. Thus, use of the Hankel kernel in conjunction with the MSM is preferred but requires the source to exhibit radial symmetry. Figure 2 compares the SFD reflectance computed via the MSM to that obtained from a discrete Fourier/Hankel transform using two different spatial sampling widths (Δ) in the discrete sum. As Δ decreases, the results of the discrete transform methods approach that of the MSM derived from the RTE. The DHT (f ρ ) outperforms the DFTs (f x , f y ) for the same spatial sampling width in comparison to the MSM. While this figure shows that the accuracy of the discrete transform methods improves as Δ is reduced, the MSM effectively represents the limit Δ → 0, as each detected photon contributes uniquely to the frequency-dependent reflectance without spatial binning.
Finally, we consider the case of oblique illumination of a two-layer medium and use the MSM to predict the amplitude and spatial phase-shift of the SFD reflectance [8]. For obliquely oriented 1D modulated sources propagating in the x-z plane, the Hankel kernel cannot be used and the f x Fourier kernel must be utilized instead. This kernel is applied to the results of a conventional MC simulation that employs an oblique δ-function source directed in the same orientation. We evaluated the SFD reflectance and phase-shift for a two-layer medium mimicking an epithelial tissue with a 200 μm superficial layer illuminated at a 45° angle. The medium has uniform optical absorption, anisotropy, and refractive index of μ a = 0.16 mm −1 , g = 0.9, and n = 1.4, respectively. We examined changes in the scattering coefficient of the superficial layer in the range μ s = 5-20 mm −1 in increments of 5 mm −1 while holding the scattering in the "stromal" layer constant at μ s = 20 mm −1 . Figure 3 shows that while scattering changes in the superficial layer have a minimal effect on the SFD reflectance, the effect on the phase-shift is significant and easily discernible. Moreover, such phase-shifts cannot be caused by changes in optical absorption [8]. This case shows how the MSM can be used for a priori design of optical diagnostic measurements in layered media.
In conclusion, we have introduced an MC method to obtain RTE solutions natively in the SFD for homogeneous and layered tissues. This was achieved by the identification of a complex photon weight via analysis of the FT of the RTE over the transverse dimensions. We derived an MSM to efficiently provide gold standard SFD radiative transport estimates for multiple spatial frequencies simultaneously. These estimates are obtained by processing the detected photon weights and exit locations from a conventional CAW-MC simulation. The method's performance was analyzed relative to discrete transform methods currently used to predict the SFD reflectance. We also provided results for a two-layer medium that shows sensitivity of the SFD phase-shift to small scattering changes in the thin superficial layer when subject to oblique illumination. These results represent the first transport rigorous solutions for oblique sources in the SFD. Moreover, the ability to obtain SFD reflectance measures using conventional CAW-MC simulations makes this method amenable for use with perturbation MC methods to effectively address important forward and inverse problems in optical imaging and diagnostics [9]. (Color online) Reflectance versus spatial frequency computed using f x , f y , and f ρ kernels in the 1D MSM. Relative differences between the f x and f y kernel results relative to the f ρ kernel are also shown. (Color online) SFD reflectance computed by 1D discrete transform using f x , f y , and f ρ kernels versus the MSM using the Hankel kernel. Colors denote the spatial widths Δ used in the discrete transforms. (Color online) Reflectance amplitude and phase versus spatial frequency for a two-layered "epithelial" tissue illuminated with an obliquely oriented 1D modulated source. Arrow indicates the effect of increasing top-layer scattering on the reflectance and phase-shift.