OSA's Digital Library

Biomedical Optics Express

Biomedical Optics Express

  • Editor: Joseph A. Izatt
  • Vol. 2, Iss. 5 — May. 1, 2011
  • pp: 1069–1081
« Show journal navigation

Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography

Ivan T. Lima, Jr., Anshul Kalra, and Sherif S. Sherif  »View Author Affiliations


Biomedical Optics Express, Vol. 2, Issue 5, pp. 1069-1081 (2011)
http://dx.doi.org/10.1364/BOE.2.001069


View Full Text Article

Acrobat PDF (1151 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We developed an importance sampling based method that significantly speeds up the calculation of the diffusive reflectance due to ballistic and to quasi-ballistic components of photons scattered in turbid media: Class I diffusive reflectance. These components of scattered photons make up the signal in optical coherence tomography (OCT) imaging. We show that the use of this method reduces the computation time of this diffusive reflectance in time-domain OCT by up to three orders of magnitude when compared with standard Monte Carlo simulation. Our method does not produce a systematic bias in the statistical result that is typically observed in existing methods to speed up Monte Carlo simulations of light transport in tissue. This fast Monte Carlo calculation of the Class I diffusive reflectance can be used as a tool to further study the physical process governing OCT signals, e.g., obtain the statistics of the depth-scan, including the effects of multiple scattering of light, in OCT. This is an important prerequisite to future research to increase penetration depth and to improve image extraction in OCT.

© 2011 OSA

1. Introduction

Optical Coherence Tomography (OCT) is rapidly becoming an important imaging technique for numerous medical and biological applications [1

1. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]

]. It is a sub-surface imaging technique that uses either a low-coherence light source (time-domain systems) or a wavelength-swept laser source (frequency-domain systems). It has one to two orders of magnitude higher resolution than ultrasound imaging and an imaging depth that can reach up to 3 mm, depending on the optical properties of the tissue, and can produce images inside the body when integrated with optical fiber probes. The use of infrared and visible light is safer to most biological samples than ionizing radiation like X-rays or gamma rays, and it also allows for spectroscopic characterization of an object, e.g., a tumor in tissue [2

2. C. Xu, C. Vinegoni, T. S. Ralston, W. Luo, W. Tan, and S. A. Boppart, “Spectroscopic spectral-domain optical coherence microscopy,” Opt. Lett. 31(8), 1079–1081 (2006). [CrossRef] [PubMed]

].

Frequency-domain OCT (FD-OCT) has been receiving increasing attention because of its higher imaging speed and higher signal-to-noise ratio compared to time-domain OCT (TD-OCT) [3

3. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef] [PubMed]

]. However TD-OCT systems can attain higher imaging resolution and dynamic range than their Fourier-domain counterparts [4

4. B. Liu and M. E. Brezinski, “Theoretical and practical considerations on detection performance of time domain, Fourier domain, and swept source optical coherence tomography,” J. Biomed. Opt. 12(4), 044007 (2007). [CrossRef] [PubMed]

]. In this paper we focus on TD-OCT imaging, but we are currently generalizing our results to FD-OCT systems. Monte Carlo simulation of FD-OCT imaging is more challenging than TD-OCT because the image of an object is not obtained through a relatively simple reflectance based model. Due to the nature of the Fourier transform, the image of every point in the object is obtained from many measurements involving the whole object.

