Fast calculation of multipath diffusive reflectance in optical coherence tomography.

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.

1. Introduction

Optical coherence tomography (OCT) is a sub-surface imaging technique that produces tissue images 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 [1], and can be integrated with optical fiber probes and also allow spectroscopic characterization of tissue [2]. Computer simulations of light transport in multi-layered turbid media are an effective way to theoretically investigate light transport in organic tissues [1], which can also be applied in the development of optical coherence tomography (OCT) systems [2].

Yao and Wang [3] proposed an importance sampling implementation to speed up Monte Carlo simulations in time-domain OCT (TD-OCT) systems. In [4], we reviewed previous methods proposed to speed up Monte Carlo simulations and we showed how one can extend the method in [3] using importance sampling to significantly speed up the calculation of the OCT signal from the entire depth range imaged that results from the ballistic and the quasi-ballistic scattered photons, the Class I reflectance, without significant residual bias in the statistical result. That importance sampling method consists of applying multiple biases towards the source after the first biased scattering is applied to the photon packets, and using a photon splitting procedure. The Class I reflectance is produced by the photons that only undergo ballistic or quasi-ballistic propagation, besides undergoing back-scattering at the layer being imaged. Importance sampling, tailored to each particular application, is a statistical technique that has also been applied in optical communications [5-7], confocal microscopy [8], atmospheric optics [9], and in diffuse optical tomography [10].

In this paper we show that the implementation used in [4] can be enhanced to also efficiently calculate the OCT signal due to the detected photons that undergo multiple diffusive scattering in tissue, the Class II reflectance. The Class II reflectance interferes with the Class I reflectance and, as a consequence, limits the ability of OCT systems to image deeper into the tissue [11]. The importance sampling method in [4] consists of randomly applying multiple biases from the first biased scattering towards the tip of the collecting optical fiber. In the method proposed here, the application of additional biased scatterings towards the collecting optics, after the first biased scattering, are not carried out automatically as in the previous method. Instead, the bias is applied with a probability 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.

In Sec. 2, we described the variance reduction method that we developed to speed up the calculation of the reflectance using Monte Carlo simulations. In Sec. 3, we show numerical results for the standard Monte Carlo method and for our importance sampling implementations with different parameters. In that section, we validated the importance sampling method by comparing it against extensive standard Monte Carlo simulations, and we demonstrate the effectiveness of our importance sampling implementation.

2. Importance sampling applied to OCT

The reflectance in turbid media, such as organic tissue, obtained with OCT systems can be easily calculated by modeling light as electromagnetic fields using the first Born approximation [13], without which computer simulations using the full field model would not be practical. However, this approximation only includes the contribution of the ballistic photons. Therefore, this simplified model cannot account for the photons that undergo multiple diffusive scattering in tissue, which can limit the imaging depth of OCT systems. One alternative to this method is the use of the Monte Carlo method do model light transport in tissue, in which the light is modeled as photon particles that do not interact with one another. This method is applicable when the coherence length of the light is short when compared to both the mean path length in tissue and the widths of the tissue layers, which corresponds to most practical cases in which OCT imaging is exploited.

We implemented our importance sampling method for Monte Carlo modeling of light transport in multi-layered tissues (MCML) by modifying a C-language software package that is available for download from the web site of the Oregon Medical Laser Center [14,15]. MCML can be used to simulate an ensemble of photon packets launched in a steady-state one-dimensional beam, normal to the surface of the topmost tissue layer. The scattering events that arise from the simulated tissue are characterized by two random angles that determine the future direction of the photon packet in three-dimensional space after the scattering. As in the MCML, to account for the photon packet scattering with arbitrary anisotropy factor g of the tissue, in the unbiased scatterings we use the Henyey-Greenstein probability density function that is defined as

where θ is the angle between the photon packet propagation direction prior to the scattering and the new scattering direction . After rotating away from the previous propagation direction by the angle , so that , the scattering direction is rotated around by an angle ϕ that is randomly picked from a uniform probability density function from 0 to 2π.

