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.
- optical coherence tomography
- Monte Carlo simulation
- light transport in turbid media
- importance sampling
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 . However, multiply scattered photons do not carry practically useful information about the imaged tissue; therefore, they result in a degradation of the OCT signal . In addition, it has been shown that Class II diffuse reflectance represents a fundamental limit related to the imaging depth of OCT in tissue . 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 , atmospheric optics , and diffuse optical tomography .
To improve the computational efficiency of the MC-based simulation of OCT systems , Yao and Wang proposed the first importance sampling technique to simulate OCT signals from a multilayered turbid medium . 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 , 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 . 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 . 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 . 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) . 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 μa and the scattering μs 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, g, of the tissue, MCML uses the Henyey-Greenstein probability density function that is given by
where is the angle between the propagation direction of the photon packet before the current scattering and is the direction of the photon packet after the current scattering. The angle between the previous propagation direction and the new propagation direction is . Therefore, cos() = ⋅. To ensure that the new propagation direction is statistically correct, the provisional scattering direction is rotated around by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction after the current scattering. At the end of each scattering event, the photon packet weight W is reduced according to the step size and the local absorption coefficient . The weight 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/m or is allowed to continue propagating with probability 1−1/m and a new weight equal to m⋅W once the weight reaches . We use the value m = 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.
The Class I diffuse reflectance at depth equal to z is obtained by calculating the mean value of an indicator function that represents a spatial and temporal filter of Class I diffuse reflectance for all simulated photon packets. The indicator function I1 of such spatial and temporal filter for the ith photon packet is defined as
where lc is the optical source’s coherence length; ri is the distance to the origin of the ith reflected photon packet; dmax and θmax are the maximum photon packet collecting diameter and angle, respectively; is the angle with the z-axis (normal to the tissue interface); Δsi is the optical path; and z is the photon packet’s maximum depth.
At any depth, the diffuse reflectance R1 is the expected value of I1 at this depth, and its standard deviation σR,1 could be estimated by
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 .
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 of a scattered photon toward the tip of the light-collecting optical fiber, , as the photon packet reaches the depth range of interest. By defining the origin of a Cartesian coordinate system at the center of the tip of the light-collecting optical fiber, the bias direction in which this fiber is located is defined as , where is the position vector of the scattering location in the tissue.
All photon packets propagating in a direction close to will contribute to the simulated Class I diffuse reflectance with a higher probability. Therefore, this bias direction is more efficient than biasing only in the backward direction, which may not be consistent with the direction of the light-collecting optical fiber. This choice of the bias direction is particularly effective for photon packets propagating deep in the tissue, where such photon packets experience one or more scattering events before they are diffusively reflected.
3.1 Scattering angle due to first event of backscattering
As the photon packet reaches the depth range targeted, the propagation direction of the scattered photon packet is biased toward the bias direction , as opposed to being most likely scattered close to the previous propagation direction as in the practical case with anisotropy g close to 1 and different from the bias toward , the opposite of the direction of propagation, as it is done in . To randomly select the biased angle between the new scattering direction and the biased direction , we use the same probability density function in Eq. (1). However, the bias coefficient does not have to correspond to the anisotropy factor g. Therefore, the probability density function of the biased angle is given by
where a is a bias coefficient. After randomly sampling a biased angle away from the biased direction , so that cos() = ⋅, the provisional scattering direction is rotated around by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction after the current scattering. This procedure ensures a more accurate model of the light scattering in tissue. The difference in the rotation by ϕ 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 , while the rotation in the standard model in is done around the direction prior to the current scattering. After the first biased scattering, this procedure produces the new propagation direction of the photon packet. Afterward, the scattered photon packet is associated with a likelihood ratio as discussed in other applications of this method [10, 11]. Using our biased angle’s probability density function, the likelihood ratio of the photon packet, Eq. (5), is given by
where cos() = ⋅is determined, after the biased scattering, from the randomly sampled values of and ϕ. 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 , the likelihood ratio also depends also on . Figure 1 shows a schematic drawing of these vectors and the angles used in this direction biasing procedure. Note that the choice of bias distribution only affects the speed of convergence of the simulation. Therefore, other biased probability function could also be used to randomly generate the biased scattering toward the bias direction .
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 toward the direction , pointing to the apparent position of the optical collection system, at every scattering point until the photon packet is removed. These additional biases still use both Eqs. (7) and (8). Since the random values drawn for the angle between the scattering direction and the biased direction are independent of each other and are also independent of the previous scattering events, the overall likelihood ratio of a collected photon packet results from the multiplication of all the likelihood ratios of all the biased scattering in that particular simulation.
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 its standard deviation could be calculated with
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  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 . 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 as shown in Section 3.1. To ensure that there is no statistical bias associated with the forward-propagating photon packet that was split, it will be assigned a likelihood ratio , which is a complement to the likelihood of the biased backward scattered photon packet such that to this second photon packet. This second photon packet, only generated if , also undergoes biased backscattering in the tissue at the end of the next step, which could result in another photon packet split, and successive additional biased scatterings toward the tip of the collecting optical fiber until the photon packet propagates beyond the simulation domain. In cases that we investigated, this procedure increased the computational time of each photon packet by five times when compared with a photon packet computed using the standard Monte Carlo method with the same number of launched photon packets 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 μa = 1.5 cm−1, scattering coefficient μs = 60 cm−1, and refractive index n = 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: μa = 1.5 cm−1, μs = 30 cm−1, and n = 1. After the third layer, the medium was assumed to be air: μa = 0 cm−1, μs = 0 cm−1, and n = 1. The anisotropy factor was assumed g = 0.9 for the three diffusive layers.
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 .
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 is randomly selected, away from the direction of the center of the OCT probe , where , the provisional biased scattering direction is rotated around by an 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 is restricted to a maximum deviation from the biased angle by 90°. This ensures that there would not be packets with very large likelihood ratio that could reduce the efficiency of our importance sampling. The likelihood ratio of the photon packet that uses the biased probability density function in Eq. (9) is given by
where . We note that is obtained using the probability density function in Eq. (9), where it is used to obtain the new propagation direction .
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 is a random number distributed uniformly between 0 and 1, a random value for that satisfies Eq. (9) with bias coefficient a could be generated with the following conversion formula
This conversion formula was derived using probability theory .
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 with probability 0 ≤ p ≤ 1. That contrasts with the technique presented in Section 3, in which p = 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
If the biased function is selected to sample a random value of , which is an event with probability p, is a function of that is statistically sampled from the probability density function in Eq. (9). Otherwise, in the case of the complementary event with probability 1 − p, the unbiased function is used to sample a random value of and depends on the value of . Since the two random angles associated to each scattering do not depend on the random angles selected in the previous scatterings, the likelihood ratio of each collected photon packet results from the multiplication of all the likelihood ratios of all the biased scatterings in that simulated photon packet.
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 μa = 1.5 cm−1 and a scattering coefficient μs = 60 cm−1, and also contains five thin layers with absorption coefficient μa = 3 cm−1 and a scattering coefficient μs = 120 cm−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 n = 1 and an anisotropy factor g = 0.9. We note that our method is robust in the presence of refractive index mismatch along layer boundaries of the tissue . 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.
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 a = 0.925, and an additional bias probability p = 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.
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 p = 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).
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 a = 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.
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.
Conflict of interest
The authors declare no conflict of interest.