OCT imaging considers the ballistic and quasi-ballistic photons reflected from a target layer inside the tissue, which we denote Class I diffuse reflectance, as the OCT signal that produces the image [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. However multiple scattered photons reflected from the tissue, which we denote Class II diffuse reflectance, do not carry useful information about the target layer and are considered a source of noise in OCT [6

6. S. S. Sherif, C. C. Rosa, C. Flueraru, S. Chang, Y. Mao, and A. G. Podoleanu, “Statistics of the depth-scan photocurrent in time-domain optical coherence tomography,” J. Opt. Soc. Am. A 25(1), 16–20 (2008). [CrossRef] [PubMed]

]. It has been demonstrated that Class II diffuse reflectance is the fundamental limitation in increasing the imaging depth of OCT in tissue [7

7. M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple scattering in optical coherence microscopy,” Appl. Opt. 34(25), 5699–5707 (1995). [CrossRef] [PubMed]

]. Understanding the physical process governing both Class I and Class II diffusive reflectance is an important prerequisite to any effort to increase penetration depth and to enhance the quality of the images obtained with OCT. As this physical process is complicated, Monte Carlo simulation of light transport in tissue [8

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

11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

] has been used to obtain the TD-OCT signal [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

].

We use Monte-Carlo simulation with importance sampling to calculate the TD-OCT signal because the probability that a photon propagating in typical biological tissue will undergo single-scattered backreflection is very low. This very low probability of events of interest would require an unacceptably large computational time if standard Monte Carlo simulation were used. Importance Sampling is an advanced statistical method [12

12. R. Y. Rubinstein, Simulation and the Monte Carlo Method (Wiley, 1981).

] that consists of biasing random events in such a way that the events of interest, which are often rare, appear more often in Monte Carlo simulations [13

13. G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(3), 310–312 (2002). [CrossRef]

19

19. I. T. Lima, Jr., “Advanced Monte Carlo methods applied to Optical Coherence Tomography” (invited), presented at the 2009 SBMO/IEEE MTT-S International Microwave and Optoelectronics Conference, Belém, Brazil, 3–6 Nov. 2009.

]. To produce the correct statistical result in Monte Carlo simulation with importance sampling, each biased sample is weighted by the likelihood in which this sample would have been observed in the unbiased simulation. This procedure, which has to be tailored to each particular application, can reduce the computational time of Monte Carlo simulations by several orders of magnitude.

Importance sampling, tailored to each particular application, has been applied in optical communications [13

13. G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(3), 310–312 (2002). [CrossRef]

16

16. I. T. Lima, A. O. Lima, G. Biondini, C. R. Menyuk, and W. L. Kath, “A comparative study of single section polarization-mode dispersion compensators,” J. Lightwave Technol. 22(4), 1023–1032 (2004). [CrossRef]

], confocal microscopy [17

17. J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A 13(5), 952–961 (1996). [CrossRef] [PubMed]

], atmospheric optics [18

18. H. Iwabuchi, ““Efficient Monte Carlo method for radiative transfer modeling,” J. of the Atmosph,” Science 63, 2324–2339 (2006).

] diffuse optical tomography (DOT) [11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

] and TD-OCT [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. In the applications involving Monte Carlo simulation of light transport in three-dimensional space, one limitation that we addressed is the statistical bias that increases with depth as the photons are preferentially scattered in the backward direction by importance sampling based implementations.

2. Previous modeling and simulation parameters

We started the Monte Carlo modeling of light transport in multi-layered tissues (MCML) with a C-language software package that is available for download from the web site of the Oregon Medical Laser Center [20

20. “Monte Carlo simulations,” Oregon Medical Laser Center, accessed January 1, 2009, http://omlc.ogi.edu/software/mc/

]. MCML allows the simulation of an ensemble of photon packets launched in a steady-state pencil beam, normal to the surface of the topmost layer. Each photon packet produces a random walk whose step size is determined by an exponentially distributed random variable defined by the interaction coefficient, which is equal to the sum of the absorption µa and the scattering µs coefficients. The scattering events, which take place at the end of the random steps, are produced by two random angles that determine the future direction of the photon packet scattering in three-dimensional space . To account for the photon packet scattering with arbitrary anisotropy factor, g, we use the same Henyey-Greenstein probability density function used in the MCML software package that is defined as

fHG(cosθs)=1g22(1+g22gcosθs)3/2,
(1)

The Class I diffuse reflectance at the depth equal to z is obtained by calculating the mean value of the indicator function I 1 of the spatial and temporal filter of the Class I diffuse reflectance for all the photon packets (samples) in the ensemble. The indicator function I 1 of the spatial and temporal filter for the ith photon packet is defined as

I1(z,i)={1  ri<dmax,  θs,i<θmax,  |Δsi2z|<lc0  otherwise,
(2)

where lc is the coherence length of the source, ri is the distance of the ith reflected photon packet to the origin along the plane z = 0, where the collecting optical system is located, d max and θ max are the maximum collecting diameter and angle, respectively, θz,i is the angle with the z-axis, which is the axis normal to the tissue interface, Δsi is the optical path, and z is the maximum depth reached by the photon packet. The diffuse reflectance R 1 at any depth, which is the expected value of I 1 at that corresponding depth, and its corresponding standard deviation σR ,1 can be estimated by the expressions

R1(z)=1Ni=1NI1(z,i)W(i)
(3)

and

σR,12(z)=1N1i=1N[I1(z,i)W(i)R1]2,
(4)

where N is the initial number of photon packets in the Monte Carlo simulations [21

21. S. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory (Prentice-Hall, 1993).

].

3. Novel Importance sampling to simulate Class I diffuse reflectance in TD-OCT

First, we modified the MCML software package to implement the spatial filter, the angular filter, and the time gating for the diffusely reflected photons according to the optical fiber model specified. For simplicity, we also used a square time gating as in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. Then, we implemented our importance sampling method that we describe in this section. Since most tissue generally are highly forward scattering, the value of their anisotropy factor is close to 1. Thus there is a very small probability that a simulated photon packet at a given depth undergoes backscattering towards the tip of the collecting fiber. This already small probability of collecting Class I photons decreases rapidly with the depth from which the photon is backscattered. To speed up the simulation of collected Class I photons, we developed a novel importance sampling method for Monte Carlo simulations that bias the scattered photon packet direction u^' preferentially towards the apparent position of the center tip of the collecting optical system v^ once the photon packet reaches the depth range of interest. If a focusing lens is used, the bias direction v^ will be defined as the direction of the apparent position of the center of the lens. We define the center of the Cartesian coordinate system at the tip of the collecting optical system. The bias direction in which the collecting optical system is located is defined as v^=R/|R|, where R=xx^+yy^+zz^ is the vector position vector of the scattering point in the tissue.

The photon packets propagating in a direction closer to v^ will have a higher probability of contributing to the collected Class I diffuse reflectance. Therefore our bias direction is more efficient than biasing only in the backward direction as previously done in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

], since that direction may not be consistent with the apparent direction of the collecting fiber. This reasoning is particularly true deeper in the tissue, where photon packets suffer one or more forward scattering events that deviate them from their original direction of propagation.

3.1 Calculation of the apparent position of the collecting optics

To determine the bias direction of the photon scattering in any of the tissue layers, we first calculate the apparent position of the detector in each layer using the paraxial approximation. The apparent position of the detector changes with the layer because of refraction in layers of tissue with different refractive indices. In the topmost layer, in which the detector is located, the apparent position of the detector is equal to the actual position of the detector: Rd,1=0x^+0y^+0z^. We calculate the apparent position in each layer following a recursive process starting from the second layer until the last layer. The apparent z-coordinate of the detector in the jth layer, zd,j', is given by

zd,j'=zL,jnjnj1(zL,jzd,j1'),
(5)

where zL,j is the topmost coordinate of the jth layer and nj is the refractive index of the jth layer. We obtained Eq. (5) from the law of refraction using the paraxial approximation: sinθtanθ [22

22. E. Hecht, Optics, 4th ed. (Pearson Addison Wesley, 2003).

]. In the paraxial approximation, a photon packet located in the layer j would arrive at the origin (0,0,0) in the absence of additional scattering if it is directed towards the coordinate (0,0, zd,j'). The bias direction produced by a scattering located at the position R located at the jth layer is defined as

v^=zd,j'z^R.
(6)

3.2 Calculation of scattering angle of the first backscattering event

fB(cosθB)=1a22(1+a22acosθB)3/2,
(7)

where a is a bias coefficient. After randomly picking a biased angle θB away from the biased directionv^, so that cosθB=v^u^', the resultant biased scattering direction u^' is rotated around the biased direction v^ by an angle ϕ that is randomly picked from a uniform probability density function from 0 to 2π. This last procedure is equivalent to the one used in MCML software package to enable a full three-dimensional scattering, except that the scattered angle θB in Eq. (7) is defined with respect to the biased direction v^, as opposed to the direction u^ in which the photon packet was propagating prior to that scattering event. This procedure produces the new propagation direction u^' of the photon packet after the first biased scattering. Then the scattered photon packet is associated with a quantity that is defined as the likelihood ratio in the importance sampling formalism [13

13. G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(3), 310–312 (2002). [CrossRef]

16

16. I. T. Lima, A. O. Lima, G. Biondini, C. R. Menyuk, and W. L. Kath, “A comparative study of single section polarization-mode dispersion compensators,” J. Lightwave Technol. 22(4), 1023–1032 (2004). [CrossRef]

]. The likelihood ratio of the photon packet using our probability density function of the biased angle, Eq. (7), is given by

L(cosθB)=fHG(cosθS)fB(cosθB)=1g21a2(1+a22acosθB1+g22gcosθS)3/2,
(8)

where cosθS=u^u^' is a calculated variable, since it is a function of both the random values of θB and ϕ. The likelihood ratio in Eq. (8) is the ratio of the probability that this biased scattering angle would have been observed in an unbiased simulation divided by the probability of the biased scattering angle specified by the biased distribution. While θBis the angle between the scattered direction u^' and the bias direction v^, the likelihood ratio also depends on the angle θS between the scattered direction u^' and the incoming propagation direction u^, which is the direction of the forward propagation towards which the photon packet was much more likely to have been scattered. Figure 1
Fig. 1 Schematic representation of the vectors and the angles used in the bias procedure.
shows a schematic representation of these vectors and the angles used in this bias procedure. Since the actual choice of the bias distribution only affects the speed of convergence of the calculation, other biased probability function can also be used to randomly generate the biased scattering towards the bias direction v^. The Henyey-Greenstein probability density function that we used in Eq. (8), whose bias strength determined by the coefficient a, enabled a rapid convergence of the Importance Sampling method for the cases that we reported in Sec. 4 as a result of the strong bias used.

3.3 Calculation of scattering angles of further backscattering events

Once a photon packet is biased towards the apparent position of the collecting optical system at any given depth in which the photon packet interacts with the tissue, the photon packet becomes most likely to be collected at the tip of the fiber. However, the photon packet can be scattered several times in the tissue after the first backscatter bias before reaching the collecting optical system. These additional scatterings, while most likely being in the forward direction, according to the standard scattering function in Eq. (1), decrease the correlation between the biased direction and the event in which the photon packet is collected by the tip of the fiber. We overcome this lack of correlation by continuing to bias the scattering direction u^' towards the direction v^, which points to the apparent position of the collecting optical system, at every scattering point until the photon packet is terminated. These additional biases still use the probability density function of the biased angle and the corresponding likelihood ratio given in Eqs. (7) and (8), respectively. Because the values picked for the angle between the scattering direction and the biased direction are independent from one scattering event to the next, the total likelihood ratio is equal to the product of all the likelihood ratios of all biased scattering point until the photon packet is terminated.

After undergoing the first biased scattering, each simulated photon packet is biased at every additional scattering point until it is terminated, which happens when the photon packet either leaves the tissue or undergoes so many scatterings that it becomes eliminated by a Russian roulette procedure, described in Sec. 2. At the end of the ensemble of N Monte Carlo simulations, with our importance sampling, the diffuse reflectance R 1 and its corresponding standard deviation σ R,1 can be estimated by the following expressions

R1(z)=1Ni=1NI1(z,i)L(i)W(i)
(9)

and

σR,12(z)=1N1i=1N[I1(z,i)L(i)W(i)R1]2.
(10)

Equations (9) and (10) are similar to the Eqs. (3) and (4), except that the indicator function at any sample is multiplied by its corresponding likelihood ratio. This importance sampling procedure applied at a certain depth range, such as from 0 µm to 600 µm, or at any other specified depth range, allows a much larger number of photon packets to be scattered from that depth range towards the collecting optical system than what is obtained using a standard Monte Carlo method that is also based on the MCML. At the end of the simulation, each biased photon packet is weighted by its likelihood ratio, which adjusts the contribution of each sample to the estimation of the Class I diffuse reflectance, so that the estimated diffuse reflectance converges towards the true value several order of magnitudes faster than using the standard Monte Carlo method.

3.4 Comparison with a previously proposed bias method

This importance sampling implementation is comparable to the bias procedure described in [10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

], except in the way in which the biasing of the scattering angle is implemented. In our method, the angles are deterministically biased towards the preferred direction, and the scattered photon packets are weighted by the corresponding likelihood ratio that eliminates any residual statistical bias. In the method in [10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

], the scattering direction is chosen using the standard Henyey-Greenstein probability function in Eq. (1), but the actual scattering direction is rejected as many times as necessary until an acceptable direction closer to the direction of the source is chosen. Therefore, the method in [10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

] requires considerable additional computational time per biased scattering event. We note that this computational time increases with the bias strength, since the probability of rejection of the scattering directions increases with the amount of bias. Moreover, in this method a discrete, one-dimensional, lookup table of bias scatter weights has to be generated prior to the simulations to enable weighting the scattered photon packets that were biased. That process adds statistical bias to the results due to the approximation to the nearest value in the table of weights. Moreover, in [10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

] no procedure was provided to address the reduction in the number of photon packets that are penetrating the tissue due to the backward bias scattering, which limits the use of that technique to forward bias scattering.

3.5 Importance sampling statistics as a function of the depth

One limitation that results from the use of previously existing bias methods is the underestimation of the diffuse reflectance beyond the start of the target depth range in which the diffusive scattering is biased. The application of the first bias in the backward direction implies that there is a smaller probability that the photon packets will be forward scattered beyond the early part of the target depth range. This causes a systematic statistical bias that also affects the performance of the angle biasing procedure used in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

] and also the bias procedure used in [10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

], which limits their application to a thin target layer. For these reasons, these bias procedures are limited to model systems in which only forward bias scattering is applied.

