Analysis of single Monte Carlo methods for prediction of reflectance from turbid media

Michele Martinelli; Adam Gardner; David Cuccia; Carole Hayakawa; Jerome Spanier; Vasan Venugopalan

doi:10.1364/OE.19.019627

1. Introduction

Real-time optical analysis of turbid media requires rapid radiative transport solvers to predict and analyze remitted light signals for the non-invasive determination of optical scattering and absorption characteristics. While analytical techniques are often preferred for their simplicity and speed, they are often based on the diffusion approximation to the radiative transport equation (RTE). The limitations of the diffusion approximation are well known and result in accurate radiative transport estimates only when the reduced scattering coefficient, μ′_s = μ_s(1 – g), is much larger than the absorption coefficient μ_a and over spatial scales that are larger than the transport mean free path, l^* = (μ_a + μ′_s)⁻¹ [1, 2].

These limitations have motivated extensive efforts to solve the RTE, using either stochastic solvers such as Monte Carlo [3] or deterministic solvers via analytical [4] or numerical [5] approaches. RTE solutions provide accurate radiative transport predictions regardless of optical properties and on spatial scales approaching a single scattering length; ℓ_s = 1/μ_s ≈ 10–50μm for biological tissues in the red and near-infrared spectral regions. Unfortunately, these RTE solution approaches are computationally intensive and currently impractical for real-time application. For the case of homogeneous turbid media, several groups have proposed methods that utilize the results of a single MC simulation to provide rapid transport-rigorous RTE solutions for a class of related radiative transport problems.

Graaff and co-workers first investigated the use of a single MC (sMC) simulation to obtain results for turbid slabs with different sets of optical properties [10]. They demonstrated that the trajectories of unabsorbed photons in a single Monte Carlo simulation can be used to estimate the total reflectance and transmittance of turbid slabs for various single-scattering albedos, a = μ_s/(μ_a + μ_s) by recording the number of photon interactions within the simulated medium. Using this approach they predicted the spatially-resolved reflectance for different albedo values in cases for which the interaction coefficient, μ_t = μ_a + μ_s, was held constant.

Kienle and Patterson (KP) extended the capabilities of sMC simulations to provide the spatially- and temporally-resolved reflectance R(r,t) for a homogeneous semi-infinite medium with arbitrary optical properties [11]. This sMC method is based on scaling and re-weighting a discrete representation of the time-resolved reflectance provided by a reference MC simulation in a non-absorbing semi-infinite medium. In this approach the evaluation of R(r,t) with arbitrary optical properties, requires the interpolation of the discrete representation of R(r,t) in both the spatial and temporal dimensions. While powerful, the discrete representation of the reference R(r,t) and interpolation both introduce errors that are often amplified when the scaling is applied. Most recently, Alerstam and co-workers improved upon the KP sMC method by applying the scaling relationships on a photon-by-photon basis prior to spatial and temporal binning [13]. This approach requires storage of the radial exit position and total path length of each detected photon, which are individually reprocessed to evaluate the photon weights for media with different optical properties.

The motivation and objectives of this manuscript are two-fold. While sMC is widely used by the biomedical community [12, 14–17], a theoretical analysis that provides a clear line of derivation from the RTE to the scaling relations that form the basis of sMC has not been presented. Thus it is unclear whether differences between results obtained from conventional MC vs. sMC methods arise from stochastic variations, shortcomings of the scaling relationships and/or inaccuracies associated with the binning and interpolation used in the sMC implementation. To address these issues, our first goal is to determine to what extent the sMC scaling relations are transport-rigorous and to provide any ancillary conditions necessary for their validity. Second, in order to establish ‘best practices’, we aim to examine the impact of binning and interpolation strategies on the relative error/variance of the sMC predictions.

2. Theory

The temporally- and spatially-resolved reflectance R(r,t) is sufficient to derive all reflectance metrics typically used to characterize a turbid system, i.e., steady-state, temporal frequency-domain, and/or spatial frequency-domain reflectance. While conventional Monte Carlo simulations can be performed without explicit discretization of the phase space, their interpretation usually requires a mesh to be imposed on the phase space. Thus, the detected photons are placed into finite bins to visualize the simulation output. For a homogeneous half space occupying z ≥ 0 and possessing cylindrical symmetry, the spatially- and temporally-resolved reflectance, R(r,t)|_z=0⁻ at a given location (r ₀,t ₀) on the tissue surface can be estimated by

R (r_{0}, t_{0}) \approx \frac{1}{Δ A Δ t} \int_{t_{0}}^{t_{0} + Δ t} \int_{r_{0}}^{r_{0} + Δ r} \int_{𝕊^{2 -}} L (r, \hat{Ω}, t) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} 2 π r d r d t,

where L(r,

\hat{Ω}

,t) is the radiance,

\hat{Ω}

the unit direction vector, n̂ the outward pointing unit normal vector,

Δ A = π {(r_{0} + Δ r)}^{2} - π r_{0}^{2}

the surface area from which the photons are collected, and 𝕊²⁻ the solid angle interval corresponding to the upward-pointing hemisphere. Many photons will, in general, exit a given cylindrical ring within the same time interval. Thus, when these photons are gathered to obtain the value of R(r ₀,t ₀), a discretization error occurs because the quadrature interval has been represented by the nominal values (r ₀,t ₀). This error vanishes as the number of simulated photons increases without bound and as the integration intervals → 0.

Two scaling relationships are used in the sMC method to obtain spatially- and temporally-resolved reflectance estimates for a homogeneous half space. The first uses the reflectance of a ‘reference’ medium R_r(r,t) with μ_a = 0 and μ_s = μ_s,r to determine the reflectance R(r,t) from a second non-absorbing medium with arbitrary scattering coefficient μ_s [11]:

R (r, t) |_{μ_{s} \neq μ_{s, r}; μ_{a} = 0} = \frac{μ_{s}^{3}}{μ_{s, r}^{3}} R_{r} (\frac{μ_{s}}{μ_{s, r}} r, \frac{μ_{s}}{μ_{s, r}} t) .

