Abstract
We describe two importance sampling techniques for a standard Monte Carlo (MC) method that could enable fast simulation of signals from optical coherence tomography (OCT) imaging systems. These OCT signals are generated due to diffusive reflections from either multilayered or arbitrary shaped, turbid media, for example, tissue. Such signals typically consist of ballistic and quasi-ballistic components, of scattered photons inside the medium, in addition to photons that undergo multiple scattering. We show that MC simulation of these OCT signals using importance sampling reduces its computation time on a serial processor by up to three orders of magnitude compared to its corresponding standard implementation. Therefore, these importance sampling techniques enable practical simulation of OCT B-scans of turbid media, for example, tissue, using commonly available workstations.
Keywords
- optical coherence tomography
- Monte Carlo simulation
- light transport in turbid media
- importance sampling
1. Introduction
Optical coherence tomography (OCT) is a non-invasive sub-surface imaging technique that has experienced significant growth in biomedical applications [1, 2]. OCT systems could be implemented with a low-coherence light source and a mechanical scanning sub-system (time-domain OCT). More advanced systems use a low-coherence light source with a spectrometer or a wavelength-swept source (frequency-domain OCT). OCT has an imaging depth that could reach up to 3 mm, depending on the optical properties of the tissue, and it also has one to two orders of magnitude higher resolution than ultrasound imaging. OCT could also produce images inside the body when it is implemented using optical fiber probes. Unlike X-ray or gamma-ray imaging, OCT is safe for biological tissues because it utilizes non-ionizing radiation mainly in the infrared spectrum.
1.1 OCT signal simulation using a Monte Carlo method
The signal obtained by an OCT imaging system consists of ballistic and quasi-ballistic (Class I diffuse reflectance) photons, in addition to multiply scattered photons (Class II diffuse reflectance), that are reflected from tissue [3]. However, multiply scattered photons do not carry practically useful information about the imaged tissue; therefore, they result in a degradation of the OCT signal [4]. In addition, it has been shown that Class II diffuse reflectance represents a fundamental limit related to the imaging depth of OCT in tissue [5]. Therefore, it is important to account for both Class I and Class II diffusive reflectance in any practical simulation of OCT signals.
Since it is not practical to simulate light transport in turbid media, for example, tissue, using electromagnetic waves, especially due to diffusive scattering, a Monte Carlo (MC) method of simulating light transport in tissue has been typically used [6, 7, 8, 9]. However, the computational cost of this MC-based simulation of OCT systems could be very high, as the probability of detecting diffusively reflected photons from tissue is very low [4, 5].
To reduce the computational cost, thereby accelerate, this MC simulation, importance sampling could be used to speed up simulations by orders of magnitude. Importance sampling has been applied earlier to optical communications [10, 11], confocal microscopy [12], atmospheric optics [13], and diffuse optical tomography [8].
To improve the computational efficiency of the MC-based simulation of OCT systems [6], Yao and Wang proposed the first importance sampling technique to simulate OCT signals from a multilayered turbid medium [3]. However, their method only enabled the simulation of OCT signals from a thin shallow layer in tissue, as the results obtained from deeper tissue regions were underestimated. In [14], we, the authors of this chapter, developed another more advanced importance sampling technique by implementing multiple biased scatterings per photon packet, and by developing a photon splitting procedure. Our advanced importance sampling resulted in more accurate and computationally efficient simulations of OCT signals due to ballistic and quasi-ballistic photons. However, it still underestimated OCT signals due to multiply scattered photons. To enable accurate simulation of OCT signals due to both Class I and Class II diffusive reflectances, we further developed our importance sampling technique to accurately and efficiently simulate diffusive reflectance due to photons that undergo multiple scattering [15]. In this method, additional biased scatterings were randomly applied, which enabled accurate simulation of both Class I and Class II diffusive reflectances, with a speed-up of three orders of magnitude compared to the standard MC method.
Our advanced importance sampling techniques above were implemented to simulate OCT of tissues with planar geometry [6]. To enable simulation of OCT of arbitrarily shaped turbid media, we used the mesh-based MC method of light transport in tissue proposed by Fang [16]. This method uses a Plücker coordinate system to efficiently calculate intersections between paths of light propagation with interfaces of the object regions that are modeled using tetrahedrons. We combined this mesh-based MC method with our importance sampling techniques to simulate OCT signals from tissue with arbitrarily shaped regions. However, since it was still computationally costly to simulate a full OCT B-scan using this method, we also developed a parallel implementation of this simulator that exploited the massively parallel computing capabilities of Graphics Processing Units (GPUs) to accelerate this simulator by two additional orders of magnitude [17, 18]. This GPU-based implementation enabled simulation of OCT B-scans of arbitrarily shaped turbid media in a few minutes using commonly available workstations.
In Section 2 of this chapter, we present a standard MC method for simulating OCT signals. In Section 3, we present our first importance sampling implementation that enables the simulation of OCT signals from higher depths inside turbid media. In Section 4, we present our more advanced importance sampling implementation that accurately calculates both Class I and Class II diffusive reflectances, and is three orders of magnitude faster than the standard MC simulation.
2. Standard MC method for simulating OCT signals
Our implementation of the MC method to simulate OCT signals is based on Monte Carlo simulation of light transport in multilayered tissues (MCML) [6]. MCML simulates an ensemble of photon packets that are launched as a steady-state pencil beam, normal to the top surface of the medium. Within the tissue, each such photon packet undergoes a random walk whose step size is determined by an exponentially distributed random variable that is parameterized by the
where
The Class I diffuse reflectance at depth equal to
where
At any depth, the diffuse reflectance
and
where
3. Importance sampling for simulation of Class I OCT signal
Our first importance sampling technique to simulate OCT signals aimed at increasing the number of photons collected at the detector. This algorithm uses the same method described in Section 2, where we also use the same square time gating given by [3].
Since most tissues are highly forward-scattering, their anisotropy factor is close to 1. Therefore, there is a very small probability that a simulated photon packet at any given depth in the tissue would undergo scattering in the backward direction toward the OCT probe. The probability of collecting Class I photons drops rapidly with depth in the tissue from which the photon is scattered in the backward direction. To allow faster simulation of Class I photons, we designed an importance sampling method that biases the direction
All photon packets propagating in a direction close to
3.1 Scattering angle due to first event of backscattering
As the photon packet reaches the depth range targeted, the propagation direction
where
where cos(
3.2 Scattering angles of additional events of backscattering
As a photon packet is biased toward the apparent position of the collecting optical fiber, at any given depth in the tissue, the photon packet becomes more likely to be collected at the tip of the fiber. However, the photon packet could be scattered several times after the first backscatter bias before reaching the optical collection system. These additional scatterings, according to Eq. (1), reduce the correlation between the biased direction and the event in which the photon packet is collected. We overcome this reduction in correlation by continuing to bias the scattering direction
Once a photon packet experiences the first biased scattering, that photon packet is biased at all additional scattering points until it is removed from the simulation, which can occur when the photon packet is removed by Russian roulette, as described in Section 2, or it leaves the tissue. After simulating
and
Eqs. (7) and (8) are similar to Eqs. (3) and (4), except that the indicator function is multiplied by its corresponding likelihood ratio. Using this method, a significantly larger number of photon packets are scattered from a specific depth range toward the collecting optical system than the number obtained using a standard MCML implementation. At the end of this biased simulation, each photon packet is weighted by its likelihood ratio, which adjusts the contribution of each packet to the estimation of the Class I diffuse reflectance. The estimated diffuse reflectance converges toward its true value faster, by several orders of magnitude, when compared to the standard Monte Carlo method.
3.3 Importance sampling effectiveness and depth of tissue
One drawback of previously existing bias methods, for example, [3, 7, 8] is an underestimation of the diffuse reflectance beyond the targeted depth range. The application of the first backward bias reduces the probability that this photon packet would propagate beyond that portion of the tissue. This would lead to a statistical bias to this importance sampling method similar to that in the angle biasing procedure used in [3] and the method used in [7, 8], which limits the effectiveness of those methods to a thin target layer.
We make sure we obtain correct statistics by splitting the photon packet into two photon packets before the first biased scattering [14]. The first of these two photon packets is the one biased toward the collecting optical system. The second photon packet starts propagating from the location in which the biased backscattering occurred, where its initial direction calculated by applying the standard procedure to the previous direction
3.4 Numerical results
We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. As shown in Figure 2, light is emitted by an optical fiber probe that is reflected by a prism.
The shown optical system has a focusing lens with a numerical aperture (NA) that allows collecting light at an angle of up to 4° and a diameter of 0.5 mm. Similar to the setups in [5, 10, 11], we assume a point source that emits in the vertical direction. Air is present between the center of the probe and the first layer of tissue, which is placed 2.12 mm from the center of the fiber. We simulate a three-layer turbid medium with refractive-index mismatch at its interfaces. The first layer, extending from 2.12 to 2.22 mm from the tip of the fiber, has absorption coefficient
From Figure 3, we note an excellent correspondence between results obtained with our new importance sampling method and results obtained using MCML, that is, standard Monte Carlo simulations. However, our results were obtained in one-thousandth of the time required by the standard method.
4. Importance sampling for simulation of Class I and Class II OCT signals
In this section, we further improve the importance sampling technique that was described in Section 3, so we can simulate Class II OCT signals more accurately and more efficiently [19].
4.1 Scattering angle of first backscattering event
In the MC simulation described in Section 3, we note that the bias function in (5) produces large values of the likelihood ratio (>>1) when photon packets are scattered in the then unlikely forward direction. These photon packets contributed to a slow decrease in the relative variation, which corresponds to relative error, with the increase in the number of photon packets launched for the calculation of the OCT signal. Referring to Figure 1, we could reduce this relative variation by choosing a distribution function for the scattering angle that limits it to the backward direction. This modified distribution is given by:
where
where
To sample angles according to the biased probability density function in (8), one could use any uniform pseudo-random number generator that would be typically available in scientific software libraries. For example, if
This conversion formula was derived using probability theory [20].
4.2 Scattering angles of additional biased backscatterings
A second enhancement that could be made to the importance sampling technique, described in Section 3, is to bias the additional scatterings toward the center of the OCT probe
If the biased function
4.3 Numerical results
We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. We consider a tissue model that consists of multiple layers that could be imaged with an OCT system, as shown schematically in Figure 4. The modeled tissue extends from 0 to 1 mm, and consists primarily of a turbid layer with absorption coefficient
In Figures 5 and 6, we show results obtained with 108 Monte Carlo photon packets with importance sampling, which has a computational cost of simulating about 9 × 108 photon packets using standard Monte Carlo. The computational cost of applying this importance sampling technique depends on the target depth range, and on the average photon mean free path in the given tissue. The target depths in the shown simulations were set from 0 to 1 mm. Therefore, every single photon scattering that occurs in the depth range from 0 to 1 mm would be biased. We used a bias coefficient
We note that the results obtained with the MCML have confidence intervals that are noticeably larger than those of the results obtained with importance sampling shown in Figure 6, even though the standard Monte Carlo simulations have a computational cost 113 times larger than those obtained with importance sampling. In Figure 6, we also note that our importance sampling technique reduces the computational cost of calculating the Class II reflectance by more than two orders of magnitude.
In Figure 7 we show the relationship between the relative error in calculating Class I and the Class II reflectances at two different depths: 400 and 670 μm and the bias coefficient
We note that Class I reflectance has a minimum relative error at 400 μm with
5. Conclusions
We described two importance sampling techniques for a standard Monte Carlo (MC) method that could enable fast simulation of signals from optical coherence tomography (OCT) imaging systems. These OCT signals are generated due to diffusive reflections from either multilayered or arbitrarily shaped, turbid media, for example, tissue. Such signals typically consist of ballistic and quasi-ballistic components, of scattered photons inside the medium, in addition to photons that undergo multiple scattering. We showed that MC simulation of these OCT signals using our importance sampling reduced its computation time on a serial processor by up to three orders of magnitude compared to its corresponding standard implementation. Therefore, our importance sampling techniques enable practical simulation of OCT B-scans of turbid media, for example, tissue, using commonly available workstations.
References
- 1.
Brezinski M, Fujimoto J. Optical coherence tomography: High-resolution imaging in non transparent tissue. IEEE Journal of Selected Topics in Quantum Electronics. 1999; 5 :1185-1192 - 2.
Drexler W, Fujimoto JG. Optical coherence tomography technology and applications. Berlin, Heidelberg: Springer; 2008. DOI: 10.1108/02640470510611517 - 3.
Yao G, Wang LV. Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media. Physics in Medicine and Biology. 1999; 44 :2307-2320 - 4.
Sherif SS, Rosa CC, Flueraru C, et al. Statistics of the depth-scan photocurrent in time-domain optical coherence tomography. Journal of the Optical Society of America. A. 2008; 25 :16-20 - 5.
Yadlowsky MJ, Schmitt JM, Bonner RF. Multiple scattering in optical coherence microscopy. Applied Optics. 1995; 34 :5699-5707 - 6.
Jacques SL, Wang L. Monte Carlo modeling of light transport in tissues. Optical-Thermal Response of Laser-Irradiated Tissue. 2013; 2607 :73-100 - 7.
Chen NG, Bai J. Estimation of quasi-straightforward propagating light in tissues. Physics in Medicine and Biology. 1999; 44 :1669-1676 - 8.
Chen N. Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry. Applied Optics. 2007; 46 :1597 - 9.
Wilson BC, Adam G. A Monte Carlo model for the absorption and flux distributions of light in tissue. Medical Physics. 1983; 10 :824-830 - 10.
Biondini G, Kath WL, Menyuk CR. Importance sampling for polarization-mode dispersion. IEEE Photonics Technology Letters; 2002; 14 (3):310-312 - 11.
Lima IT et al. Efficient computation of outage probabilities due to polarization effects in a WDM system using a reduced stokes model and importance sampling. IEEE Photonics Technology Letters. 2003; 15 (1):45-47 - 12.
Schmitt JM, Ben.-Letaief K. Efficient Monte Carlo simulation of confocal microscopy in biological tissue. Journal of the Optical Society of America. A. 1996; 13 :952-961 - 13.
Iwabuchi H. Efficient Monte Carlo method for radiative transfer modeling. Journal of the Atmospheric Sciences. 2006; 63 :2324-2339 - 14.
Lima IT Jr, Kalra A, Sherif SS. Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography. Biomedical Optics Express. 2011; 2 (5):1069-1081. DOI: 10.1364/BOE.2001069 - 15.
Lima IT, Kalra A, Hernández-Figueroa HE, Sherif SS. Fast calculation of multipath diffusive reflectance in optical coherence tomography.Biomedical Optics Express. 2012; 3 :692-700 - 16.
Fang Q. Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates. Biomedical Optics Express. 2010; 1 :165-175 - 17.
Malektaji S, Lima IT Jr, Sherif SS. Monte Carlo simulation of optical coherence tomography for turbid media with arbitrary spatial distributions. Journal of Biomedical Optics. 2 April 2014; 19 (4):046001. DOI: 10.1117/1.JBO.19.4.046001 - 18.
Malektaji S, Lima IT, Escobar I, et al. Massively parallel simulator of optical coherence tomography of inhomogeneous turbid media. Computer Methods and Programs in Biomedicine. 2017; 150 :97-105 - 19.
Lima IT, Kalra A, Hernández-Figueroa HE, et al. Fast calculation of multipath diffusive reflectance in optical coherence tomography. Biomedical Optics Express. 2012; 3 :692 - 20.
Papoulis A. Probability, Random Variables, and Stochastic Processes. New York, NY: McGraw-Hill; 1984