We ensure the correct statistics in the importance sampling method by splitting the photon packet that is biased towards the apparent position of the collecting optics in two photon packets prior to the first biased scattering [19

19. I. T. Lima, Jr., “Advanced Monte Carlo methods applied to Optical Coherence Tomography” (invited), presented at the 2009 SBMO/IEEE MTT-S International Microwave and Optoelectronics Conference, Belém, Brazil, 3–6 Nov. 2009.

]. One of these two photon packets is the one biased towards the collecting optical system, which is associated with the likelihood ratio specified in Eq. (6), as described earlier in this section, and we apply the successive biased scatterings until the photon packet is terminated. Then the other photon packet continues a forward propagation starting at the location of the biased backward scattering point, having the initial direction u^ and scattered by the standard procedure that is described in Sec. 2. We ensure that there is no systematic statistical bias in the statistical result of the forward propagating split photon packet by assigning to this second photon the likelihood ratio L'(i), which is the complement of the likelihood ratio of the biased backward scattered photon packet L(i), so that L'(i)=1L(i). This second photon packet, which is only created if L(i)<1, is also allowed to undergo one biased backscattering, which may lead to another photon packet split, and successive additional biased scatterings towards the tip of the collecting optical system until the photon packet propagates beyond the simulation domain. In the cases that we studied, this procedure increased the computational time of each sample by a factor of 5 when compared with a sample computed using the standard Monte Carlo. The increase of computational time of importance sampling compared with the standard method depends on the average value of the mean free path and on the width of the target depth range. It is important to point out that each split photon packet is not counted as an additional photon packet when determining the value of N in Eqs. (9) and (10), since the use of the likelihood ratio associated to each photon packet in these equations will produce the correct statistical result. Once a photon packet exceeds the region within the depth target layer, it will no longer be biased and will likely be terminated after exceeding the boundary of the last layer while propagating in the forward direction. Then a new photon packet will be created at the origin as in the standard MCML method, and a new Monte Carlo sample will be simulated. Despite the higher computational cost per photon packet, the computation of the Class I diffuse reflectance in the case that we studied using Monte Carlo simulations with importance sampling required as little as one-thousandth of the computational time required to achieve the same accuracy in the diffuse reflectance calculation using the standard Monte Carlo method.