The second scaling relationship relates the reflectance from a non-absorbing ‘reference’ medium (μ_a,r = 0) to a second medium with identical scattering but non-zero absorption μ_a:

R (r, t) |_{μ_{s} = μ_{s, r}; μ_{a} \neq 0} = R_{r} (r, t) exp (- μ_{a} c t)

where c is the velocity of light propagation in the medium. The analysis that follows will establish a rigorous foundation for these two scaling relationships.

We begin our analysis by considering the source-free RTE for a non-absorbing medium:

\frac{1}{c} \frac{\partial L (r, \hat{Ω}, t)}{\partial t} + \nabla L (r, \hat{Ω}, t) \cdot \hat{Ω} = - μ_{s} L (r, \hat{Ω}, t) + μ_{s} \int_{𝕊^{2}} p ({\hat{Ω}}^{'} \to \hat{Ω}) L (r, {\hat{Ω}}^{'}, t) d Ω^{'}

where L(r,

\hat{Ω}

,t) is the radiance at a given location r at at given time t in a given direction

\hat{Ω}

. c is the velocity of light propagation in the medium and is equivalent to the ratio between the speed of light in vacuum and the medium’s refractive index (= c ₀/n). μ_s and μ_a are the scattering and absorption coefficients, respectively, and p(

\hat{Ω}

′ →

\hat{Ω}

) is the single-scattering phase function, a probability density function for scattering from

\hat{Ω}

′ to

\hat{Ω}

. For boundary conditions, we assume a vanishing radiance at remote locations in the medium

L (r \to \infty, \hat{Ω}, t) \to 0,

and allow a possible refractive-index mismatch at the surface:

L (r, {\hat{Ω}}_{r}, t) = ℛ_{F} (\hat{n} \cdot {\hat{Ω}}_{i}) L (r, {\hat{Ω}}_{i}, t), r \in S, \hat{n} \cdot {\hat{Ω}}_{i} > 0,

where ℛ_F is the angular distribution of the Fresnel reflectance for unpolarized light and

\hat{Ω}

_r is the unit vector in the direction of the reflected ray for an incident direction

\hat{Ω}

_i. Thus ℛ_F (n̂ ·

\hat{Ω}

_i) represents the probability for internal reflection for a photon incident on an interface with a refractive index discontinuity.

We first express the RTE and boundary conditions in dimensionless form using the following variables:

{\begin{array}{l} τ = μ_{s} c t \\ ρ = μ_{s} r \\ \tilde{Λ} (ρ, \hat{Ω}, τ) = L (r, \hat{Ω}, t) / L_{0}, \end{array}

where τ is a dimensionless time, ρ a dimensionless position and

\tilde{Λ}

(ρ,

\hat{Ω}

,τ) a dimensionless radiance defined as the ratio between the radiance, L (r,

\hat{Ω}

, t), and a characteristic radiance L ₀. Substitution of these variables into the RTE yields:

\frac{\partial \tilde{Λ} (ρ, \hat{Ω}, τ)}{\partial τ} + \nabla \tilde{Λ} (ρ, \hat{Ω}, τ) \cdot \hat{Ω} = - \tilde{Λ} (ρ, \hat{Ω}, τ) + \int_{𝕊^{2}} p ({\hat{Ω}}^{'} \to \hat{Ω}) \tilde{Λ} (ρ, {\hat{Ω}}^{'}, τ) d Ω^{'} .

The two boundary conditions, Eqs. (5) and (6) become

\tilde{Λ} (ρ \to \infty, \hat{Ω}, τ) \to 0, and

\tilde{Λ} (ρ, {\hat{Ω}}_{r}, τ) = ℛ_{F} (\hat{n} \cdot {\hat{Ω}}_{i}) \tilde{Λ} (ρ, {\hat{Ω}}_{i}, τ), ρ \in S, \hat{n} \cdot {\hat{Ω}}_{i} > 0 .

Note that the non-dimensionalization does not alter the form of the boundary conditions. Having defined the boundary conditions, the diffuse reflectance is given by

R (r, t) = \int_{𝕊^{2 -}} L (r, \hat{Ω}, t) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} .

By analogy, a dimensionless reflectance is defined as

\tilde{R} (ρ, τ) = \int_{𝕊^{2 -}} \tilde{Λ} (ρ, \hat{Ω}, τ) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} .

To derive the first sMC scaling relation, Eq. (2), we consider the dimensionless analog of Eq. (1)

\tilde{R} (ρ_{0}, τ_{0}) = \frac{1}{Δ \tilde{A} Δ τ} \int_{τ_{0}}^{τ_{0} + Δ τ} \int_{ρ_{0}}^{ρ_{0} + Δ ρ} \int_{𝕊^{2 -}} \tilde{Λ} (ρ, \hat{Ω}, τ) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} 2 π ρ d ρ d τ

where

Δ \tilde{A} = π {(ρ_{0} + Δ ρ)}^{2} - π ρ_{0}^{2} .

Since the ratio between the radiance L (r,

\hat{Ω}

,t) and the dimensionless radiance

\tilde{Λ}

(ρ,

\hat{Ω}

,τ) is defined as L ₀, one can easily verify that the ratio between R(r,t) and R̃(ρ,τ) is also equal to L ₀. The triple integral in Eq. (1) represents a characteristic energy E ₀ that is emitted in the space-time bin considered. This characteristic energy E ₀ is equal to the product of the number photons N emitted in the bin multiplied by the energy of each photon, h ν. Because we integrate L(r,

\hat{Ω}

,t) and

\tilde{Λ}

(ρ,

\hat{Ω}

,τ) over equivalent ranges of solid angle, space, and time, we can write

\int_{t_{0}}^{t_{0} + Δ t} \int_{r_{0}}^{r_{0} + Δ r} \int_{𝕊^{2 -}} L (r, \hat{Ω}, t) d \hat{Ω} 2 π r d r d t = E_{0} \int_{τ_{0}}^{τ_{0} + Δ τ} \int_{ρ_{0}}^{ρ_{0} + Δ ρ} \int_{𝕊^{2 -}} \tilde{Λ} (ρ, \hat{Ω}, τ) d \hat{Ω} 2 π ρ d ρ d τ .

