## Fast calculation of multipath diffusive reflectance in optical coherence tomography |

Biomedical Optics Express, Vol. 3, Issue 4, pp. 692-700 (2012)

http://dx.doi.org/10.1364/BOE.3.000692

Acrobat PDF (1192 KB)

### Abstract

We show how to efficiently calculate the signal in optical coherence tomography (OCT) systems due to the ballistic photons, the quasi-ballistic photons, and the photons that undergo multiple diffusive scattering using Monte Carlo simulations with importance sampling. This method enables the calculation of these three components of the OCT signal with less than one hundredth of the computational time required by the conventional Monte Carlo method. Therefore, it can be used as a design tool to characterize the performance of OCT systems, and can also be used in the development of novel signal processing techniques that can extend the imaging range of OCT systems. We investigate the parameter dependence of our importance sampling method and we validate it by comparison to an existing method.

© 2012 OSA

## 1. Introduction

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]

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

4. I. T. Lima Jr, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express **2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

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

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

7. I. T. Lima Jr, A. M. Oliveira, 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]

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

9. H. Iwabuchi, “Efficient Monte Carlo method for radiative transfer modeling,” J. Atmos. Sci. **63**(9), 2324–2339 (2006). [CrossRef]

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

4. I. T. Lima Jr, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express **2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

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

4. I. T. Lima Jr, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express **2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

*p*specified for the ensemble of Monte Carlo simulations [12]. In the event in which a bias towards the collecting optics is not applied, the scattering will follow the unbiased model. Moreover, the biased scatterings in this method can only produce scatterings in the half-sphere along the apparent direction of the collecting optics, as opposed to scattering in all directions as in the previous method. These procedures enable a significant increase in the number of collected photon that undergo multiple diffusive scattering, along with a decrease in the range of values of the likelihood ratio that accelerate the convergence of the calculation of the Class II reflectance. The combination of these two factors leads to faster convergence of the Monte Carlo simulations for the calculation of both the Class I and the Class II reflectances.

## 2. Importance sampling applied to OCT

13. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. **1**(4), 153–156 (1969). [CrossRef]

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

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

*g*of the tissue, in the unbiased scatterings we use the Henyey-Greenstein probability density function that is defined as

*θ*is the angle between the photon packet propagation direction

_{s}*ϕ*that is randomly picked from a uniform probability density function from 0 to 2π.

### 2.1. First biased scattering

