1. Introduction

Numerical studies of light propagation in complex optical media such as biological tissue [1]-[5] are often performed through the Monte-Carlo simulation approach due to lack of analytic solutions to the radiative transfer equation [6]. Although Monte-Carlo modeling has a broad generality and can cope with arbitrary media complexity, it suffers from poor efficiency; substantial computational effort is required in order to obtain acceptable statistical precision. The poor efficiency of classical Monte-Carlo (CMC) simulations is particularly detrimental in problems involving back-reflected light from biological media having multiple-scattering characteristics but low backscattering. Approaches to combat this drawback involve the use of variance reduction techniques which reduce the estimator variance of the scored physical quantity, thus allowing a given statistical precision to be achieved using fewer light propagation paths. Among these methods are implicit photon capture [5] and importance sampling (IS) [2].

In the implicit photon capture technique, many photons (i.e., a packet of photons) rather than single photons are propagated simultaneously along the tissue, therefore improving the efficiency of the CMC simulation by eliminating the all-or-none nature of the single photon absorption events [5]. For the importance sampling method, the trajectories of photons propagating in the random medium are sampled from an importance proposal distribution, which tends to coerce the event of interest to occur more frequently. Typical events include the remitted radial distance of the photon from the laser source or the location from which the photon emerges given that it has escaped from the tissue surface. The fundamental difficulty with IS is the design of a “good” importance proposal distribution. Choosing a “good” distribution can result in significant run-time savings; the penalty for selecting a “bad” distribution can be longer run times than for the CMC simulation. Common “good” importance proposal distributions which force incident photons to penetrate deeper into the medium, while simultaneously enhancing backscattering are based on physical intuition and involve biasing the scattering lengths and the polar angles at certain steps throughout the simulation [2]. However, as the complexity of the light transport model in the biological medium increases, it becomes harder to design suitable photon trajectory distributions a priori.

In this paper, we apply the multicanonical Monte-Carlo (MMC) technique [7] to simulations of light propagation in biological media. The MMC method was originally designed to study first-order phase transitions [7], and it has been recently applied to compute also polarization mode dispersion statistics [8], bit error rates in optical communication systems [9] and noise statistics of semiconductor optical amplifiers [10]. Like ordinary IS, the MMC algorithm increases the number of rare events of interest by sampling the photon trajectories from an importance proposal distribution. However, the advantage of MMC is that it approaches the optimal importance distribution iteratively with little (if any) a priori knowledge of how to bias. The sampling from the sub-optimal importance distribution used within a specific MMC iteration is performed through random walks. To demonstrate the power of the MMC simulation method, we highlight simulations for which the CMC techniques are computationally inefficient. For brevity in the remainder of the manuscript, we will designate the acronym CMC to represent the classical Monte-Carlo simulations that propagate photon packets rather than single photons. Our results include estimating the rare event probabilities of the radial distance of photons exiting a semi-infinite homogeneous tissue as well as the rare event probabilities of the maximum penetration depth of photons propagating in an inhomogeneous tissue. The improved efficiency of the MMC algorithm suggests that multicanonical sampling will be useful in more complex optical forward and inverse problems for which a large space of feasible photon trajectories is to be explored and measures of resolution and adequate statistical precision of the solutions are desirable.

This paper is organized as follows. In Section 2, we outline the theoretical framework of the MMC algorithm and in Section 3 we describe its applicability to problems related to light propagation in biological media. Section 4 is devoted to the simulation results and discussion. Finally, conclusions are drawn is Section 5.

2. Multicanonical Monte-Carlo (MMC) algorithm

The MMC method [7] is a numerical technique that like IS algorithms increases the number of samples in the tail region of the computed probability density function (pdf) by sampling more frequently the values that have higher impact on the parameter being estimated. This biased sampling of the important values reduces the variance of the calculated pdf estimator. However, in contrast to standard IS where one has to intuitively choose an appropriate bias that will encourage sampling of the important values, the MMC algorithm uses an adaptive iterative procedure to select this bias with little (if any) a priori knowledge of how to bias. The iterative procedure updates the bias for the subsequent iteration by drawing samples trough a random walk which is performed in the important regions of the state space. As the iterations evolve the bias histogram tends to equilibrate at the optimal bias. Next, we describe the mathematical framework of the MMC method.