Using this equivalence, as well as the definitions for R(r ₀,t ₀) and R̃(ρ ₀,τ ₀) in Eq. (1) and Eq. (13), we find that L ₀ is the ratio between the reflectance and the dimensionless reflectance

L_{0} = \frac{R (r_{0}, t_{0})}{\tilde{R} (ρ_{0}, τ_{0})} = N h ν \frac{Δ \tilde{A} Δ τ}{Δ A Δ t} = E_{0} c μ_{s}^{3} .

With this definition for L ₀, the reflectance can be written in terms of the dimensionless reflectance as

R (r, t) = c μ_{s}^{3} E_{0} \tilde{R} (ρ, τ) = c μ_{s}^{3} E_{0} \int_{𝕊^{2 -}} \tilde{Λ} (ρ, \hat{Ω}, τ) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} .

This equation shows that the diffuse reflectance for a specific value of the scattering coefficient is uniquely related to the dimensionless reflectance multiplied by a reference radiance L ₀.

Thus Eq. (2) can be formulated by considering a reference reflectance R_r(r,t) for a non-absorbing medium with scattering coefficient μ_s = μ_s,r. Using Eq. (17), to express the reference reflectance within the KP sMC method in terms of the dimensionless reflectance, we get:

R_{r} (r, t) = c μ_{s, r}^{3} E_{0} {\tilde{R}}_{r} (ρ, τ) = c μ_{s, r}^{3} E_{0} {\tilde{R}}_{r} (μ_{s, r} r, c μ_{s, r} t),

while the diffuse reflectance for an arbitrary value of μ_s becomes

R (r, t) = c μ_{s}^{3} E_{0} \tilde{R} (ρ, τ) = c μ_{s}^{3} E_{0} \tilde{R} (\frac{μ_{s}}{μ_{s, r}} ρ_{r}, \frac{μ_{s}}{μ_{s, r}} τ_{r})

where ρ_r = μ_s,rr and τ_r = μ_s,r ct. By rewriting the dimensionless reflectance in the last equation as

\tilde{R} (\frac{μ_{s}}{μ_{s, r}} ρ_{r}, \frac{μ_{s}}{μ_{s, r}} τ_{r}) = R_{r} (\frac{μ_{s}}{μ_{s, r}} r, \frac{μ_{s}}{μ_{s, r}} t) \frac{1}{c μ_{s, r}^{3} E_{0}}

and substituting back into Eq. (19), we obtain Eq. (2)

R (r, t) = \frac{μ_{s}^{3}}{μ_{s, r}^{3}} R_{r} (\frac{μ_{s}}{μ_{s, r}} r, \frac{μ_{s}}{μ_{s, r}} t),

which demonstrates the consistency of the first sMC scaling relationship with the RTE.

The second scaling relation used within the sMC method, Eq. (3), relates the reflectance from a non-absorbing reference medium with the reflectance from a second medium possessing identical scattering properties but finite absorption. To evaluate the impact of absorption, consider that every photon emitted from the target surface is characterized by a photon exit position r_i, time of flight t_i, and weight w_i. While photon absorption is an analog phenomena, an effective technique to reduce the variance of the resulting photon weights is to diminish the photon weight continuously along its trajectory according to Beer’s law i.e., exp(−μ_act_i) [21]. This approach, known as continuous absorption weighting, provides a rigorous unbiased stochastic estimator of the radiative transport process [3]. Consider a reference simulation with zero absorption and calculate N_r(r ₀,t ₀) as the sum of the weights of the photons emitted from the top surface within a specific radial and temporal bin (r ₀,t ₀)

N_{r} (r_{0}, t_{0}) = Σ w_{i} = Σ w_{0}

where w ₀ = 1 – R_sp is the initial weight of the photon reduced by the specular reflection given by the refractive index mismatch between air and tissue. Thus, following Eq. (1), we can express the reference reflectance for the bin represented by r ₀ and t ₀ as

R_{r} (r_{0}, t_{0}) = \frac{1}{Δ A Δ t} \int_{t_{0}}^{t_{0} + Δ t} \int_{r_{0}}^{r_{0} + Δ r} \int_{𝕊^{2 -}} h ν c 𝒩_{r} (r_{0}, t_{0}, Ω) |_{z = 0^{-}} (\hat{n} \cdot \hat{Ω}) d \hat{Ω} 2 π r d r d t,

where 𝒩_r(r ₀,t ₀, Ω) [m⁻³ sr⁻¹] represents the volumetric photon density per unit solid angle at location r ₀ and time t ₀ in the reference MC simulation.

For a MC simulation with μ_a ≠ 0 the weight of each photon must be decreased exponentially according to Beer’s law, i.e., N_r (r ₀,t ₀) = Σw_i = Σw ₀ exp (−μ_act ₀). Thus the reflectance is

\begin{array}{l} R (r_{0}, t_{0}) & \approx & \frac{1}{Δ A Δ t} \int_{t_{0}}^{t_{0} + Δ t} \int_{r_{0}}^{r_{0} + Δ r} \int_{𝕊^{2 -}} h ν c 𝒩_{r} (r_{0}, t_{0}, Ω) |_{z = 0^{-}} exp (- μ_{a} c t_{0}) (\hat{n} \cdot \hat{Ω}) d \hat{Ω} 2 π r d r d t \\ = & R_{r} (r_{0}, t_{0}) exp (- μ_{a} c t_{0}) . \end{array}

If we consider continuous spatial and temporal variables we obtain the expression

R (r, t) = R_{r} (r, t) exp (- μ_{a} ct)

which recovers the second sMC scaling relationship Eq. (3). Equation (25) is rigorous only when applied individually to each emitted photon. That is, the value of t used must correspond to the photon’s time-of-flight rather than to a nominal time of the temporal bin in which it is tallied.

3. Materials and Methods