*g*of tissue is close to 1, and the probability of reflection also decreases with the depth, the convergence of the calculation of the Class I and the Class II reflectances using standard Monte Carlo simulations is very slow, as shown in [3

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

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

*a*is a bias coefficient in the range (0,1). After randomly picking a biased 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 the MCML software package to enable a full three-dimensional scattering. Then, the scattered photon packet is associated with a quantity that is defined as the likelihood ratio [5

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

7. I. T. Lima Jr, A. M. Oliveira, 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]

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

### 2.2. Additional biased scatterings

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

*p*≤ 1, as opposed to

*p*= 1 as it is the previous method. If a biased scattering is not applied at a given point in the tissue, the photon packet will undergo an unbiased scattering at that location according to the probability density function in Eq. (1). The likelihood ratio of this scattering, whether the biased or the unbiased probability density function is randomly selected, is calculated according to the expression

*p*,

*p*, the unbiased function

### 2.3. Calculation of the reflectance

*z*are obtained by calculating the mean value of the indicator functions

*I*

_{1}and

*I*

_{2}of the spatial and temporal filter of the Class I reflectance and the Class II reflectance, respectively, for all the photon packets (samples) in the ensemble. The indicator function

*I*

_{1}and

*I*

_{2}of the spatial and temporal filter for the

*i*th photon packet is defined as

*z*is the depth imaged,

*l*is the coherence length of the source,

_{c}*r*is the distance of the

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

*θ*is the angle with the

_{z,i}*z*-axis, which is the axis normal to the tissue interface, Δ

*s*is the optical path,

_{i}*z*

_{max}is the maximum depth reached by the photon packet. A detected photon packet is considered a Class II photon packet if it does not reach a depth that is consistent with its optical path, so that it interferes constructively with corresponding detected Class I photons packets without bringing any information from the depth in which those Class I photons packets were reflected. For simplicity, the indicator functions in Eq. (5) and in Eq. (6) were defined with a square time gating as in [3

**44**(9), 2307–2320 (1999). [CrossRef] [PubMed]

*W*(

*i*) is the weight of the

*i*th photon packet in MCML, which is a quantity affected by the absorption coefficient at the scattering points, and L(

*i*) is the product of the likelihood rations of all the biased scatterings that affected the

*i*th photon packet. Using the Monte Carlo method with the importance sampling implementation described in this section, the calculation of the reflectances in Eq. (7) converge two to three orders of magnitude faster with the number of samples

*N*than the standard Monte Carlo method used in MCML.

### 2.4. Generation of random biased angles

16. The Gnu Project, “Gnu Scientific Library,” accessed June 15, 2011, http://www.gnu.org/s/gsl/.

*u*is sampled from the random number generator of the Gnu Scientific Library. Equation (9) was derived based on the probability theory in [17].

_{i}## 3. Numerical results

*µ*= 1.5 cm

_{a}^{−1}and a scattering coefficient

*µ*= 60 cm

_{s}^{−1}, but also contains five thin layers with absorption coefficient

*µ*= 3 cm

_{a}^{−1}and a scattering coefficient

*µ*= 120 cm

_{s}^{−1}. These five thin layers with higher scattering coefficient are located from 200 µm to 215 µm, from 365 µm to 395 µm, from 645 µm to 660 µm, from 760 µm to 775 µm, and from 900 µm to 915 µm. We assume that this tissue has the same refractive index

*n*= 1 and an anisotropy factor

*g*= 0.9, as in [3

**44**(9), 2307–2320 (1999). [CrossRef] [PubMed]

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

**44**(9), 2307–2320 (1999). [CrossRef] [PubMed]

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

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

^{8}Monte Carlo simulations with importance sampling, which has a computational cost of about 9 × 10

^{8}standard Monte Carlo simulations in this particular simulation setup because of the computational cost of the photon splitting procedure [4

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

*a*= 0.925 and additional bias probability

*p*= 0.5. The results shown in Figs. 2 and 3 indicate that our new importance sampling procedure reduces the computational cost to obtain the Class I diffuse reflectance by about three orders of magnitude when compared to the standard Monte Carlo method.

*a*for

*p*= 0.5. The depths at 400 µm and 670 µm correspond to the tissue regions near the second and the third peaks of the reflectance. The relative error is defined as the ratio between the standard deviation, which is the square root of the variance in Eq. (8), divided by the estimated value of the reflectance in Eq. (7). The Class I reflectance has its minimum relative error at 400 µm close to

*a*= 0.925, but the minimum shifts to about

*a*= 0.95 µm at 670 µm because a stronger bias is necessary to collect more samples from deeper regions in the tissue. However, as the bias coefficient is increased towards 1, larger variations in the likelihood ratio due to very strong bias leads to an increase in the relative error with the bias coefficient. The Class II reflectance has its minimum relative error at 400 µm close to

*a*= 0.91, and shifts to only about

*a*= 0.925 µm at 670 µm. That is lower than the optimum bias coefficient observed in the Class I reflectance because strong bias leads to an increase in the number of ballistic and quasi-ballistic photons and a decrease in the number of collected photons that were scattered multiple times in the tissue. Figure 4 shows that there is a region between 0.9 and 0.95 for the bias parameter

*a*that enables a rapid convergence of the calculation of the Class I and Class II reflectance with the Monte Carlo method with this importance sampling implementation.

*p*for

*a*= 0.9. At the depth of 400 µm, the relative error is minimized when

*p*= 0.5 for both the Class I and the Class II reflectances. At the depth of 670 µm the relative error of the Class I reflectance is minimized at

*p*= 0.6, while the Class II reflectance has its relative error minimized at

*p*= 0.55. Nevertheless, the relative error of the numerical results is not very sensitive to the value of

*p*when

*p*is between 0.3 and 0.7. At

*p*= 1 the bias used is comparable to the one presented in [4

**2**(5), 1069–1081 (2011). [CrossRef] [PubMed]

*p*= 0.55, which implies into a convergence in the new method that is four times faster than the previous method. At

*p*= 0, which corresponds to the parameter regime in which only one bias scattering is applied, this method is the least effective in the reflectance calculation.

## 4. Conclusion

## Acknowledgments

## References and links

1. | A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. |

2. | C. Xu, C. Vinegoni, T. S. Ralston, W. Luo, W. Tan, and S. A. Boppart, “Spectroscopic spectral-domain optical coherence microscopy,” Opt. Lett. |

3. | G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. |

4. | I. T. Lima Jr, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express |

5. | G. Biondini, W. L. Kath, and C. R. Menyuk, “Importance sampling for polarization-mode dispersion,” IEEE Photon. Technol. Lett. |

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

7. | I. T. Lima Jr, A. M. Oliveira, G. Biondini, C. R. Menyuk, and W. L. Kath, “A comparative study of single section polarization-mode dispersion compensators,” J. Lightwave Technol. |

8. | J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A |

9. | H. Iwabuchi, “Efficient Monte Carlo method for radiative transfer modeling,” J. Atmos. Sci. |

10. | N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. |

11. | M. J. Yadlowsky, J. M. Schmitt, and R. F. Bonner, “Multiple scattering in optical coherence microscopy,” Appl. Opt. |

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

13. | E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. |

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

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

16. | The Gnu Project, “Gnu Scientific Library,” accessed June 15, 2011, http://www.gnu.org/s/gsl/. |

17. | A. Papoulis, |

**OCIS Codes**

(110.4500) Imaging systems : Optical coherence tomography

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

**ToC Category:**

Optical Coherence Tomography

**History**

Original Manuscript: January 3, 2012

Revised Manuscript: February 10, 2012

Manuscript Accepted: February 10, 2012

Published: March 12, 2012

**Citation**

Ivan T. Lima, Anshul Kalra, Hugo E. Hernández-Figueroa, and Sherif S. Sherif, "Fast calculation of multipath diffusive reflectance in optical coherence tomography," Biomed. Opt. Express **3**, 692-700 (2012)

http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-3-4-692

Sort: Year | Journal | Reset

### References

- 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]
- 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]
- 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]
- I. T. Lima, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express2(5), 1069–1081 (2011). [CrossRef] [PubMed]
- 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]
- 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]
- I. T. Lima, A. M. Oliveira, 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]
- J. M. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A13(5), 952–961 (1996). [CrossRef] [PubMed]
- H. Iwabuchi, “Efficient Monte Carlo method for radiative transfer modeling,” J. Atmos. Sci.63(9), 2324–2339 (2006). [CrossRef]
- N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt.46(10), 1597–1603 (2007). [CrossRef] [PubMed]
- 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]
- 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.
- E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun.1(4), 153–156 (1969). [CrossRef]
- 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]
- “Monte Carlo simulations,” Oregon Medical Laser Center, accessed Jan. 1, 2009, http://omlc.ogi.edu/software/mc/ .
- The Gnu Project, “Gnu Scientific Library,” accessed June 15, 2011, http://www.gnu.org/s/gsl/ .
- A. Papoulis, Probability, Random Variables, and Stochastic Processes (McGraw-Hill, New York, 1984).

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