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

Optics Express, Vol. 19, Issue 20, pp. 19627-19642 (2011)

http://dx.doi.org/10.1364/OE.19.019627

Acrobat PDF (1615 KB)

### Abstract

Starting from the radiative transport equation we derive the scaling relationships that enable a single Monte Carlo (MC) simulation to predict the spatially- and temporally-resolved reflectance from homogeneous semi-infinite media with arbitrary scattering and absorption coefficients. This derivation shows that a rigorous application of this single Monte Carlo (sMC) approach requires the rescaling to be done individually for each photon biography. We examine the accuracy of the sMC method when processing simulations on an individual photon basis and also demonstrate the use of adaptive binning and interpolation using non-uniform rational B-splines (NURBS) to achieve order of magnitude reductions in the relative error as compared to the use of uniform binning and linear interpolation. This improved implementation for sMC simulation serves as a fast and accurate solver to address both forward and inverse problems and is available for use at

© 2011 OSA

## 1. Introduction

*μ′*=

_{s}*μ*(1 –

_{s}*g*), is much larger than the absorption coefficient

*μ*and over spatial scales that are larger than the transport mean free path,

_{a}*l*= (

^{*}*μ*+

_{a}*μ′*)

_{s}^{−1}[1

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]

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]

*a*=

*μ*/(

_{s}*μ*+

_{a}*μ*) 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,

_{s}*μ*=

_{t}*μ*+

_{a}*μ*, was held constant.

_{s}*R*(

*r,t*) for a homogeneous semi-infinite medium with arbitrary optical properties [11

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]

*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

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

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]

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]

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]

## 2. Theory

*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*

_{0},

*t*

_{0}) on the tissue surface can be estimated by where

*L*(

*r*,

*t*) is the radiance,

*n̂*the outward pointing unit normal vector,

^{2−}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*

_{0},

*t*

_{0}), a discretization error occurs because the quadrature interval has been represented by the nominal values (

*r*

_{0},

*t*

_{0}). This error vanishes as the number of simulated photons increases without bound and as the integration intervals → 0.

*R*(

_{r}*r,t*) with

*μ*= 0 and

_{a}*μ*=

_{s}*μ*to determine the reflectance

_{s,r}*R*(

*r,t*) from a second non-absorbing medium with arbitrary scattering coefficient

*μ*[11

_{s}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]

*μ*= 0) to a second medium with identical scattering but non-zero absorption

_{a,r}*μ*: where

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

*L*(

**r**,

*t*) is the radiance at a given location

**r**at at given time

*t*in a given direction

*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*

_{0}/

*n*).

*μ*and

_{s}*μ*are the scattering and absorption coefficients, respectively, and

_{a}*p*(

*is the angular distribution of the Fresnel reflectance for unpolarized light and*

_{F}*is the unit vector in the direction of the reflected ray for an incident direction*

_{r}*. Thus ℛ*

_{i}*(*

_{F}*n̂*·

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

_{i}*τ*is a dimensionless time,

*ρ*a dimensionless position and

*ρ*,

*τ*) a dimensionless radiance defined as the ratio between the radiance,

*L*(

**r**,

*t*), and a characteristic radiance

*L*

_{0}. Substitution of these variables into the RTE yields: The two boundary conditions, Eqs. (5) and (6) become 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 By analogy, a dimensionless reflectance is defined as

*R*(

_{r}*r,t*) for a non-absorbing medium with scattering coefficient

*μ*=

_{s}*μ*. Using Eq. (17), to express the reference reflectance within the KP sMC method in terms of the dimensionless reflectance, we get: while the diffuse reflectance for an arbitrary value of

_{s,r}*μ*becomes where

_{s}*ρ*=

_{r}*μ*and

_{s,r}r*τ*=

_{r}*μ*

_{s,r}*ct*. By rewriting the dimensionless reflectance in the last equation as and substituting back into Eq. (19), we obtain Eq. (2) which demonstrates the consistency of the first sMC scaling relationship with the RTE.

*r*, time of flight

_{i}*t*, and weight

_{i}*w*. 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(−

_{i}*μ*) [21

_{a}ct_{i}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]

*N*(

_{r}*r*

_{0},

*t*

_{0}) as the sum of the weights of the photons emitted from the top surface within a specific radial and temporal bin (

*r*

_{0},

*t*

_{0}) where

*w*

_{0}= 1 –

*R*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

_{sp}*r*

_{0}and

*t*

_{0}as where 𝒩

*(*

_{r}*r*

_{0},

*t*

_{0}, Ω) [m

^{−3}sr

^{−1}] represents the volumetric photon density per unit solid angle at location

*r*

_{0}and time

*t*

_{0}in the reference MC simulation.

*μ*≠ 0 the weight of each photon must be decreased exponentially according to Beer’s law, i.e.,

_{a}*N*(

_{r}*r*

_{0},

*t*

_{0}) = Σ

*w*= Σ

_{i}*w*

_{0}exp (−

*μ*

_{a}ct_{0}). Thus the reflectance is

*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 third approach is similar to the KP method where the reflectance for any value of the optical scattering and absorption,*

_{p}*R*(

*r,t*), is obtained by scaling a reference signal,

*R*(

_{r}*r,t*), obtained from a MC simulation with

*μ*=

_{s}*μ*and

_{s,r}*μ*= 0. As this latter method employs interpolation of the reference reflectance prior to application of the scaling equations, we denote it as sMC

_{a}*. 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.*

_{i}### 3.1. Monte Carlo Simulations

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]