The flowcharts in Fig. 1 depict three approaches used here to obtain and compare MC estimates for spatially- and temporally-resolved reflectance from a homogeneous semi-infinite medium. In the first approach, a conventional Monte Carlo (cMC) simulation is used to provide “gold standard” reflectance estimates. The second method is based on the single Monte Carlo approach applied on a photon-by-photon basis (sMC_p). The third approach is similar to the KP method where the reflectance for any value of the optical scattering and absorption, R(r,t), is obtained by scaling a reference signal, R_r(r,t), obtained from a MC simulation with μ_s = μ_s,r and μ_a = 0. As this latter method employs interpolation of the reference reflectance prior to application of the scaling equations, we denote it as sMC_i. Below, we describe the detailed procedure used in each of these methods and also describe alternate binning and interpolation schemes that we have examined to potentially improve upon the KP sMC_i method.

Fig. 1 Flow chart representing the processes used to obtain the spatially and temporally resolved reflectance from the different Monte Carlo approaches examined in this study.

Download Full Size | PDF

3.1. Monte Carlo Simulations

We modified the well-known MCML MC simulation package developed by Jacques and Wang [21] to provide time-resolved output. All simulations run for this study inject photons in a semi-infinite homogeneous medium with n = 1.4 using a collimated ‘pencil’ beam oriented perpendicular to the top surface. The Henyey-Greenstein phase function with an anisotropy coefficient of g = 0.8 is used to model the angular distribution of single-scattering events. When a photon encounters the z = 0 surface, we select a pseudorandom number uniformly distributed in [0,1) and compare it to the Fresnel reflectance to determine whether the photon should exit or be internally reflected. To control the computation time, simulated photons are terminated before exiting the medium if they experience 3 × 10⁵ collisions. To achieve good statistics, the propagation of 10⁸ photons was simulated for each set of optical properties.

3.1.1. Conventional Monte Carlo (cMC)

In cMC simulations, the step size for photon propagation is sampled using the probability density function s = −ln(ξ)/μ_s, where μ_s is the scattering coefficient and ξ is a pseudorandom number uniformly distributed in the interval [0,1). For each photon simulated, the exit position r_i and total time of flight t_i are recorded. The weight of each photon w_i is then calculated using continuous absorption weighting according to Beer’s law for a specific value of μ_a. To reduce the computational burden we terminated all photons that traveled to a depth exceeding 1 m. After storing exit position, time, and weight (r_i,t_i, w_i) of each of the reflected photons in a database, the recorded photon weights were averaged over a set of spatial and temporal bins to obtain the reflectance signal [13].

3.1.2. Single Monte Carlo on a photon-by-photon basis (sMC_p)

The application of sMC on a photon-by-photon basis requires simulation of the photon trajectories and storage of the photon exit locations and times in a database. We perform a dimensionless MC simulation with zero absorption, set μ_s equal to a unitless value of one and sample the photon step size from s = −ln(ξ). The maximum depth for photon propagation was adjusted to match the termination condition of the cMC simulations. The dimensionless position and time (ρ_i,τ_i) of each exiting photon was recorded. Since μ_a = 0 in this simulation, all the emitted photons had the same weight w_i = 1 – R_sp. To obtain the temporally- and spatially-resolved reflectance for any value of μ_a and μ_s we first apply Eq. (7) to determine the exit position and time, r_i and t_i, and apply the scaling relationships, Eqs. (2) and (3) to each exiting photon to determine the proper photon weight. Once the exiting position, time, and weight of each photon is obtained, the individual photon weights are averaged over a set of spatial and temporal bins to obtain the reflectance signal as in the cMC simulations.

3.1.3. Interpolated single Monte Carlo (sMC_i)

sMC_i refers to the method established by Kienle and Patterson which uses a cMC simulation to form the reference reflectance R_r(r,t) using reference optical properties μ′_s,r = 1 mm⁻¹ and μ_a = 0. The resulting photon database is then processed to obtain the reference reflectance for a set of nominal positions, R_r(r_k,t_l). The reference values sampled through the MC simulation are smoothed by fitting to an analytical function and then interpolated. The reflectance for any combination of optical properties is obtained by applying the scaling relations Eqs. (2) and (3) to the interpolated representation of the reflectance.

3.2. Error Analysis

Our objective is to compare both sMC_p and sMC_i relative to the conventional MC simulations. The relative error for the spatially- and temporally- resolved reflectance is given by:

ɛ_{rel (p, i)} (r_{k}, t_{l}) = \frac{| R_{cMC} (r_{k}, t_{l}) - R_{sMC (p, i)} (r_{k}, t_{l}) |}{R_{cMC} (r_{k}, t_{l})} .

The cMC values chosen as gold standard are stochastic estimates of the RTE solution. The relative error estimates are meaningful only if the relative uncertainty associated with the R _cMC(r_k,t_l) estimates is at least an order of magnitude smaller. This relative uncertainty is calculated as

σ_{rel} (r_{k}, t_{l}) = \frac{σ_{cMC} (r_{k}, t_{l})}{R_{cMC} (r_{k}, t_{l})} .

3.3. Binning

To represent the MC data, the detected photon weights must be tallied in a set of radial and temporal bins of finite width. For each reflected photon emitted at radial position r_i and time t_i, a tally of w_i is added to the reflectance R(r_k,t_l) where k is the index associated with the radial ring that contains r_i and l is the index for the time interval that includes t_i. This process is necessary to obtain R(r_k,t_l) using either the cMC or sMC_p approaches. For the KP sMC_i approach, the reference MC simulation is first binned to obtain the nominal values R(r_k,t_l). The binned reference MC simulation is then processed to obtain R_r(r,t) as illustrated in Fig. 1.

Three different discretization techniques have been examined in this work. The first is commonly known as equal width discretization (EWD) [25] and is based on the subdivision of the physical range into m_r equispaced radial bins and m_t equispaced time bins resulting in a total of m_r × m_t bins. This binning strategy is most commonly adopted for the evaluation of the reflectance signal from the output of a cMC simulation [17]. A slight variation presented in the literature utilizes multiple sets of uniformly-sized bins [11].