4. Validation of our new importance sampling method

We validate the calculation of the Class I diffuse reflectance using Monte Carlo simulations with the importance sampling method that we developed by comparing it against extensive standard Monte Carlo simulations. We simulate a turbid media that consists of three layers, shown schematically in Fig. 2
Fig. 2 Schematic representation of a simulation setup similar to [5].
. The first layer extends from 0 µm to 330 µm, and the third layer extends from 360 µm to 1.2 mm, have an absorption coefficient µa = 1.5 cm−1 and a scattering coefficient µs = 60 cm−1. The second layer extends from 330 µm to 360 µm and has an absorption coefficient µa = 3 cm−1 and a scattering coefficient µs = 120 cm−1. We assume the three layers to have the same refractive index n = 1 and an anisotropy factor g = 0.9, as in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. We simulate a TD-OCT system that is delivered and collected by the tip of an optical fiber with a radius of 10 µm and an acceptance angle of 5 degrees as in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. For simplicity, a point source of light is assumed along the vertical direction as in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

,10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

], since the objective of this study is to validate and demonstrate the effectiveness of the importance sampling implementation that we developed.

In Figs. 3
Fig. 3 The Class I diffuse reflectance as a function of the depth for the first simulation setup described in Sec. 4, whose schematic representation is shown in Fig. 2. The black solid line are results obtained with 2 × 106 Monte Carlo simulations with the importance sampling method described in Sec. 4. The green long-dashed line are results obtained with 1010 standard Monte Carlo samples. The black double-dashed line are results obtained using the angle bias procedure described in [5]. The pink short-dashed lines are results of plus and minus one standard deviation of Monte Carlo simulations with importance sampling with 2 × 106 samples that were estimated using an ensemble with 64 × 106 simulations.
and 4
Fig. 4 The same Class I diffuse reflectance results shown in Fig. 3, except that the Class I diffuse reflectance is shown in linear scale for the depth interval from 310 µm to 380 µm. The error bars shown for the solid and for the dashed line were estimated in each respective Monte Carlo simulation.
, we show results obtained with 2 × 106 Monte Carlo simulations with importance sampling, which has a computational cost slightly smaller than 107 standard Monte Carlo simulations. In Figs. 3 and 4, we also show results obtained with the angle biasing method described in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

] when applied to the entire range from 0 µm to 600 µm. Because of the systematic bias due to the artificial reduction in the number of photons penetrating the tissue as pointed out in Sec. 3.4, this method can only be applied to calculate the diffusive reflectance produced by a very narrow layer as described in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

]. We run the Monte Carlo simulations with importance sampling with the bias coefficient a = 0.9. The choice of a = g is the value for the bias coefficient that enabled the fastest conversion of the statistical results in (9) and (10) with respect to the number of simulated photon packets because it produces the best combination of strong bias with a limited range of variation of the likelihood ratio in Eq. (8). Therefore, the results shown in Figs. 3 and 4 indicate that our new importance sampling procedure reduces the computational cost to obtain the Class I diffuse reflectance by about three orders of magnitude.