*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

^{5}collisions. To achieve good statistics, the propagation of 10

^{8}photons was simulated for each set of optical properties.

#### 3.1.1. Conventional Monte Carlo (cMC)

*s*= −ln(

*ξ*)/

*μ*, where

_{s}*μ*is the scattering coefficient and

_{s}*ξ*is a pseudorandom number uniformly distributed in the interval [0,1). For each photon simulated, the exit position

*r*and total time of flight

_{i}*t*are recorded. The weight of each photon

_{i}*w*is then calculated using continuous absorption weighting according to Beer’s law for a specific value of

_{i}*μ*. To reduce the computational burden we terminated all photons that traveled to a depth exceeding 1 m. After storing exit position, time, and weight (

_{a}*r*,

_{i}*t*,

_{i}*w*) 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

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

#### 3.1.2. Single Monte Carlo on a photon-by-photon basis (sMC_{p})

_{p}

*μ*equal to a unitless value of one and sample the photon step size from

_{s}*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}*τ*) of each exiting photon was recorded. Since

_{i}*μ*= 0 in this simulation, all the emitted photons had the same weight

_{a}*w*= 1 –

_{i}*R*. To obtain the temporally- and spatially-resolved reflectance for any value of

_{sp}*μ*and

_{a}*μ*we first apply Eq. (7) to determine the exit position and time,

_{s}*r*and

_{i}*t*, 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.

_{i}#### 3.1.3. Interpolated single Monte Carlo (sMC_{i})

_{i}

*refers to the method established by Kienle and Patterson which uses a cMC simulation to form the reference reflectance*

_{i}*R*(

_{r}*r,t*) using reference optical properties

*μ*′

*= 1 mm*

_{s,r}^{−1}and

*μ*= 0. The resulting photon database is then processed to obtain the reference reflectance for a set of nominal positions,

_{a}*R*(

_{r}*r*,

_{k}*t*). 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.

_{l}### 3.2. Error Analysis

*and sMC*

_{p}*relative to the conventional MC simulations. The relative error for the spatially- and temporally- resolved reflectance is given by: 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*

_{i}*R*

_{cMC}(

*r*,

_{k}*t*) estimates is at least an order of magnitude smaller. This relative uncertainty is calculated as

_{l}### 3.3. Binning

*r*and time

_{i}*t*, a tally of

_{i}*w*is added to the reflectance

_{i}*R*(

*r*,

_{k}*t*) where

_{l}*k*is the index associated with the radial ring that contains

*r*and

_{i}*l*is the index for the time interval that includes

*t*. This process is necessary to obtain

_{i}*R*(

*r*,

_{k}*t*) using either the cMC or sMC

_{l}*approaches. For the KP sMC*

_{p}*approach, the reference MC simulation is first binned to obtain the nominal values*

_{i}*R*(

*r*,

_{k}*t*). The binned reference MC simulation is then processed to obtain

_{l}*R*(

_{r}*r,t*) as illustrated in Fig. 1.

*m*equispaced radial bins and

_{r}*m*equispaced time bins resulting in a total of

_{t}*m*×

_{r}*m*bins. This binning strategy is most commonly adopted for the evaluation of the reflectance signal from the output of a cMC simulation [17

_{t}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]

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]

*m*×

_{r}*m*intervals of irregular width so that each interval contains approximately the same number of samples. Instead of fixing the boundary of the bins

_{t}*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*) photons where

_{t}*N*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

_{out}*m*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*(

*r*,

_{k}*t*) 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

_{l}*R*(

*r*,

_{k}*t*) from the photon database generated with sMC

_{l}*.*

_{p}*. In AWD the number of bins is not imposed*

_{i}*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*= 5 × 10

_{r}^{5}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}_{+1}–

*r*) ≤ Δ