The second approach, known as equal frequency discretization (EFD) [26], divides the sorted values into m_r × m_t intervals of irregular width so that each interval contains approximately the same number of samples. Instead of fixing the boundary of the bins a priori, the bins are built dynamically based on the radial and temporal cumulative distribution functions of the collected photons. Thus each interval contains N_out/(m_r × m_t) photons where N_out is the number of photons that exited the phase space and were not terminated during the simulation. This binning process requires sorting of the photon database according to the radial position. Once the data set is split into m_r bins in space with the same number of photons in each bin, the photons in each spatial bin are sorted according to the photon time of flight. The reflectance values are then computed by using EFD along the time dimension. Both EFD and EWD were used to provide a discrete representation of the reflectance R(r_k,t_l) from the cMC simulation output. For all the combinations of optical properties under investigation we save the bin limits obtained using EFD and use them to obtain R(r_k,t_l) from the photon database generated with sMC_p.

We developed a third binning strategy based on an adaptive width discretization (AWD) algorithm to construct the reference reflectance used for the sMC_i. In AWD the number of bins is not imposed a priori and the bin limits are dynamically built to ensure that each bin has a sufficient sample size of photons to reduce variance while also spacing the bins so that they accurately capture any rapid spatial and/or temporal dynamics of the reflectance signal. To accomplish this, we established the following binning criteria: (i) a minimum of n_r = 5 × 10⁵ photons must occupy each radial bin; (ii) the number of photons occupying each time bin, n_t (r), increases from 2500 to 5000 over the radial distance range; (iii) the interval between two adjacent radial bins cannot exceed 0.2 mm i.e., (r_k ₊₁ – r_k) ≤ Δr_max = 0.2 mm; (iv) the interval between two adjacent time bins cannot exceed 20 ps, i.e., (t_l ₊₁ – r_l) ≤ Δt_max = 20 ps; and (v) the relative change in the reflectance value between two adjacent time bins cannot exceed 10% i.e., (R_k ₊₁ – R_k)/R_k ≤ ΔR_rel,max = 0.1. While other values for n_r, n_t (r), Δr_max, Δt_max and ΔR_min were tested, the values above provided a smaller variance for the reflectance values R(r_k,t_l) and reduced noise in the data while accurately capturing sharp features of the reflectance signal.

Independent from any specific binning algorithm, the representation of all the photons collected within a bin with nominal coordinates (r_k,t_l) introduces a discretization error. Typically each bin is represented by its midpoint. To reduce this error we used the mean of the radial position and time of flight of all the photons collected within each bin. While this error reduction strategy has not been examined carefully, we noted that when large bins are present in EFD or EWD, the use of the midpoint leads to artifacts in the measured reflectance. Regardless of the binning approach used, the value of the reflectance calculated over an annular bin of area ΔA_k and time span Δt_l is computed as [21]

R (r_{k}, t_{l}) = \frac{1}{Δ A_{k}} \frac{1}{Δ t_{l}} \frac{1}{N} \sum_{t = 1}^{N} w_{i},

where N is the total number of injected photons. The weight w_i is set to zero for all the photons collected outside the bin under consideration. The standard deviation of the reflectance measured in the specific bin is calculated as [27]:

σ (r_{k}, t_{l}) = \frac{1}{\sqrt{N - 1}} {[\sum_{i = 1}^{N} \frac{w_{i}^{2}}{Δ A_{k}^{2} Δ t_{l}^{2} N} - R^{2} (r_{k}, t_{l})]}^{1 / 2} .

3.4. Reference Reflectance

In principle, the use of Eqs. (2) and (3) enables a continuous representation of R(r,t) for any absorption and scattering coefficient. However when utilizing the results of a cMC simulation with a finite number of photons, it is possible to measure the reflectance only for a discrete set of positions. Thus, it is necessary to interpolate between these points to provide a continuous representation of R(r,t) for an arbitrary set of optical properties. To represent the reflectance signal over the physical domain of interest, we establish a native set of reflectance values, R_r(r_k,t_l) over a time range from 0 to 20 ns and a spatial range from 0 to 100 mm.

The adaptive binning technique described above is used to provide a representation of the reference simulation output which must subsequently be interpolated. Formation of an interpolating surface requires a support mesh of points. Since the number of the reflectance values and their temporal locations are different for every radial bin, we resample each reflectance curve R(r_k = constant,t_l) over a uniform time vector T to obtain a regular grid that is more suitable for interpolation. These time points are chosen in a way that reflects the cumulative distribution function of all the photons, using a decreasing density of sampled points per unit time in the range [0–20] ns.

3.5. Interpolation and Resampling

The AWD approach had one limitation in that small numbers of photons were typically collected at late times. For this reason the last time coordinate value that satisfied condition (ii), R(r_k,t_last), occurs at a time smaller than 20 ns for radial distances r_k < 15 mm. To obtain nominal values for the reflectance in this time range we defined a set of bins 100 ps in width from this valid time location to 20 ns. For each radial bin that exhibited this problem, the photons collected within its limits are tallied in these bins. The variance associated with the photon collected in these 100-ps-width bins is still too large to use them directly for the creation of the reference surface. Thus we fit these points to the following nonlinear function, f_t (t), with three free parameters (a,b,c)

f_{t} (t) = \frac{a}{t^{3 / 2}} [exp (- \frac{b}{t}) - exp (- \frac{c}{t})] .

This function has the form of a source-image Green’s function without absorption decay within the standard diffusion approximation [2]. This functional form fits the data very well most likely because the photons collected at these remote distances and long times have undergone a very large number of collisions. Figure 2(a) shows the result of this approach to represent the reference reflectance for r_k = 1 mm.

Fig. 2 (a) Results of the fitting approach used to extrapolate points on the tail of the time-resolved reflectance R(t) for r_k = 1 mm. (b) Values of R_r(r_k,t_l) for r_k ∈ [5 – 100] mm obtained after binning, extrapolating and resampling.

Download Full Size | PDF