In Fig. 5
Fig. 5 The number of Class I photon packets diffusely reflected as a function of the depth in the simulations whose results are shown in Fig. 3. The black solid line are results obtained with 2 × 106 Monte Carlo simulations with the importance sampling method described in Sec. 3. The green long-dashed line are results obtained with 1010 standard Monte Carlo simulations. The pink short-dashed line are results obtained with 107 previously implemented Monte Carlo method samples.
we show a clear evidence of the effectiveness of our importance sampling procedure that we described in Sec. 3. While an ensemble with 107 standard Monte Carlo simulations produces a very small number of Class I diffusely reflected photon packets, an ensemble with 2 × 106 Monte Carlo simulations with importance sampling produces a number of Class I diffusely reflected photon packets that exceeds the number of Class I diffusely reflected photon packets produced by as many as 1010 standard Monte Carlo simulations at depths greater than 132 µm. The agreement between the Class I diffuse reflectance, obtained with our importance sampling and with the standard Monte Carlo method with a much larger number of samples, indicates that the biased samples with importance sampling are weighted in a way that eliminates any residual statistical bias. Therefore, not only does our importance sampling method converge to the true statistical result, but it does so at a much faster rate than the standard Monte Carlo method.

We also validate our method to simulate an TD-OCT imaging system in the presence multiple layers with different scattering properties and refractive indices by comparing it against extensive standard Monte Carlo simulations. In this case, the apparent position of the detector is calculated recursively in all the tissue layers using Eqs. (5) and (6). The light is emitted by a fiber optic probe that is radially reflected by a prism, as shown in Fig. 6
Fig. 6 Schematic representation of a simulation setup with multiple layers with different refractive indices.
.