_{k}*r*= 0.2 mm; (

_{max}*iv*) the interval between two adjacent time bins cannot exceed 20 ps, i.e., (

*t*

_{l}_{+1}–

*r*) ≤ Δ

_{l}*t*= 20 ps; and (

_{max}*v*) the relative change in the reflectance value between two adjacent time bins cannot exceed 10% i.e., (

*R*

_{k}_{+1}–

*R*)/

_{k}*R*≤ Δ

_{k}*R*= 0.1. While other values for

_{rel,max}*n*,

_{r}*n*(

_{t}*r*), Δ

*r*, Δ

_{max}*t*and Δ

_{max}*R*were tested, the values above provided a smaller variance for the reflectance values

_{min}*R*(

*r*,

_{k}*t*) and reduced noise in the data while accurately capturing sharp features of the reflectance signal.

_{l}*r*,

_{k}*t*) 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 Δ

_{l}*A*and time span Δ

_{k}*t*is computed as [21

_{l}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]

*N*is the total number of injected photons. The weight

*w*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]:

_{i}### 3.4. Reference Reflectance

*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*) over a time range from 0 to 20 ns and a spatial range from 0 to 100 mm.

_{l}*R*(

*r*= constant,

_{k}*t*) over a uniform time vector

_{l}*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

*ii*),

*R*(

*r*,

_{k}*t*), occurs at a time smaller than 20 ns for radial distances

_{last}*r*< 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,

_{k}*f*(

_{t}*t*), with three free parameters (

*a,b,c*) This function has the form of a source-image Green’s function without absorption decay within the standard diffusion approximation [2

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]

*r*= 1 mm.

_{k}*t*

_{0}. The time span from

*t*

_{0}to the nominal position of the first bin,

*t*

_{1}, 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*

_{0}, and the time of the first reflectance value obtained with the binning process,

*t*

_{1}.

*T*, we fit a 3rd degree Non Uniform Rational B-spline (NURBS) [28,29] curve to each

*R*(

*r*= const,

_{k}*t*) curve. After calculating the approximating NURBS curve we resample it to obtain

_{l}*R*(

_{r}*r*= const,

_{k}*t*=

_{l}*T*). The values of

*R*(

_{r}*r*,

_{k}*t*) obtained after resampling is the final set used in the sMC

_{l}*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*

_{i}*R*(

_{r}*r*,

_{k}*t*) as native points for the interpolation necessary to calculate

_{l}*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*) as a 2D surface.

_{l}## 4. Results

### 4.1. cMC Simulations and Effects of Discretization

*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

*μ*∈ [10

_{a}^{−3}, 0.3]mm

^{−1}and

*μ′*∈ [0.5, 2]mm

_{s}^{−1}

*σ*(

_{rel}*r,t*) for the case of

*μ′*= 1mm

_{s}^{−1}and

*μ*= 0.1mm

_{a}^{−1}. Figure 3 shows that the values of

*σ*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

_{rel}*μ*≥ 0.1 mm

_{a}^{−1}. 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%.

*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

*σ*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.

_{rel}### 4.2. Comparison of single Monte Carlo Techniques

*and sMC*

_{p}*approaches established in §3 as compared to the reflectance provided by a cMC simulation*

_{i}*R*

_{cMC}(

*r*,

_{k}*t*) represented by a set of 200 × 200 radial/temporal bins formed using EFD.

_{l}*data set to evaluate the values (*

_{p}*r*,

_{i}*t*,

_{i}*w*) from (

_{i}*ρ*,

_{i}*τ*) using the scaling imposed by the value of

_{i}*μ*and modifying the photon weight in media with

_{s}*μ*≠ 0 by applying Beer’s law using the values of the photon exit time

_{a}*t*. To calculate

_{i}*R*

_{sMCp}(

*r*,

_{k}*t*) 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

_{l}*R*

_{sMCi}(

*r*,

_{k}*t*). 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).

_{l}*and (c) sMC*

_{p}*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*

_{i}*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*

_{p}*ɛ*

