## 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 * interaction coefficient* of this tissue. This interaction coefficient is equal to the sum of the absorption

*and the scattering*μ

_{a}

*coefficients of this tissue. The scattering events that take place at the end of the random steps are characterized by two random angles that determine the next direction of the photon packet. To account for the photon packet scattering, given an anisotropy factor,*μ

_{s}

*, of the tissue, MCML uses the Henyey-Greenstein probability density function that is given by*g

where * ϕ*, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction

*is reduced according to the step size and the local absorption coefficient*W

*, which is initially set at 1, is proportional to the number of photons in the photon packet. The photon packet is either removed with probability 1/*W

*or is allowed to continue propagating with probability 1−1/*m

*and a new weight equal to*m

*once the weight reaches*m⋅W

*= 10 in this work. This procedure, denoted Russian roulette, is an unbiased technique to end simulation of photon packets that have a negligibly low contribution to the Monte Carlo simulation, so that a new photon packet can be initiated and simulated.*m

The Class I diffuse reflectance at depth equal to * z* is obtained by calculating the mean value of an indicator function

I

_{1}of such spatial and temporal filter for the

*th photon packet is defined as*i

where * l* is the optical source’s coherence length;

_{c}

*is the distance to the origin of the*r

_{i}

*th reflected photon packet;*i

d

_{max}and

θ

_{max}are the maximum photon packet collecting diameter and angle, respectively;

*-axis (normal to the tissue interface); Δ*z

*is the optical path; and*s

_{i}

*is the photon packet’s maximum depth.*z

At any depth, the diffuse reflectance _{1} is the expected value of _{1} at this depth, and its standard deviation _{R,1} could be estimated by

and

where * N* is the number of photons packets used in MC-based simulations.

## 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 * g* close to 1 and different from the bias toward

*. Therefore, the probability density function of the biased angle is given by*g

where * a* is a bias coefficient. After randomly sampling a biased angle

⋅

*, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction*ϕ

*between the model with importance sampling and the standard model is that the rotation in the model with importance sampling is done around the biased direction*ϕ

where cos(* ϕ*. The ratio of the probability of the scattering angle appearing in the biased case with the standard case is the likelihood ratio that is shown in Eq. (6). In addition to depending on

### 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 * N* launched photon packets using importance sampling, the diffuse reflectance

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 * N*. The increase in the computational time of our importance sampling-based implementation, compared to the standard method with the same number of launched photon packets, depends on the average value of the mean free path, and on the width of the target depth range. We note that in our importance sampling implementation, we do not count split photon packets as additional photon packets when determining the value of N in Eqs. (7) and (8), as the use of their corresponding likelihood ratios will generate the correct result. As a photon packet propagates beyond the target region, the packet will propagate using the standard scattering procedure until it is terminated. Once this photon packet is terminated, a new photon packet will be simulated from the OCT probe, as it is the case in the standard MCML. Even though the splitting procedure implies that the cost of simulating a launched photon with this importance sampling method is higher than the computational cost of simulating launched photos using the standard MCML, the computation cost of the Class I diffuse reflectance in our Monte Carlo simulations with importance sampling required as little as one-thousandth of the computational cost required by the standard Monte Carlo method to achieve the same accuracy in the calculated diffuse reflectance.

### 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 * μ* = 1.5 cm

_{a}

^{−1}, scattering coefficient

*= 60 cm*μ

_{s}

^{−1}, and refractive index

*= 1. The second layer, extending from 2.32 to 2.42 mm from the tip of the fiber, has the same absorption and scattering coefficients as the first layer, but its refractive index is*n

*= 1.33. The third layer, extending from 2.42 to 2.62 mm from the tip of the fiber, has the following parameters:*n

*= 1.5 cm*μ

_{a}

^{−1},

*= 30 cm*μ

_{s}

^{−1}, and

*= 1. After the third layer, the medium was assumed to be air:*n

*= 0 cm*μ

_{a}

^{−1},

*= 0 cm*μ

_{s}

^{−1}, and

*= 1. The anisotropy factor was assumed*n

*= 0.9 for the three diffusive layers.*g

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 * a* is the bias coefficient that can be selected between 0 and 1. Once a biased angle

*randomly sampled from a uniform distribution in the range from 0 to 2π. These parameters are defined in the same manner as those used in the biased distribution presented in Section 3. The only difference is that the domain of*ϕ

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 * a* could be generated with the following conversion formula

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 * p* ≤ 1. That contrasts with the technique presented in Section 3, in which

*= 1 (all the additional scatterings were biased). An unbiased scattering is performed in case a bias scattering is not applied in a given point where scattering takes place. The likelihood ratio associated with this scattering is calculated according to the formula*p

If the biased function * p*,

*, the unbiased function*p

### 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 * μ* = 1.5 cm

_{a}

^{−1}and a scattering coefficient

*= 60 cm*μ

_{s}

^{−1}, and 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 to 215 μm, from 365 to 395 μm, from 645 to 660 μm, from 760 to 775 μm, and from 900 to 915 μm. We assume that this tissue has the same refractive index

*= 1 and an anisotropy factor*n

*= 0.9. We note that our method is robust in the presence of refractive index mismatch along layer boundaries of the tissue [19]. We simulate an OCT system where the light is delivered/collected by the tip of an optical fiber having a radius of 10 μm and an acceptance angle of 5°. For simplicity, the light source is assumed to be a one-dimensional light beam propagating along the vertical direction as in [3, 8], since the purpose of this example is to validate and demonstrate the effectiveness of our second importance sampling technique when it is applied to the standard MC simulation.*g

In Figures 5 and 6, we show results obtained with 10^{8} Monte Carlo photon packets with importance sampling, which has a computational cost of simulating about 9 × 10^{8} 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 * a* = 0.925, and an additional bias probability

*= 0.5, to run the Monte Carlo simulations with importance sampling. The results shown in Figures 5 and 6 show that our new importance sampling technique reduces the computational cost for obtaining the Class I diffuse reflectance by approximately three orders of magnitude when compared to the standard Monte Carlo method. This algorithm is optimum when the additional bias probability is equal to*p

*= 0.5. Since only half of the back-scatterings are biased, this choice contributes toward enabling an optimum number of Class II photons to be collected by the tip of the optical fiber.*p

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

*= 0.5. The depths at 400 and 670 μm correspond to tissue regions near the second and third regions with high local reflectance due to the higher local scattering coefficient. The relative variation in the results is the ratio between the square root of the variance, shown in Eq. (4), and the reflectance in Eq. (3).*p

We note that Class I reflectance has a minimum relative error at 400 μm with * a* = 0.925, but the minimum error at 670 μm occurs at

*= 0.95 μm at 670 μm. The deeper the tissue region, the stronger the required bias because of the increase in the number of scatterings with the depth. However, as the bias coefficient is increased toward 1, larger variations in the likelihood ratio lead to an increase in the relative error. We also note that Class II reflectance has its minimum relative error at 400 μm with*a

*= 0.91, while its minimum relative error at 670 μm increased to only about*a

*= 0.925 μm. The optimum amount of bias required by the Classs II diffusive reflectance in both wavelengths is lower than the optimum bias coefficient observed in the Class I reflectance because the number of ballistic and quasi-ballistic photons increases with the bias, which leads to a decrease in the number of collected photon packets that undergo multiple scatterings. Figure 7 also shows that there is a range for the bias parameter*a

*between 0.9 and 0.95 that enables fast calculation of both Class I and Class II reflectances using our importance sampling-based implementation.*a

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