The optical system has a focusing lens that has a numerical aperture that enables collecting light at an angle of up to 4 degrees and a diameter of 0.5 mm. As in [5

5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

,10

10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

,11

11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

], a point source that emits in the vertical direction is assumed. There is air between the center of the probe and the first layer, which is located at 2.12 mm from the center of the fiber. We simulate a turbid media that consists of three layers with refraction occurring at their interfaces. The first diffusive layer, which extends from 2.12 mm to 2.22 mm from the center of the fiber, has the following specifications: refractive index n = 1, absorption coefficient µa = 1.5 cm−1, and the scattering coefficient µs = 60 cm−1. The second diffusive layer, which is extends from 2.32 mm to 2.42 mm from the center of the fiber, has the same absorption and scattering coefficient as the first layer, but its refractive index is n = 1.33. The third layer, which extends from 2.42 mm to 2.62 mm from the center of the fiber, has the following specifications: refractive index n = 1, absorption coefficient µa = 1.5 cm−1, and the scattering coefficient µs = 30 cm−1. We assume there is air at the end of the third layer. We consider the three diffusive layers have anisotropy factor g = 0.9. In Fig. 7
Fig. 7 The Class I diffuse reflectance as a function of the distance from the center of the optical fiber for the second simulation setup described in Sec. 4, whose schematic representation is shown in Fig. 6. The black solid line are results obtained with 2 × 105 Monte Carlo simulations with the importance sampling method described in Sec. 4. The green long-dashed line are results obtained with 109 standard Monte Carlo samples. The blue dots are results obtained with 106 standard Monte Carlo simulations. The pink short-dashed lines are results of plus and minus one standard deviation of Monte Carlo simulations with importance sampling with 2 × 105 samples that were estimated using an ensemble with 64 × 105 simulations.
we show the simulation results of this TD-OCT setup probing a tissue with multiple layers with different refractive indices. We observed an excellent convergence between the results obtained with our new importance sampling method and the results obtained using standard Monte Carlo simulations using MCML. Our results, however, were obtained in one-thousandth of the time required by standard method.