Let f_X (x) represent a given pdf defined on a high-dimensional space, χ, that describes the set of possible configurations of a system. In light propagation simulations this set corresponds to the possible trajectories of photons traveling in the random medium. Also, let the random variable Y(x) describe the scored physical quantity and let L_i represent the event {y_i < Y < y_i + d_y }. Then, the probability of the event L_i is simply

where f_Y (y_i ) is the pdf of Y evaluated at y_i and 1L_i (x) is the indicator function of the event L_i

Hence, f_Y (y_i ) can be estimated using the CMC estimator which reads as

where {x_n${\}}_{n=1}^N$ are N independent identically distributed (iid) samples drawn from f_X (x). Eq. (3) indicates that the CMC estimator simply computes the percentage of samples (out of the N drawn samples) corresponding to the event L_i . In order to reduce the variance of ${\widehat{f}}_Y^{\mathit{\text{CMC}}}$ (y_i ), the MMC method employs the importance sampling (IS) estimation principle [11], [12] and introduces an importance proposal distribution, q(x). This allows Eq. (1) to be written as

where w(x)= f_X(x) / q(x) is the so called importance weight [11], [12] and the corresponding IS estimator of f_Y(y_i) is expressed as

with x_n being generated now from q(x). Notice that the IS estimator use the importance weight w(x_n) in order to correct for the use of samples x_n drawn from the biased distribution q(x) -rather than from the original distribution f_X(x) - therefore ensuring that the pdf estimation is unbiased alike the CMC estimator. Also, note that Eqs. (5) and (3) are identical when q(x)= f_X(x) and that the optimal importance distribution which minimizes the variance of the IS estimator reads as [11], [12]

However, the optimal importance distribution is not very useful since it depends on f_Y(y_i), i∈ N which is unknown. Hence, while the aim in IS is to design a “good” approximation of q ^opt (x) based on physical intuition, the idea behind the MMC algorithm is to iterate over sub-optimal importance proposal distributions that converge toward q ^opt (x). Specifically, the sub-optimal importance proposal distribution used in the k-th iteration of the MMC method, q ^opt(k) (x) is given by

where $f_Y^{\mathit{(}k\mathit{-}\mathit{1}\mathit{)}}$ (y_i) is the estimation of the pdf of Y obtained in the previous iteration. The sub-optimal importance proposal distribution q ^opt(k) (x) tends to equilibrate at the optimal bias q ^opt (x) as the iteration number increases. $f_Y^{\mathit{\left(}\mathit{0}\mathit{\right)}}$ (y_i), i∈ N can be a uniform distribution or can be calculated by the CMC estimator (Eq. (3)) using a small number of samples. Within the k-th MMC iteration, the Markov chain Monte Carlo (MCMC) strategy [12] and specifically the Metropolis algorithm [13] is employed to generate samples x_n from q ^opt(k) (x). A Metropolis step involves proposing a candidate sample $x_n^{*}$= x_n + δx where x_n is the current sample and δx is sampled from a symmetric distribution (e.g., uniform distribution with zero mean). The candidate sample is then accepted as the next sample with probability min{1, q ^opt(k) (${\mathbf{x}}_n^{*}$) / q ^opt(k) (x_n)=[f_X($x_n^{*}$) $f_Y^{(k-1)}$ (y(x_n))] / [f_X(x_n) $f_Y^{(k-1)}$ (y(${\mathbf{x}}_n^{*}$))]} otherwise the next sample remains x_n. Notice that by using the Metropolis method we only need to know q ^opt(k) (x) up to a constant of proportionality. Essentially, the Metropolis algorithm generates samples from q ^opt(k) (x) by exploring the state space χ using a random walk, or equivalently a homogeneous, irreducible aperiodic Markov chain which employs the detailed balance condition [12] to converge toward the invariant distribution q ^opt(k) (x) independently of the starting state. Consequently, the success or failure of the algorithm hinges on the choice of distribution of δx as well as on the number of Metropolis steps (N) and the number of MMC iterations that are simulated. Note that the selection of these parameters is performed empirically such that each MMC iteration covers a larger values range of the simulated random variable with adequate standard deviation (including those values located in the pdf’s tail) and that the last MMC iteration achieves the desired probability level. Finally, the updated estimation of the pdf of Y evaluated at y_i , i∈ N which results from the k-th MMC iteration is computed using Eq. (5) and reads as