2.1. First biased scattering

Considering that the probability that a photon packet is scattered towards the apparent position of the collecting optics at any given layer is very low, since the anisotropy factor 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]. In [4] we proposed an importance sampling implementation for the Class I reflectance that was based on biased scattering towards the apparent position of the collecting optics and a photon splitting procedure followed by successive biased scatterings towards the apparent position of the collecting optics, whose direction we defined as the unit vector . Because that method was not as efficient in the calculation of the Class II reflectance, we significantly enhanced that method by implementing two modifications. The first modification consists of using a biased probability density function for the first biased scattering that produces random scattering with an angle no lesser than 90 degrees away from the direction to the apparent position of the collecting optics. This biased probability density function, which was based on the Henyey-Greenstein probability density function in Eq. (1), is given by

where a is a bias coefficient in the range (0,1). After randomly picking a biased angle away from the direction of the apparent position of the collecting optics, the biased direction, so that, the resultant biased scattering direction is rotated around 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 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-7], which ensure converge of the statistical result towards the expected value. The likelihood ratio of the photon packet using the biased probability density function in Eq. (2) is given by

where is a function of , which is statistically drawn from the probability density function in Eq. (2) that is used to define the new propagation direction of the photon packet at the first biased scattering. A schematic representation of the angles and vectors used in the biased and in the unbiased scatterings is shown in Fig. 1 of [4].

Fig. 1

Schematic representation of a simulation setup similar to [5].

Other probability density functions can also effectively speed up the calculation using this method, provided that they significantly increase the probability that the photon packet is scattered towards the apparent position of the collecting optics.

2.2. Additional biased scatterings

The second enhancement that we made to the variance reduction method in [4] consists of applying additional biased scatterings towards the apparent direction of the collecting fiber with probability 0 ≤ 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

If the biased function is selected to draw a random value of , which is an event with probability p, is a function of that is statistically drawn from the probability density function in Eq. (2). Otherwise, which is a complementary event with probability 1 – p, the unbiased function is selected to draw a random value of and is a function of . Because each random scattering is independent from the other scatterings, the likelihood ratio of each photon packet is equal to the product of all the likelihood ratios of all the biased scatterings of that photon packet.

2.3. Calculation of the reflectance

The Class I and the Class II reflectances at depth z are obtained by calculating the mean value of the indicator functions I1 and I2 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 I1 and I2 of the spatial and temporal filter for the ith photon packet is defined as

and

where z is the depth imaged, l is the coherence length of the source, r is the distance of the ith reflected photon packet to the origin along the plane z = 0, where the collecting optical system is located, dmax and θmax are the maximum collecting diameter and angle, respectively, θ is the angle with the z-axis, which is the axis normal to the tissue interface, Δs is the optical path, zmax 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]. The estimated values of the Class I and Class II reflectances and their respective variances are given by the following expressions

and

where W(i) is the weight of the ith 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 ith 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

To generate random angles according to the biased probability density function in Eq. (2) we use the pseudo-random number generator of the Gnu Scientific Library [16]. That random number generator produces pseudo-random numbers uniformly distributed from 0 to 1, which we convert to the probability density function in Eq. (2) according to the following conversion formula

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

3. Numerical results

We validate the importance sampling method for Monte Carlo simulations that we developed by comparing it against extensive standard Monte Carlo simulations. We simulate a turbid media that consists of multiple layers, shown schematically in Fig. 1. The tissue extends from 0 to 1 mm, consisting primarily of a turbid layer with absorption coefficient µ = 1.5 cm−1 and a scattering coefficient µ = 60 cm−1, but also contains five thin layers with absorption coefficient µ = 3 cm−1 and a scattering coefficient µ = 120 cm−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]. In [4] we showed that our method is robust in the presence of refractive index mismatch along the tissue, in which the apparent position of the collecting optics is different in each layer due to refraction at the interfaces. 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. For simplicity, the light source is assumed to be a one-dimensional light beam propagating along the vertical direction as in [3,4,10], since the purpose of this study is to validate and demonstrate the effectiveness of the importance sampling implementation that we developed.