5. Conclusion

We developed and validated a new importance sampling method that reduces the computation time of obtaining Class I diffuse reflectance in TD-OCT simulations using Monte Carlo simulations by three orders of magnitude. Practically, this amounts to reducing the computation time of the diffuse reflectance from several hours, which may be prohibitive for many practical applications, to seconds. We also showed how our method extends and improves existing methods that were proposed to speed up Monte Carlo simulations of general light transport in tissue. We are currently developing another new importance sampling method to obtain Class II diffuse reflectance, which limits the imaging depth in OCT systems, in a computationally efficient way. Our fast Monte Carlo calculation of TD-OCT signals could be used to further study the physical process governing both Class I and Class II signals, e.g., obtain the statistics of the depth scan, including the effects of multiple scattering of light, in TD-OCT. This paper focused on TD-OCT, but we are currently generalizing our results to FD-OCT. This is an important prerequisite to many efforts to increase penetration depth and to better image extraction in OCT systems.

References and links

1.

A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]

2.

C. Xu, C. Vinegoni, T. S. Ralston, W. Luo, W. Tan, and S. A. Boppart, “Spectroscopic spectral-domain optical coherence microscopy,” Opt. Lett. 31(8), 1079–1081 (2006). [CrossRef] [PubMed]

3.

M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef] [PubMed]

4.

B. Liu and M. E. Brezinski, “Theoretical and practical considerations on detection performance of time domain, Fourier domain, and swept source optical coherence tomography,” J. Biomed. Opt. 12(4), 044007 (2007). [CrossRef] [PubMed]

5.

G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]

6.

S. S. Sherif, C. C. Rosa, C. Flueraru, S. Chang, Y. Mao, and A. G. Podoleanu, “Statistics of the depth-scan photocurrent in time-domain optical coherence tomography,” J. Opt. Soc. Am. A 25(1), 16–20 (2008). [CrossRef] [PubMed]

7.

M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple scattering in optical coherence microscopy,” Appl. Opt. 34(25), 5699–5707 (1995). [CrossRef] [PubMed]

8.

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

9.

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

10.

N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]

11.

N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

12.

R. Y. Rubinstein, Simulation and the Monte Carlo Method (Wiley, 1981).

13.

G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(3), 310–312 (2002). [CrossRef]

14.

S. L. Fogal, G. Biondini, and W. L. Kath, “Multiple importance sampling for first- and second-order polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(9), 1273–1275 (2002). [CrossRef]

15.

I. T. Lima Jr, A. O. Lima, J. Zweck, and C. R. Menyuk, “Efficient computation of outage probabilities due to polarization effects in a WDM system using a reduced Stokes model and importance sampling,” IEEE Photon. Technol. Lett. 15(1), 45–47 (2003). [CrossRef]

16.

I. T. Lima, A. O. Lima, G. Biondini, C. R. Menyuk, and W. L. Kath, “A comparative study of single section polarization-mode dispersion compensators,” J. Lightwave Technol. 22(4), 1023–1032 (2004). [CrossRef]

17.

J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A 13(5), 952–961 (1996). [CrossRef] [PubMed]

18.

H. Iwabuchi, ““Efficient Monte Carlo method for radiative transfer modeling,” J. of the Atmosph,” Science 63, 2324–2339 (2006).

19.