To improve the accuracy in the representation of the initial rise of the reference reflectance signal we introduce additional time points at earlier times via extrapolation. For every radial bin we record the arrival time of the first photon t ₀. The time span from t ₀ to the nominal position of the first bin, t ₁, was split into ten uniform bins in which the photons were tallied to provide additional reflectance values. The logarithm of the reflectance values sampled with this method is fitted to a third degree polynomial. For each radial bin, this polynomial is used to extrapolate values for the points of T between the arrival time of the first photon, t ₀, and the time of the first reflectance value obtained with the binning process, t ₁.

To resample the reflectance curves at the temporal locations specified by T, we fit a 3rd degree Non Uniform Rational B-spline (NURBS) [28,29] curve to each R(r_k = const,t_l) curve. After calculating the approximating NURBS curve we resample it to obtain R_r(r_k = const,t_l = T). The values of R_r(r_k,t_l) obtained after resampling is the final set used in the sMC_i method and are shown in Fig. 2(b) for radial distances spanning 5 to 100 mm. The noise reduction obtained by resampling allowed the use of the reference points R_r(r_k,t_l) as native points for the interpolation necessary to calculate R(r,t) for continuous values of r and t for an arbitrary set of optical properties using Eqs. (2) and (3). We use a 3rd degree NURBS interpolation in both position and time to represent R_r(r_k,t_l) as a 2D surface.

4. Results

4.1. cMC Simulations and Effects of Discretization

We first consider the relative standard deviation of the reflectance provided by our “gold standard” cMC simulations as these are the results against which the sMC approaches will be compared. We calculate the relative standard deviation using Eq. (27) after binning the reflected photons provided by the cMC simulation using either EFD or EWD. We used a set of 200×200 bins to calculate the reflectance values over spatial and temporal ranges of r ∈ [0, 50] mm and t ∈ [0, 5] ns, respectively. We performed cMC simulations for 15 sets of optical properties where the absorption and reduced scattering coefficients were in the range μ_a ∈ [10⁻³, 0.3]mm⁻¹ and μ′_s ∈ [0.5, 2]mm⁻¹

Figure 3(a) provides a colored contour plot of spatial and temporal dependence of the relative standard deviation of the reflectance estimates σ_rel(r,t) for the case of μ′_s = 1mm⁻¹ and μ_a = 0.1mm⁻¹. Figure 3 shows that the values of σ_rel are remarkably uniform over the entire domain. Across all 15 sets of optical properties examined, the mean of the relative standard deviation across all space and time using EFD was less than 3%. For larger values of the absorption coefficient, σ_rel increases for combinations of short distances and late times and for long distances and early times. This is because the reflectance signal is weaker in these regions and the span covered by the EFD bins is large at late times in order to collect the required number of photons. The relative standard deviation at any given position/time pair exceeds 5% only for μ_a ≥ 0.1 mm⁻¹. By contrast, in Fig. 3(b) we see that the relative standard deviation obtained with EWD is not uniform and increases strongly for longer times and larger distances with a mean value larger than 5%. Moreover, for all the sets of optical properties examined, there were position/time pairs for which the relative standard deviation in the reflectance estimates exceeded 100%.

Fig. 3 Relative standard deviation of the reflectance provided by discretizing the results of a cMC simulation for μ_a = 0.1mm⁻¹ and μ′_s = 1mm⁻¹ using (a) 200 × 200 EFD bins, (b) 200 × 200 EWD bins, (c) 250 × 1000 EWD bins.

Download Full Size | PDF

A comparison of the uncertainty obtained using the same number of bins using the two techniques may lead one to underappreciate the improvement obtained with EFD. It is true that in certain regions, the uncertainty in the cMC estimates obtained when using EWD with a 200 × 200 set of uniform bins is smaller than that obtained using EFD. This is because the uniform bins have radial and temporal width of Δr = 0.25 mm and Δt = 25 ps. Thus for combinations of small radial positions and times, the number of photons tallied within these limits is very large and the statistical noise is reduced. However bins that are 25 ps in width are too large to capture the initial peak of the reflectance signal. For a more fair comparison we use EWD with 250 radial intervals of width Δr = 0.2 mm and 1000 time bins of width Δt = 5 ps. These results are given in Fig. 3(c) and shows that the performance of the EWD at large times and radial positions is adversely effected when the number of bins is increased. In fact for this set of 250×1000 bins the mean σ_rel across all the time and spatial bins exceeds 15% for all the sets of optical properties examined. Thus while EWD is simple to implement, the relative standard deviation will exhibit significant variability since the number of photons tallied in each of the bins varies widely. EFD, on the other hand, reduces the variation in the number of tallied photons across all bins by imposing finer sampling in the regions of phase space where number density of reflected photons is large while also imposing coarser discretization in regions where the reflected signal is sparse. This has the overall effect of not only reducing the relative standard deviation but also also minimizing its variation over the phase space.

4.2. Comparison of single Monte Carlo Techniques

Now with the relative standard deviation of the “gold-standard” cMC simulations characterized, we can determine the relative error provided by sMC_p and sMC_i approaches established in §3 as compared to the reflectance provided by a cMC simulation R _cMC(r_k,t_l) represented by a set of 200 × 200 radial/temporal bins formed using EFD.

For all the combinations of optical properties under investigation, we process the sMC_p data set to evaluate the values (r_i,t_i, w_i) from (ρ_i,τ_i) using the scaling imposed by the value of μ_s and modifying the photon weight in media with μ_a ≠ 0 by applying Beer’s law using the values of the photon exit time t_i. To calculate R _{sMC_p} (r_k,t_l) for the same nominal spatial and temporal locations, we use the bins defined by the cMC simulation for the specific combination of the optical properties to tally the photon weights and build up the reflectance values. We use the NURBS surface obtained by the procedure described in §3.4 to calculate R _{sMC_i} (r_k,t_l). First, the scaling given by the scattering coefficient is imposed using Eq. (2). Once established, the effects of non-zero absorption is calculated using Eq. (3).