_{rel,sMCp}(

*r*,

_{k}*t*) over the entire physical range (

_{l}*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.*

*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.*

_{i}*μ*values smaller than 0.1 mm

_{a}^{−1}across all

*μ′*values, the mean relative error of the sMC

_{s}*approach tallied over the entire temporal and radial range excluding the first temporal bin is smaller than that obtained with sMC*

_{i}*method. This reflects the advantage gained by interpolating the reference signal which reduces the impact of statistical noise. For*

_{p}*μ*values equal to 0.1 mm

_{a}^{−1}and greater, we see that while the mean relative error for the reflectance values provided by the sMC

*approach is larger than that provided by the sMC*

_{i}*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*

_{p}*μ*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.

_{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

_{a}^{−1}and approaches 100% for

*μ*= 0.3 mm

_{a}^{−1}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

*μ*values of 0.1 mm

_{a}^{−1}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.

*method for large distances and long times decreases with increasing absorption. For*

_{i}*μ*= 0.001 mm

_{a}^{−1}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

*μ*values, the relative error increases for large distances and late times. Larger

_{a}*μ′*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.

_{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

_{s}*μ*= 0.1 mm

_{a}^{−1}and

*μ′*= 0.5 mm

_{s}^{−1}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

*approach typically reduces the relative error by an order of magnitude and provides accuracy comparable to that obtained using the sMC*

_{i}*approach. This provides a reflectance signal with excellent accuracy even when the optical properties chosen represent a large perturbation from the reference reflectance.*

_{p}## 5. Conclusion

*μ*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.

_{a}*method can be used to derive*

_{i}*R*(

*r,t*) for a wide range of optical properties and source detector separations. Given its speed, the sMC

*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/.*

_{i}## Acknowledgments

## 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 |

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. |

3. | J. Spanier and E. M. Gelbard, |

4. | A. D. Kim, “Transport theory for light propagation in biological tissue,” J. Opt. Soc. Am. A |

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. |

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 |

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. |

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. |

9. | B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and distributions of light in tissue,” Med. Phys. |

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. |

11. | A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol. |

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. |

13. | E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. |

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. |

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. |

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. |

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. |

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. |

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. |

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. |

21. | S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in multi-layered tissue,” Comput. Methods Programs Biomed. |

22. | W. Cheong and S. A. Prahl. “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. |

23. | L. G. Henyey and J. L. Greenstein, “Diffuse radiation of the galaxy,” Astrophys. J. |

24. | W. Hines, C. D. Montgomery, M. D. Goldsman, and M. C. Borror, |

25. | J. Dougherty, R. Kohavi, and S. Mehran, “Supervised and unsupervised discretization of continuous features,” in |

26. | M. Boulle, “Optimal bin number for equal frequency discretization,” Intell. Data Anal. |

27. | M. H. Kalos and P. A. Whitlock, |

28. | W. Tiller and L. Piegl, |

29. | F. D. Rogers, |

**OCIS Codes**

(170.3660) Medical optics and biotechnology : Light propagation in tissues

(170.5280) Medical optics and biotechnology : Photon migration

(170.7050) Medical optics and biotechnology : Turbid media

**ToC Category:**

Medical Optics and Biotechnology

**History**

Original Manuscript: August 4, 2011

Revised Manuscript: August 23, 2011

Manuscript Accepted: August 25, 2011

Published: September 22, 2011

**Virtual Issues**

Vol. 6, Iss. 10 *Virtual Journal for Biomedical Optics*

**Citation**

Michele Martinelli, Adam Gardner, David Cuccia, Carole Hayakawa, Jerome Spanier, and Vasan Venugopalan, "Analysis of single Monte Carlo methods for prediction of reflectance from turbid media," Opt. Express **19**, 19627-19642 (2011)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-20-19627

Sort: Year | Journal | Reset

### References

- 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. A14, 246–254 (1997). [CrossRef]
- 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]
- J. Spanier and E. M. Gelbard, Monte Carlo Principles and Neutron Transport Problems (Dover Publications, 2008).
- A. D. Kim, “Transport theory for light propagation in biological tissue,” J. Opt. Soc. Am. A21, 820–827 (2004). [CrossRef]
- 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]
- 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. A11, 2727–2741 (1994). [CrossRef]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- E. Alerstam, S. Andersson-Engels, and T. Svensson, “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt.13, 041301 (2008). [CrossRef]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- W. Cheong and S. A. Prahl. “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron.26, 2166–2185 (1990). [CrossRef]
- L. G. Henyey and J. L. Greenstein, “Diffuse radiation of the galaxy,” Astrophys. J.93, 70–83 (1941). [CrossRef]
- W. Hines, C. D. Montgomery, M. D. Goldsman, and M. C. Borror, Probability and Statistics in Engineering (John Wiley & Sons Inc., 2003).
- 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).
- M. Boulle, “Optimal bin number for equal frequency discretization,” Intell. Data Anal.9, 175–188 (2005).
- M. H. Kalos and P. A. Whitlock, Monte Carlo Methods Volume 1: Basics (Wiley-Interscience, 1995).
- W. Tiller and L. Piegl, The NURBS Book (Springer, 1995).
- F. D. Rogers, An Introduction to NURBS With Historical Perspective (Morgan Kaufmann, 2004).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.