In Figs. 2 and 3 , we show results obtained with 108 Monte Carlo simulations with importance sampling, which has a computational cost of about 9 × 108 standard Monte Carlo simulations in this particular simulation setup because of the computational cost of the photon splitting procedure [4]. The computational cost increase of this importance sampling method depends on the target depth range and on the photon mean free path in the tissue. The target depth range of these simulations was set from 0 mm to 1 mm. Therefore, every single photon scattering that occurs in the depth range from 0 mm to 1 mm is biased. We run the Monte Carlo simulations with importance sampling with the bias coefficient 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.

Fig. 2

The Class I reflectance, shown with thick solid black curve, and the Class II reflectance, shown with thin solid red curve, as a function of the depth for the importance sampling implementation described in Sec. 2 with 108 samples. The pink short dashed curve and the blue long dashed curve are results of 1011 standard Monte Carlo simulations of the Class I reflectance and the Class II reflectance, respectively.

Fig. 3

The reflectance results shown in Fig. 2 for the depth interval from 640 µm to 680 µm. The error bars shown for every other point were estimated in the same ensemble of simulations.

We used a computer with the AMD Opteron 246 processor with a clock of 2 GHz and 4GB of RAM to run all the simulations presented. The simulation using the standard Monte Carlo method, whose results are shown in Figs. 2 and 3, required eight days, 22 hours, and 7 minutes of computer time, while the simulation using the importance sampling method for Monte Carlo required only 1 hour and 53 minutes of computer time. We observed that the confidence intervals of the results obtained using the standard Monte Carlo method are significantly larger than the ones obtained with importance sampling for the results shown in Fig. 3, even though the standard Monte Carlo simulations required 113 times the computational time of the simulations with importance sampling. It would have been necessary to increase the number of samples in the standard Monte Carlo simulations by one order of magnitude (89 days of computer simulation) to obtain confidence intervals of the Class I reflectance comparable to those obtained using importance sampling. In Fig. 3, we also observed that this method reduces the computational cost of calculating the Class II reflectance by more than two orders of magnitude.

In Fig. 4 we show the dependence of the relative error of the calculation of the Class I and the Class II reflectances at two different depths: 400 µm and 670 µm with respect to the bias coefficient 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.

Fig. 4

The relative error in the calculation of the reflectance using importance sampling as a function of the bias coefficient a at 400 µm and at 670 µm of depth for p = 0.5.

In Fig. 5 we show the dependence of the relative error of the calculation of the Class I and the Class II reflectances at two different depths: 400 µm and 670 µm, with respect to the probability of additional bias 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], in which all the additional scatterings (after the first biased scattering) are biased. In that case, the Class II reflectance at 670 µm of depth has a relative error that is nearly twice the minimum value obtained when 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.

Fig. 5

The relative error in the calculation of the reflectance using importance sampling as a function of the probability of additional bias p at 400 µm and at 670 µm of depth for a = 0.9.

4. Conclusion

We developed a new importance sampling method that enables a reduction of the computation time of the Class I and the Class II reflectances in TD-OCT by as much as three orders of magnitude, which we validated by comparing it against a large ensemble of standard Monte Carlo simulations based on MCML. This enables a computation time reduction from several days down to tens of minutes. This fast Monte Carlo calculation of TD-OCT signals can be a valuable tool in the investigation of the physical process governing both the Class I and the Class II reflectances that may enable the development of novel signal processing techniques to extend the imaging depth of these systems. Since the photons that undergo multiple scatterings in tissue are also responsible for the Class II reflectance in Frequency Domain OCT systems, we believe that this importance sampling algorithm can be easily extended to those systems.