Figures 4, 5, and 6 show (a) the reflectance determined from the cMC simulations on a logarithmic scale along with the relative error in the reflectance predictions provided by the (b) sMC_p and (c) sMC_i approaches for different combinations of optical properties. The spatial and temporal ranges shown extend only to 20 mm and 1 ns to provide reasonable resolution over the region where the largest contribution to the reflectance signal resides. Figure 4(b) shows that the sMC_p method has a relative error of ≲ 3% over the entire spatial and temporal domain considered. As this is roughly equal to the maximum relative standard deviation of the 200 × 200 EFD representation of the cMC, we can conclude that the reflectance estimates obtained using the cMC and the sMC_p methods are equivalent. In fact, for all sets of optical properties investigated, the mean value of the relative error ɛ _{rel,sMC_p} (r_k, t_l) over the entire physical range (r ∈ [0, 50] mm and t ∈ [0, 5] ns) never exceeds 3% which is smaller than the mean of the relative standard deviation associated with the cMC measurements. Thus there appears to be no statistically significant difference between the reflectance estimates given by cMC and sMC_p simulations when using a set of 200 × 200 radial/temporal bins formed using EFD.

Fig. 4 (a) Reflectance signal obtained with EFD for a cMC simulation, (b) relative error of the reflectance provided by the sMC_p method, (c) relative error of the reflectance provided by the sMC_i, for μ_a = 0.001 mm⁻¹ and μ′_s = 2 mm⁻¹.

Download Full Size | PDF

Fig. 5 (a) Reflectance signal obtained with EFD from the cMC simulation, (b) relative error of the reflectance provided by the sMC_p method, (c) relative error of the reflectance provided by the sMC_i method, for μ_a = 0.01 mm⁻¹ and μ′_s = 1 mm⁻¹.

Download Full Size | PDF

Fig. 6 (a) Reflectance signal obtained with EFD from the cMC simulation, (b) relative error of the reflectance provided by the sMC_p method, (c) relative error of the reflectance provided by the sMC_i method, for μ_a = 0.1 mm⁻¹ and μ′_s = 1.5 mm⁻¹.

Download Full Size | PDF

The relative error of the reflectance calculated using the sMC_i approach displays two characteristic features regardless of optical properties. Firstly, the relative error in the first temporal bin can be as large as 100%. It is important to note that this large error is measured only for the first time bin, a location where the reflectance value is much smaller then the peak value. Secondly, the relative standard deviation associated with the first time bin at large source detector separations is generally larger than the average error measured over the entire physical range.

We also see specific features of the relative error depending on optical properties. For example, for μ_a values smaller than 0.1 mm⁻¹ across all μ′_s values, the mean relative error of the sMC_i approach tallied over the entire temporal and radial range excluding the first temporal bin is smaller than that obtained with sMC_p method. This reflects the advantage gained by interpolating the reference signal which reduces the impact of statistical noise. For μ_a values equal to 0.1 mm⁻¹ and greater, we see that while the mean relative error for the reflectance values provided by the sMC_i approach is larger than that provided by the sMC_p approach, the error itself remains below 5%. From Figs. 4(c) and 5(c) one can see that the error introduced by using Eq. (30) to extrapolate the reflectance at long times is very small for small μ_a values. The values obtained in the latest time bins for short distances are obtained by scaling the reference reflectance in a range where the nominal reference values were obtained with the fitting approach, and the relative error remains very small.

For higher μ_a values, the error measured in the latest time bins at small radial locations increases as shown in Fig. 6(c). the error exceeds 50% for μ_a ≥ 0.1 mm⁻¹ and approaches 100% for μ_a = 0.3 mm⁻¹ over the physical range shown. This large error is due to the fact that Eq. (3) is approximate when applied to binned data where a median value of t is used to rescale the weights of all photons captured in a discrete time bin. It is important to note that the relative standard deviation of these bins is also large. While these relative error values are large, one must recognize that the reflectance values span a significantly larger dynamic range for μ_a values of 0.1 mm⁻¹ and larger. Specifically, within the radial and temporal range considered, the reflectance values span at least 8 orders of magnitude. When viewed in this context, we see that these large relative errors occur in regions where the reflectance is already very small.

Since the reference signal represented by the NURBS surface extends to temporal and spatial values of 20 ns and 100 mm, respectively, it is possible to overcome the dynamic range limitation previously reported. The accuracy of the results given by the sMC_i method for large distances and long times decreases with increasing absorption. For μ_a = 0.001 mm⁻¹ if we exclude the first temporal bin, the relative error does not exceed 10% across all the values of the reduced scattering coefficient examined. For increasing μ_a values, the relative error increases for large distances and late times. Larger μ′_s values increase the temporal range over which the error remains low, but decrease the spatial range where the model is accurate. This may be due to the fact that a larger scattering coefficient increases the number of collisions that a photon undergoes before exiting the tissue thus it increases the signal at later times, while it concentrates the light distribution over a smaller radial range.

As reported by Kienle and Patterson, interpolation errors are largest for short times and small distances and well as for μ′_s values far removed from the reference μ′_s value. This arises due to shortcomings of discretization methods that tend to be inaccurate in regions where rapid spatial/temporal variations of the reflectance signal exist. Figure 7 compares the the relative error measured for times ≤0.1 ns and source detector locations ≤15 mm for μ_a = 0.1 mm⁻¹ and μ′_s = 0.5 mm⁻¹ using (a) the NURBS surface built upon the reference nominal data or (b) linear interpolation of the reflectance binned in uniform bins with Δt = 5 ps and Δr = 0.2 mm. This result shows that the NURBS sMC_i approach typically reduces the relative error by an order of magnitude and provides accuracy comparable to that obtained using the sMC_p approach. This provides a reflectance signal with excellent accuracy even when the optical properties chosen represent a large perturbation from the reference reflectance.

Fig. 7 (a) Log₁₀ of the relative error of the sMC_i approach using NURBS interpolation on the reference values, (b) log₁₀ of the relative error of the sMC_i approach based on linear interpolation of the uniformly binned reflectance as done by Kienle and Patterson for μ_a = 0.1 mm⁻¹ and μ′_s = 0.5 mm⁻¹

Download Full Size | PDF

5. Conclusion