with C being a normalizing constant ensuring that ∫_{-∞ ∞} $f_Y^{\left(k\right)}$ (y) dy = 1.

To conclude, the MMC technique uses the IS principle (Eq. (5)) together with sub-optimal importance proposal distributions that are computed iteratively (Eq. (7)) to approach the optimal importance proposal distribution (Eq. (6)). An N step Metropolis algorithm is constructed to mimic samples drawn from these sub-optimal importance proposal distributions.

3. Application of the MMC method to light propagation in biological media

The set of possible photon packet trajectories in a biological medium is determined by the polar angles {θ_m}, the azimuthal angles {ϕ_m}, the scattering lengths {l_m}, and the reflection events {r_m} at the boundaries of layers having distinct refractive indices. m = 0,1, ⋯,s, where s indicates the index of the last scattering event. Let f_Θ(θ_m),f_Φ(ϕ_m), f_L(l_m) and f_R(r_m) represent the pdf’s of the iid random variables θ_m , ϕ_m , l_m and r_m , hence f_X (x) = ${\prod }_{m=0}^s$ f _Θ(θ_m )f _Φ(ϕ_m )f_L (l_m )f_R (r_m )where x = (θ₁ , ⋯,θ_s ,ϕ₁ , ⋯,ϕ_s ,l₁ , ⋯, l_s , r₁ , ⋯, r_s , and χ=(Θ,Φ,L,R). Also, let y(x) denote the event of interest for which we wish to compute the pdf. Without loss of generality, we concentrate on events which correspond to photons exiting the medium. Therefore, the computed pdf’s are essentially conditional pdf’s of y(x) given that the photon packet exits the medium. Similar principles can be applied for the events related to photons being absorbed along their propagation in the medium. Finally, notice that we denote conditional pdf’s and pdf’s similarly throughout this work in order to simplify notations and the correct interpretation of the pdf is rather understood from the context.

Simulations of photon packet trajectories in biological media are more efficient than single photon trajectories [4], [5]. Hence, we first employ the MMC algorithm to estimate f ^no_abs _Y,Wpacket (y,w _packet) (that is,f ^no_abs _{Y,Wpacket|EXIT}(y,w _packet | EXIT)), which is the joint pdf of y(x) without absorption and the corresponding photon packet weight, w _packet (x), given that the photon packet exists the medium, and then we use the relationship

to compute the pdf of y (x) with absorption. This approach significantly contributes to the overall efficiency of the MMC method and is preferable over MMC simulations of single photon trajectories which directly estimate f_y(y).

In view of the fact that f_X(x) is a separable function, each Metropolis step within the k-th MMC iteration is implemented as follows. We first perturb separately θ_m , ϕ_m , l_m and r_m , m=0,1, ⋯,s by dθ ~ U(-ε_θ / 2,ε_θ /2), dϕ ~ U(-ε_ϕ / 2,ε_ϕ /2), dl ~ U(-ε_l / 2,ε_l /2) and dr ~ U(-ε_r /2,ε_r/2), respectively, with U symbolizing uniform distribution and ε_θ , ε_ϕ , ε_l , ε_r being selected such that most of the estimated pdf modes are visited equivalently. Perturbation coefficients leading to an acceptance ratio - defined as the ratio of the accepted Metropolis steps to the total number of Metropolis steps (N) - of 30%-45% were found to be appropriate in our studies. Next, we accept or reject independently each of the proposed ${\theta }_m^{*}$ = θ_m +dθ, ${\varphi }_m^{*}$ = ϕ_m + dϕ, $l_m^{*}$ = l_m + dl and $r_m^{*}$ = r_m + dr, m=0,1,⋯,s with probability min{1, f_Θ${\mathit{(}\theta }_m^{\mathit{*}}$) / f_Θ(θ_m)}, min{1,f_Φ${\mathit{(}\varphi }_m^{\mathit{*}}$) / f_Φ (ϕ_m)}, min{1, f_L${\mathit{(}l}_m^{\mathit{*}}$) /f_L (l_m)}, and min{1, f_R${\mathit{(}r}_m^{\mathit{*}}$) / f_R (r_m)}, respectively. The fact that ϕ_m and r_m are uniformly distributed and are subjected to periodic boundary conditions results in accepting ${\varphi }_m^{*}$ and $r_m^{*}$ with probability one. Following the selection of the test trajectory x ^* , y (x^*) is then calculated by propagating the photon packet along the biological media as described by Wang et al. [5]. Finally, the entire Metropolis step from x to x ^* is accepted with probability min{1,f ^no_abs(k-1) _Y,Wpacket (y(x),w _packet (x)) /f ^no_abs(k-1) _Y,Wpacket (y(x^*),w _packet (x^*))}. Note that the Metropolis algorithm is designed to accept only steps which result in photon packets that exit the medium, hence estimating the conditional pdff ^no_abs _{Y,Wpacket|EXIT} (y,w _packet | exit).

It is worth mentioning that the described perturbation scheme explores the high dimensional space χ more efficiently than the method presented in section 2 due to the independent perturbations of the pdf‘s comprising f_X(x). However, notice that this technique is applicable only in cases for which f_X(x) is a separable function.

Following N Metropolis steps, the joint pdf f _Y,wPacket(y,w _{packe t}) is updated similarly to Eq. (8):

where L_ij represents the event {y_i < Y < y_i +dy, w _packetj < W _packet < W _packetj+dw _packet} and C stands for a normalizing constant ensuring that ${\int }_{-\infty }^{\infty }$${\int }_0^1$ f̂^no_abs(k) _Y,Wpacket(y_i , w _packet) dw _packet dy = 1. Finally, the estimated f_Y(y) obtained from the k-th MMC iteration is

Figure 1 summarizes the pseudo code of the MMC algorithm applied in computations of pdf’s emerging in simulations of light transport in biological media using the Heyney-Greenstein pdf for the cosine of the deflection angle [5].

 

Fig. 1. Multicanonical Monte-Carlo (MMC) algorithm for light transport in biological media.

Download Full Size | PDF

4. Simulation results and discussion

Two scenarios of light propagation in a biological medium are examined in this section. These scenarios serve as demonstrations for the potential use of the multicanonical sampling in stochastic simulations of cases for which CMC methods are impractical or computationally inefficient.

The results presented hereon were simulated using Matlab software running on a conventional Pentium 4 desktop computer.

4.1 Radially resolved steady state diffuse reflectance pdf for photons propagating in a semi-infinite homogeneous random medium

To illustrate the efficacy of the MMC method compared with that of CMC simulations, we first used multicanonical sampling to compute the pdf of the radial distance of photons exiting a semi-infinite homogeneous tissue with the following optical properties: μ_a = 0.25 mm^-1, μ_s = 50 mm^-1 and g = 0.9. We assumed that the tissue-ambient medium boundary is refractive-index matched. In the simulations, the photon packets were launched orthogonally onto the tissue at the origin. Tracking a photon packet was terminated either when the photon packet exited the medium, the weight of the photon packet was sufficiently low (< 10^-15) or when the number of scattering events exceeded a predefined value. For this study, the maximum number of scattering events was 2000, resulting in a state space χ of dimension (1999 polar angles) × (1999 azimuthal angles) × (2000 scattering lengths) and ensuring unbiased low probability values which were further confirmed by trial simulations using 10000 scattering events. Furthermore, the joint pdf of the radial distance without absorption and the corresponding photon packet weight is scored using a linear scale for the radial distance axis and a semi-logarithmic scale for the packet weight axis.

The MMC simulations employed six iterations which were found to be sufficient in order to estimate diffuse reflectance probability values as low as ~10^-8 with standard deviation of one half-decade. The first iteration consisted of 2187 photon packets and used a CMC simulation which introduced the first coarse estimation for the diffuse reflectance pdf. Note that alternatively one can set f ^no_abs(0) _Y,Wpacket (y(x),w _packet (x)) = 1/D with D = [0, y_max ] × [0,1] and execute an MMC iteration. These two approaches for computing the first coarse estimation of the diffuse reflectance pdf showed similar results. The remaining iterations comprised five multicanonical sampling procedures with 15000, 22500, 33750, 50625 and 75938 photon packets. Accordingly, a total number of 2·10⁵ photon packets was simulated during 80 minutes. Notice that the number of photon packets used in each iteration was increased in order to encourage the Metropolis random walk to accept samples also in the tail of the estimated conditional pdf and not only at shorter radial exit distances. Furthermore, the number of photon packets in each iteration was determined such that it yielded an improvement of approximately one and a half orders of magnitude in the computed diffuse reflectance pdf values.

Figure 2 shows four (out of six) of the MMC iterations. It includes also the analytic expression of the diffuse reflectance pdf [14] in dashed line. We observe that the estimation of the diffuse reflectance pdf resulting from the first MMC iteration is limited to values higher than 10^-3 as indicated by the straight horizontal line appearing at the pdf tail. This line places a threshold of one half-decade for the fluctuations of the estimated pdf. Note that for the tissue optical parameters used in this example, pdf values higher than 10^-3 correspond to radial exit distances lower than 2 mm. The subsequent MMC iterations cover a larger radial exit distance range which finally extends over approximately 8 mm. Notice that the last MMC iteration considerably reduces the estimation variance of the preceding MMC iterations and generates a pdf estimator which successfully follows the tail of the theoretical diffuse reflectance pdf down to values of ~10^-8. To conclude, Fig. 2 illustrates the usefulness of multicanonical sampling in increasing the dynamic range of the estimated pdf region with adequate statistical accuracy. Moreover, each MMC iteration improves the tail estimation in approximately one order of magnitude while maintaining a similar order of simulated photon packets number.

 

Fig. 2. MMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Download Full Size | PDF

To compare the performances of the MMC and CMC techniques, we performed a CMC simulation employing 2·10⁵ photon packets - a number identical to that simulated previously by the MMC method. Figure 3 presents the estimated diffuse reflectance pdf obtained from the CMC simulation (in red color) and MMC simulation (in green color) using an identical number of photon packets. Once more, the dashed line describes the theoretical pdf calculated by diffusion theory [14], and the straight horizontal lines represent the minimum achievable pdf value for which the fluctuations of the estimated pdf have not exceed a threshold of one half-decade. Clearly, the MMC estimator of the diffuse reflectance pdf is significantly superior to the CMC estimator - it covers a larger range of radial exit distances and approaches pdf values that are more than two orders of magnitude smaller while preserving an acceptable estimation variance over the entire radial exit distance range.

It is important to point out that the CMC simulation ran 2.5 times faster than the MMC simulation, even though both simulation codes employed an identical photon transport function and had an identical number of photon packets. This difference in the run time stemmed mainly from the fact that longer photon trajectories should be computed throughout the MMC iterations in order to achieve lower values of the diffuse reflectance pdf. This point is further confirmed by executing a CMC simulation using 3.9·10⁶ photon packets as shown in Fig. 3 in blue color. While both CMC and MMC techniques attained similar pdf values as low as 10^-8 with similar statistical accuracy, the CMC simulation required nine hours to complete - 6.75 times slower than the MMC algorithm. These results illustrate the improved performance of the MMC algorithm in estimating adequately the diffuse reflectance pdf over a larger range of radial exit distances.

 

Fig. 3. MMC and CMC estimators of the diffuse reflectance pdf of a semi-infinite homogeneous random medium.

Download Full Size | PDF

4.2 pdf of the maximum penetration depth of photons that propagate and exit a two-layered random medium

To further explore the performance of multicanonical sampling of photon packets propagating in biological media, we performed MMC simulations for a two-layered random medium. In these simulations, a pencil beam of photons entered the top layer of the medium along the z direction and traveled inside the inhomogeneous media. Then, the maximum penetration depth Z_max of photon packets that successfully exited the top layer of the medium (i.e., Z_max | EXIT) was recorded and its pdf was computed. Tracking a photon packet was terminated either when the weight of the photon packet was sufficiently low (< 10^-15) or when the number of scattering events exceeded a predefined value. For this study, the maximum number of scattering events was 3000 (therefore guaranteeing unbiased pdf estimators) and a semi-logarithmic grid was used for the packet weight axis of the joint pdf of Z_max without absorption and the corresponding photon packet weight.

The scattering properties of the two layers were identical, that is μ_s1 = μ_s2 = 50 mm^-1 and g₁ = g₂ = 0.9, whereas the absorption coefficients of the first and second layer were μ_a1 = 0.25 mm^-1 and μ_a2 = 0.025 mm^-1, respectively. Furthermore, the top-layer thickness was 3.75 mm and the refractive indices for the two layers as well as for the overlying ambient medium were equal.

First, we confirmed that the MMC estimation for the diffuse reflectance pdf is consistent with the analytic result [15]. Next, six MMC iterations were executed during 80 min using a total number of 2·10⁵ photon packets in order to compute the pdf of the random variable Z_max | EXIT. Figure 4 presents the last iteration of the MMC pdf estimation for Z_max | EXIT (in green color). The straight horizontal line at the pdf tail shows that pdf values lower than 10^-7 (for which the fluctuations of the estimated pdf have not exceed a threshold of one half-decade) were obtained throughout the simulation. Furthermore, two distinct slopes are clearly observed, where the slope for 0 < Z_max < 3.75 mm is larger than that for Z_max > 3.75 mm. This result is expected and stems from the lower absorption coefficient of the second layer.

 

Fig. 4. MMC and CMC estimators of the pdf of the maximum penetration depth of light in a two-layered inhomogeneous random medium.

Download Full Size | PDF

We compared the performance of the multicanonical and Monte-Carlo sampling methods by carrying out two CMC simulations. The first CMC simulation lasted 30 minutes and estimated the pdf of the maximum penetration depth using 2·10⁵ photon packets (Fig. 4, red dashed line). The second CMC simulation employed 3.75·10⁶ photon packets and was completed in eight hours (Fig. 4, in blue line). We note that with the CMC method, simulating 2·10⁵ photon packets is not adequate to predict the correct pdf for the Z_max | EXIT. Moreover, the number of photon packets should be increased by at least one order of magnitude in order to obtain an acceptable estimation of this pdf. However, the improved estimation variance comes at the expense of computational time which in this case is six-fold slower than that achieved by multicanonical sampling for a similar level of estimation variance. This study therefore demonstrates the superior efficacy of the MMC algorithm over the CMC method in sampling back reflected light in order to successfully identify tissue inhomogeneities.

5. Conclusions

In propagation simulations, where importance proposal distributions can be intuitively determined a priori, importance sampling (IS) [2] will be highly efficient, since the photon trajectories that give rise to the tail-states are predetermined. In this paper, we have presented the application of an adaptive importance sampling technique based on multicanonical Monte-Carlo (MMC) simulations to the studies of light propagation in biological media. Like IS, MMC efficiently accesses the low-probability tails of numerous pdf’s of interest in biological light transport simulations. Unlike standard IS algorithms, however, the MMC method is blind to the physics of the light propagation model. As a result, MMC may be considered more general than IS allowing it to be applied to a wider variety of simulation geometries in which the importance proposal distributions can not be intuitively advised. It is important to point out that the efficiency of MMC sampling hinges on the choice of the perturbation coefficients used in the Metropolis random walk as well as on the number of the simulated Metropolis steps and MMC iterations. These parameters are selected empirically such that (1) most of the estimated pdf modes are visited equivalently through the random walk, (2) the desired probability level is obtained with adequate standard deviation, and (3) each MMC iteration successfully increases the dynamic range in the pdf tail region.

We have shown the efficacy of the MMC method by evaluating the pdf for the radial distance of photons exiting from a semi-infinite homogeneous medium as well as the pdf for the maximum penetration depth of photons propagating in an inhomogeneous medium. The computed pdf’s were consistent with those estimated by classical Monte-Carlo (CMC) simulations. Furthermore, in these examples a six- to seven-fold improvement in computational time was achieved using multicanonical sampling. These results illuminate the potential usefulness of the MMC algorithm for an efficient sampling of pdf’s in complex biological simulations. Based on these encouraging results we feel that further investigation of the computational efficiency of MMC for other simulation parameters_as well as extension of MMC to model temporal, polarization, and coherence effects are merited.