I. T. Lima, Jr., “Advanced Monte Carlo methods applied to Optical Coherence Tomography” (invited), presented at the 2009 SBMO/IEEE MTT-S International Microwave and Optoelectronics Conference, Belém, Brazil, 3–6 Nov. 2009.

20.

“Monte Carlo simulations,” Oregon Medical Laser Center, accessed January 1, 2009, http://omlc.ogi.edu/software/mc/

21.

S. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory (Prentice-Hall, 1993).

22.

E. Hecht, Optics, 4th ed. (Pearson Addison Wesley, 2003).

OCIS Codes
(110.4500) Imaging systems : Optical coherence tomography
(170.3660) Medical optics and biotechnology : Light propagation in tissues

ToC Category:
Optics of Tissue and Turbid Media

History
Original Manuscript: February 3, 2011
Revised Manuscript: March 31, 2011
Manuscript Accepted: March 31, 2011
Published: April 4, 2011

Citation
Ivan T. Lima, Anshul Kalra, and Sherif S. Sherif, "Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography," Biomed. Opt. Express 2, 1069-1081 (2011)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-2-5-1069


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. 66(2), 239–303 (2003). [CrossRef]
  2. C. Xu, C. Vinegoni, T. S. Ralston, W. Luo, W. Tan, and S. A. Boppart, “Spectroscopic spectral-domain optical coherence microscopy,” Opt. Lett. 31(8), 1079–1081 (2006). [CrossRef] [PubMed]
  3. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef] [PubMed]
  4. B. Liu and M. E. Brezinski, “Theoretical and practical considerations on detection performance of time domain, Fourier domain, and swept source optical coherence tomography,” J. Biomed. Opt. 12(4), 044007 (2007). [CrossRef] [PubMed]
  5. G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef] [PubMed]
  6. S. S. Sherif, C. C. Rosa, C. Flueraru, S. Chang, Y. Mao, and A. G. Podoleanu, “Statistics of the depth-scan photocurrent in time-domain optical coherence tomography,” J. Opt. Soc. Am. A 25(1), 16–20 (2008). [CrossRef] [PubMed]
  7. M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple scattering in optical coherence microscopy,” Appl. Opt. 34(25), 5699–5707 (1995). [CrossRef] [PubMed]
  8. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). [CrossRef] [PubMed]
  9. L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef] [PubMed]
  10. N. G. Chen and J. Bai, “Estimation of quasi-straightforward propagating light in tissues,” Phys. Med. Biol. 44(7), 1669–1676 (1999). [CrossRef] [PubMed]
  11. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]
  12. R. Y. Rubinstein, Simulation and the Monte Carlo Method (Wiley, 1981).
  13. G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(3), 310–312 (2002). [CrossRef]
  14. S. L. Fogal, G. Biondini, and W. L. Kath, “Multiple importance sampling for first- and second-order polarization-mode dispersion,” IEEE Photon. Technol. Lett. 14(9), 1273–1275 (2002). [CrossRef]
  15. I. T. Lima, A. O. Lima, J. Zweck, and C. R. Menyuk, “Efficient computation of outage probabilities due to polarization effects in a WDM system using a reduced Stokes model and importance sampling,” IEEE Photon. Technol. Lett. 15(1), 45–47 (2003). [CrossRef]
  16. I. T. Lima, A. O. Lima, G. Biondini, C. R. Menyuk, and W. L. Kath, “A comparative study of single section polarization-mode dispersion compensators,” J. Lightwave Technol. 22(4), 1023–1032 (2004). [CrossRef]
  17. J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A 13(5), 952–961 (1996). [CrossRef] [PubMed]
  18. H. Iwabuchi, ““Efficient Monte Carlo method for radiative transfer modeling,” J. of the Atmosph,” Science 63, 2324–2339 (2006).
  19. I. T. Lima, Jr., “Advanced Monte Carlo methods applied to Optical Coherence Tomography” (invited), presented at the 2009 SBMO/IEEE MTT-S International Microwave and Optoelectronics Conference, Belém, Brazil, 3–6 Nov. 2009.
  20. “Monte Carlo simulations,” Oregon Medical Laser Center, accessed January 1, 2009, http://omlc.ogi.edu/software/mc/
  21. S. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory (Prentice-Hall, 1993).
  22. E. Hecht, Optics, 4th ed. (Pearson Addison Wesley, 2003).

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.

CrossCheck Deposited