We have demonstrated how the scaling relations used for the single Monte Carlo (sMC) method can be derived from the radiative transport equation. While the scaling associated with the scattering coefficient is rigorous, the application of Beer’s law to model the effect of absorption is only approximate when applied to binned data. As a result, the accuracy of the sMC approach degrades for larger μ_a values and/or when the discretization intervals are large. We show that when sMC is applied on a photon by photon basis it yields reflectance estimates that are statistically equivalent to those obtained using conventional MC simulations.

However a more practical implementation of the sMC requires the development of improved spatial and temporal binning and interpolation approaches. The challenge in such methods is to construct a reference reflectance in a fashion that minimizes the discretization error and error propagation due to interpolation. We demonstrate the advantages of the equal frequency discretization approach as compared to equal width discretization, which provides a low, and nearly uniform, relative error over the entire measurement domain. Even better performance can be realized using a novel adaptive binning strategy that accommodates the large dynamic range of the reflectance values while also reducing interpolation errors. The regularization of the number of photons collected in each bin leads to an almost uniform relative standard deviation across all the measured reflectance values over a wide physical range. Moreover the windowing mechanism implemented in the adaptive strategy provides a more dense network of nominal values that can be robustly interpolated.

Our implementation of NURBS interpolation provides a fast and stable algorithm that requires a small amount of information as compared to other fitting methods. The results shown here indicate that the sMC_i method can be used to derive R(r,t) for a wide range of optical properties and source detector separations. Given its speed, the sMC_i method can potentially be used as an inverse solver engine for the recovery of optical properties from experimental data as it provides transport-rigorous estimates, is based on very fast and stable algorithms and requires little data storage. The sMC_i forward solver based on adaptive binning and NURBS interpolation (NURBS Forward Solver) is available for use in the Virtual Tissue Simulator at http://www.virtualphotonics.org/.

Acknowledgments

We acknowledge support from the Laser Microbeam and Medical Program a NIH Biomedical Technology Resource ( NIH-P41-RR01192). We thank the UC Irvine Office of Information and Technology for use of the Broadcom Distributed United Cluster.

References and links

1. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]

2. A. H. Hielscher, “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol. 40, 1957–1975 (1995). [CrossRef] [PubMed]

3. J. Spanier and E. M. Gelbard, Monte Carlo Principles and Neutron Transport Problems (Dover Publications, 2008).

4. A. D. Kim, “Transport theory for light propagation in biological tissue,” J. Opt. Soc. Am. A 21, 820–827 (2004). [CrossRef]

5. A. H. Hielscher, R. E. Alcouffe, and R. E. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43, 1285–1302 (1998). [CrossRef] [PubMed]

6. R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

7. A. Kienle and M. S. Patterson, “Determination of the optical properties of semi-infinite turbid media from frequency-domain reflectance close to the source,” Phys. Med. Biol. 42, 1801–1819 (1997). [CrossRef] [PubMed]

8. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28, 2331–2336 (1989). [CrossRef] [PubMed]

9. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and distributions of light in tissue,” Med. Phys. 10, 824–830 (1983). [CrossRef] [PubMed]

10. R. Graaff, M. H. Koelink, F. F. M. De Mul, W. G. Zijistra, A. C. M. Dassel, and J. G. Aarnoudse, “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt. 32, 426–434 (1993). [CrossRef] [PubMed]

11. A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol. 41, 2221–2227 (1996). [CrossRef] [PubMed]

12. A. Pifferi, P. Taroni, G. Valentini, and S. Andersson-Engels, “Real-time method for fitting time-resolved reflectance and transmittance measurements with a Monte Carlo model,” Appl. Opt. 37, 2774–2780 (1998). [CrossRef]

13. E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. 13, 041301 (2008). [CrossRef]

14. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45, 1062–1071 (2006). [CrossRef] [PubMed]

15. G. M. Palmer, C. F. Zhu, T. M. Breslin, F. S. Xu, K. W. Gilchrist, and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Application to breast cancer,” Appl. Opt. 45, 1072–1078 (2006). [CrossRef] [PubMed]

16. H. Xu, T. J. Farrell, and M. S. Patterson, “Investigation of light propagation models to determine the optical properties of tissue from interstitial frequency domain fluence measurements,” J. Biomed. Opt. 11, 1–18 (2006). [CrossRef]

17. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. R. Ayers, and B. J. Tromberg, “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt. 14, 024012 (2009). [CrossRef] [PubMed]

18. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26, 1335–1337 (2001). [CrossRef]

19. I. Seo, J. S. You, C. K. Hayakawa, and V. Venugopalan, “Perturbation and differential Monte Carlo methods for measurement of optical properties in a layered epithelial tissue model,” J. Biomed. Opt. 12, 1–15 (2007). [CrossRef]

20. A. Kienle, M. S. Patterson, N. Dognitz, R. Bays, G. Wagnières, and H. van Den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

21. S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in multi-layered tissue,” Comput. Methods Programs Biomed. 47, 131–146 (1995). [CrossRef] [PubMed]

22. W. Cheong and S. A. Prahl. “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26, 2166–2185 (1990). [CrossRef]

23. L. G. Henyey and J. L. Greenstein, “Diffuse radiation of the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]

24. W. Hines, C. D. Montgomery, M. D. Goldsman, and M. C. Borror, Probability and Statistics in Engineering (John Wiley & Sons Inc., 2003).

25. J. Dougherty, R. Kohavi, and S. Mehran, “Supervised and unsupervised discretization of continuous features,” in Proceedings 12th International Conference on Machine Learning (Morgan Kaufmann, 1995).

26. M. Boulle, “Optimal bin number for equal frequency discretization,” Intell. Data Anal. 9, 175–188 (2005).

27. M. H. Kalos and P. A. Whitlock, Monte Carlo Methods Volume 1: Basics (Wiley-Interscience, 1995).

28. W. Tiller and L. Piegl, The NURBS Book (Springer, 1995).

29. F. D. Rogers, An Introduction to NURBS With Historical Perspective (Morgan Kaufmann, 2004).

Analysis of single Monte Carlo methods for prediction of reflectance from turbid media

Abstract

1. Introduction

2. Theory