1.

Introduction

Optical coherence tomography (OCT) is a rapidly growing noninvasive imaging technique with an increasing number of biomedical applications.^1–6 The study of physical processes underlying OCT imaging of different objects is mainly due to the experimental studies. Therefore, the development of a general and fast OCT imaging simulator would greatly facilitate these studies. This advanced simulator would also be a powerful tool to design novel OCT systems with improved performance, e.g., increased depth of imaging. Earlier, we developed fast Monte Carlo methods to compute OCT signals from a turbid multilayered object.^7,8 In this article, we generalize our previous work to include objects with arbitrary spatial distributions.

Many available OCT simulators^7–12 are based on the Monte Carlo simulation of light transport in multilayered turbid media (MCML) that was developed by Wang et al.¹³ MCML is still a popular simulator. However, its main drawback is its restriction to multilayered media.

A simulation of OCT imaging of multilayered objects whose layer boundaries are given by different mathematical functions, instead of planes, was implemented by Kirillin et al.¹⁴ This approach to modeling nonplanar multilayered media has been used to simulate OCT imaging of paper,¹⁴ human enamel,¹⁵ and polarization-sensitive OCT images of human skin.^16,17

To develop a simulator of OCT imaging of arbitrary-shaped media, we need to first consider available simulation methods of light transport in such media. Among the first efforts to model arbitrary-shaped media was the cubic voxelization method introduced by Pferer et al.¹⁸ In this method, the medium is modeled as a set of voxels with different optical properties.^19–21 Another method to model the shape of turbid media uses standard geometrical building blocks such as ellipsoids, cylinders, and polyhedrons.²² A more recent approach to modeling turbid media is based on a surface mesh method proposed by Cote and Vitkin²³ and Margallo-Balbas and French.²⁴ In this approach, the surface of every homogeneous region of the turbid media is represented by a surface mesh. But a high-computational cost is involved in locating the intersection of the path of a photon with its enclosing surface mesh. Therefore, Fang proposed a Plücker coordinate system-based, mesh-based Monte Carlo method for fast computation of the location of such intersections.²⁵

To locate the intersection of the path of a photon with the enclosing surface mesh, one needs to find intersections with all planes comprising the surface mesh. If the enclosing surface mesh is convex, then finding such intersections would be considerably simpler. Since a tetrahedron volume is convex with minimum number of planes, its usage as a building block for a surface mesh minimizes the computational cost for simulation of light in arbitrary-shaped media. A Monte Carlo simulator of light transport in arbitrary-shaped media modeled with tetrahedrons is tetrahedron-based inhomogeneous Monte Carlo optical simulator (TIM-OS) that was developed by Shen and Wang.²⁶ Our OCT simulator is based on the TIM-OS code. A comprehensive review of methods to simulate light transport in turbid media can be found in Zhu and Liu.²⁷

2.

Photon Tracing in Turbid Media with Tetrahedron-Based Mesh

In TIM-OS, each arbitrary-shaped region defined by optical parameters, scattering coefficient ${\mu }_s$, absorption coefficient ${\mu }_a$, refractive index $n$, and anisotropy factor $g$, is divided into a number of tetrahedrons. A large number of successive pencil photon packets are launched into the medium from the origin of a rectangular coordinate system ($x,y,z$), where the $z$-axis is normal to the surface of the medium. The initial weights of these photon packets are unity, $W=1$. A photon packet will travel a distance equal to the mean free path after which it undergoes absorption and scattering. The photon packet rates of absorption and scattering depend on the absorption and scattering coefficients of the regions where this packet is traveling, respectively. If a traveling photon packet enters a region with a different refractive index, then it will undergo both specular reflection and refraction at the boundary between the two regions. Therefore, at each step, we have to check if the packet path and the enclosing tetrahedron intersect. This process is the most computationally demanding part of the simulation.

Let the photon packet position be $\widehat{p}$, its direction be $\widehat{u}$, and its mean free-path length be $s$. Therefore, any point ${\widehat{p}}^{\prime }$ on the packet path satisfies ${\widehat{p}}^{\text{'}}=\widehat{p}+\widehat{u}t$, where $t\in [0,s]$. Let $\widehat{n}.\widehat{x}+d=0$ be the equation of any plane of the enclosing tetrahedron, where $\widehat{x}$ is any point on this plane, $\widehat{n}$ is its inward normal unit vector, and $d$ is its minimum distance from the origin. If the photon packet path intersects any of the enclosing tetrahedron facets, then the distance from the current packet position to the plane will be given by

The minimum positive-valued $t$ associated with intersections with each tetrahedron plane indicates the point of intersection of the photon packet path and the enclosing tetrahedron.

3.

Photon Tracing in Turbid Media with Tetrahedron-Based Mesh

Time-domain OCT signals consist of three types of photons collected by the probe: (1) ballistic photons that are single scattered, (2) quasiballistic photons that are multiple scattered within the coherence length of the optical source, and (3) multiple scattered photons beyond the coherence length of the optical source. Ballistic and quasiballistic photons contribute to Class I diffusive reflectance. The multiple scattered photons beyond the coherence length contribute to Class II diffusive reflectance. It has been shown that the Class II diffusive reflectance is a fundamental limit for OCT imaging depth.^28,29

In our simulator, we used a fiber probe with radius $d_{\max }$ and acceptance angle ${\theta }_{\max }$. To calculate the Class I diffuse reflectance, we define a spatial–temporal indicator function as

We calculate Class I diffusive reflectance, $R_1(\mathrm{z})$, and Class II diffusive reflectance, $R_2(\mathrm{z})$, from a particular depth, $z$, as the weighted mean of these indicator functions

An estimate of the variance of these estimations could be calculated as follows:

4.

Importance Sampling for Reducing Computational Time of OCT Signals

Importance sampling is a technique used to reduce the variance of results obtained by Monte Carlo methods using the same number of statistical samples. Therefore, importance sampling could be used to bias a Monte Carlo method to reduce the computational cost of a result with a required accuracy. Light propagating in turbid tissue is usually scattered in the forward direction, i.e., typical turbid tissue has a high-anisotropy factor and the probability of backscattered events is very small. Therefore, a very large number of photon packet simulations are required to estimate OCT signals in a standard, i.e., without importance sampling, Monte Carlo simulation. By properly biasing the scattering direction, we could considerably increase the number of such backscattered events.

A photon packet in turbid tissue will scatter according to the Henyey–Greenstein phase function,

The importance sampling technique used in this article is similar to the one described in Lima et al.⁸ If a photon packet is traveling away from the probe ($u_z>0$), then we bias its scattering direction toward the actual position of the probe,$\widehat{v}$, to increase the probability of its detection. For the first biased scattering event, we use the following probability density function to sample the biased longitudinal scattering angle ${\theta }_B$:

To maintain unbiased estimates of diffusive reflectance, we assign a compensating likelihood value to this biased scattering event,

After the first biased scattering, the photon packet can undergo unbiased scatterings with probability p or other biased scattering with probability1–p. For biased scattering, we sample the longitudinal scattering angle from a Henyey–Greenstein phase function that is oriented toward the actual position of the probe and with anisotropy factor a. For both biased and unbiased scattering events, we use the following likelihood function:

After the tracing of a photon packet ends, due to its collection by the probe, its exit from the medium or being absorbed by it, the total likelihood of the photon packet is calculated. This total likelihood $L$ is the product of all likelihood values assigned to all scattering events undergone by the photon packet. If $L<1$, another packet with initial likelihood value of $L^{\prime }=1-L$ will be launched from the previous first scattering event position with a direction sampled from the unbiased Henyey–Greenstein phase function.⁸

5.

Implementation and Validation of Our OCT Simulator

We implemented our OCT simulator in American National Standards Institute (ANSI)-C and used a random number generator included in the GNU scientific library (GSL).³⁰ We used an optical fiber probe with radius 1 μm, acceptance angle 5 deg, and an optical source with coherence length $l_c=0.015\mathrm{mm}$. We also implemented the importance sampling scheme described above with bias coefficient $a=0.925$ and additional bias probability $p=0.5$, respectively. We ran our simulations on a 2 GHz Intel Core i7 CPU with 4 GB of RAM. Each A-scan simulation used ${10}^7$ photon packets. Our simulation results were validated by comparing them with results obtained from the previously verified OCT simulators for multilayered media by Yao and Wang⁹ and Lima et al.⁸ Since the simulator by Yao and Wang⁹ does not use the advanced importance sampling method in Sec. 4, ${10}^9$ photon packets were used to obtain results with comparable accuracy.

5.1.

Simulation of OCT Signal from a Multilayered Object

As a multilayered object, we used a medium consisting of four layers with optical properties, as shown in Table 1, similar to objects used in Refs. 7–9. In our novel tetrahedron-based simulator, this multilayered object was divided into 9600 tetrahedrons that resulted in 2205 vertices.

Table 1

Optical parameters of the multilayered object used to validate our new tetrahedron-based OCT simulator.

Layers	Height (cm)	Scattering coefficient ${\mu }_s({\mathrm{cm}}^{-1})$	Absorption coefficient ${\mu }_a({\mathrm{cm}}^{-1})$	Anisotropy factor $g$	Refractive index $n$
1	0.0200	60	1.5	0.9	1.0
2	0.0015	120	3	0.9	1.0
3	0.0345	60	1.5	0.9	1.0
4	0.0440	120	3	0.9	1.0

Figures 1 and 2 show the A-scans obtained from our object using two previously validated OCT simulators for multilayered objects and our novel tetrahedron-based OCT simulator for arbitrary-shaped objects. As seen in Figs. 1 and 2, the results from the three simulators are in excellent agreement, which validates the results from our new simulator. One would expect the computational cost to simulate OCT signals from a multilayered object using a tetrahedron-based simulator would be higher than using a layer-based one. The layer-based simulator by Lima et al.⁸ took 24 min to calculate Class I and Class II diffusive reflectances for a single A-scan, whereas our new tetrahedron-based simulator took 43 min.

Fig. 1

A-scan of Class I diffusive reflectance from the multilayered object above. The green and blue lines are results from the previously validated layered-based OCT simulators by Yao and Wang⁹ and Lima et al.,⁸ respectively. The red line represents the Class I signal of our novel tetrahedron-based OCT simulator.

Fig. 2

A-scan of Class II diffusive reflectance from the multilayered object above. The green and blue lines are results from the previously validated layered based OCT simulators by Yao and Wang⁹ and Lima et al.⁸ respectively. The red line represents the Class II signal of our novel tetrahedron-based OCT simulator.

5.2.

Simulation of OCT Signal from a Nonlayered Object

To demonstrate the ability of our new OCT simulator to simulate signals from nonlayered objects, we simulate OCT imaging of a sphere inside a homogeneous slab. As shown in Fig. 3, we placed a sphere of radius 0.1 mm, at depth 0.2 mm, inside a slab with $3\times 3\mathrm{mm}$ lateral dimensions and 1 mm axial dimension. The optical properties of this nonlayered object are shown in Table 2.

Fig. 3

Spatial structure of our nonlayered object.

Table 2

Optical parameters of the nonlayered object used in our new tetrahedron-based OCT simulator.

Medium	Absorption coefficient ${\mu }_a({\mathrm{cm}}^{-1})$	Scattering coefficient ${\mu }_s({\mathrm{cm}}^{-1})$	Anisotropy factor $g$	Refractive index $n$
Slab	1.5	60	0.9	1
Sphere	3	120	0.9	1

We used the mesh generator NETGEN³¹ to generate 4437 tetrahedrons and 800 vertices to represent our nonlayered object. To obtain an OCT B-scan, we simulated 512 equidistant A-scans along the $x$-axis from $x=–0.15$ to $x=0.15\mathrm{mm}$. A depiction of these tetrahedrons and the imaged cross-section of the object are shown in Fig. 4. The simulation of a complete B-scan took approximately 360 h on our computer.

Fig. 4

A depiction of the tetrahedrons representing our nonlayered object. The solid black lines shown in the imaged cross-section (B-scan) depict intersections of these tetrahedrons with the imaging cross-section.

Figures 5(a) and 5(b) show the simulated B-scan OCT images from our nonlayered object, each comprised 512 A-scans along the $x$-axis. Figure 5(a) represents the Class I OCT signal resulting from single-scattered photons, while Fig. 5(b) represents the Class II OCT signal resulting from multiple-scattered photons. We note that the Class I diffusive reflectance-based image in Fig. 5(a) closely resembles the object, a sphere inside a slab. Also, as expected, being an image based on single-scattered light, its intensity is considerably reduced as the imaging depth increases. From Fig. 5(b), we note that, as expected, the intensity of the Class II diffusive reflectance increases with depth, particularly inside the sphere, but as the imaging depth increases, it gets weaker due to absorption. All these results are as one would expect from a physical OCT, which further validates our new simulator.

Fig. 5

Simulated (a) Class I and (b) Class II reflectance-based B-scan OCT image of our nonlayered object.

5.3.

Resolution, Accuracy, and Computation Time of Our OCT Simulator

In the following subsections, we discuss three important aspects of our simulator: (1) resolution, (2) accuracy, and (3) computation time.

5.3.1.

Simulation resolution

The lateral resolution of physical OCT systems depends on the used wavelength and numerical aperture of the imaging optics. Our simulated lateral resolution of a B-scan is, however, equal to the distance between our simulated A-scans.

Figures 6(a) and 6(b) show the simulated Class I and Class II B-scans of our nonlayered object using 75 equidistant A-scans, along the $x$-axis from $x=–0.15$ to $x=0.15\mathrm{mm}$. These B-scans have a lower lateral-simulated resolution than the B-scans, comprised 512 A-scans, as shown in Fig. 5. As shown in Fig. 6, the sphere is still distinguishable despite the lower lateral simulation resolution. As the computation time to simulate a B-scan is linearly related to the number of simulated A-scans, the results in Fig. 6 were obtained in approximately 53 h. Therefore, the lateral-simulated resolution, i.e., the number of A-scans to simulate, should be chosen depending on the object of interest.

Fig. 6

Simulated (a) Class I and (b) Class II reflectance-based B-scan OCT images of our nonlayered object using 75 A-scans.

The axial resolution of physical OCT systems is approximately equal to the coherence length of the optical source. Our simulated axial resolution of an A-scan is, however, determined according to the Eq. (3). If a collected photon has traveled an optical distance $\mathrm{\Delta }{s}_i$ and has reached a maximum depth, $z_{\max }$, satisfies the conditions

Eq. (10)

$$l_{\mathrm{c}}>|\mathrm{\Delta }{s}_{\mathrm{i}}-2z_{\max }|,$$

or

Eq. (11)

$$l_{\mathrm{c}}<|\mathrm{\Delta }{s}_{\mathrm{i}}-2z_{\max }|,$$

we classify it as Class I or Class II photon, respectively. Next, we assign a depth, $z$, to this photon which satisfies the following inequality according to Eq. (3):

Eq. (12)

$$\frac{\mathrm{\Delta }{s}_{\mathrm{i}}-l_{\mathrm{c}}}{2}\le z\le \frac{\mathrm{\Delta }{s}_{\mathrm{i}}+l_{\mathrm{c}}}{2}.$$

Normally, one would assign a single A-scan pixel to represent this depth range. Our simulator, however, allows over sampling of simulation depth by dividing the depth range in Eq. (12) into a number of subresolution depth ranges and then assigning an additional pixel to each of these subresolution depth ranges. Therefore, our simulated axial resolution depends on the coherence length of the source, similar to physical OCT systems, and on the defined number of subresolution depth ranges, i.e., the required oversampling ratio. We used an oversampling ratio of 6 in all our simulations above. As our simulation computation time is dominated by the number of simulated photons, rather than the coherence length or oversampling ratio, any simulated axial resolution value will have a negligible impact on computation time.

5.3.2.

Simulation accuracy

Another very important aspect of our simulator is the accuracy of its results. The accuracy of Monte Carlo simulations is generally quantified by the variance of obtained results. In Fig. 7(a), we show estimates using different number of photons of Class I signals of an A-scan of our nonlayered object, in addition to confidence intervals (CIs) of these estimates. The CI of Class I and Class II signal estimates ${\mathrm{CI}}_{1,2}(z)$, is defined as

Eq. (13)

$${\mathrm{CI}}_{1,2}(z)\equiv [R_{1,2}(z)-{\sigma }_{1,2}(z),R_{1,2}(z)+{\sigma }_{1,2}(z)].$$

Fig. 7

(a) Class I signal estimates and their confidence intervals, (b) signal-to-computational-noise ratios of Class I signal estimates, and (c) signal-to-computational-noise of Class I signal estimates using importance sampling with ${10}^7$ and ${10}^5$ photons and without importance sampling using ${10}^7$ photons.

We note from Fig. 7(a) that as the depth increases, the CIs increase relative to the signal values. Therefore, as reported by Yao and Wang,⁹ the accuracy of our signal estimates relative to signal values decreases with depth.

Another measure of the accuracy of our estimated Class I and Class II signals is the signal-to-computational-noise-ratio, which is given by

Eq. (14)

$${\mathrm{SNR}}_{1,2}(z)=\frac{R_{1,2}(z)}{{\sigma }_{1,2}(\mathrm{z})}.$$

In Fig. 7(b), we show signal-to-computational-noise-ratios, ${\mathrm{SNR}}_1(z)$, of Class I signal estimates, using different number of photons, of an A-scan ($x=0$) of our sphere inside a slab object.

We note from Fig. 7(b) that, as expected in Monte Carlo simulations, the signal-to-computational-noise-ratio is proportional to the square root of the number of simulated photons, $N$, i.e.,

Eq. (15)

$${\mathrm{SNR}}_{1,2}(z)\propto \sqrt{N}.$$

We also note that our simulation computation time is linearly proportional to $N$.

To illustrate the effect of our importance sampling technique, in Fig 7(c), we compare the signal-to-computational-noise-ratio of three A-scans (Class I) at $x=0$ of our sphere inside a slab object. Two of these A-scans were obtained with our importance sampling technique using ${10}^7$ and ${10}^5$ photons. The third A-scan was obtained using ${10}^7$ photons but without using importance sampling. As seen in Fig 7(c), our importance sampling technique has improved the accuracy of the simulation by approximately 100 times.

5.3.3.

Simulation computation time

Different approaches could be used to reduce the computation time of our simulations. First, we could simulate a smaller number of A-scans, which would result in reduction of lateral resolution in the simulated B-scan images. Second, we could use a smaller number of launched photons per A-scan, which would result in simulation results with lower accuracy. Third, because the inherently parallel nature of Monte Carlo simulations, we could reduce our simulation time by implementing our simulator on graphics processor units (GPUs). Simulations of light transportation in turbid media^32–34 and of OCT in multilayered media³⁵ have been implemented on GPUs using the compute unified device architecture (CUDA) programming language. In addition to parallel implementation, a hardware implementation on field-programmable gate array (FPGA) of simulation of light transport in turbid media has been shown to reduce computation time.³⁶ This FPGA-based approach could also be used to reduce computation time of our simulator.

6.

Conclusions

We developed a novel Monte Carlo-based simulator of OCT imaging for turbid media with arbitrary spatial distributions. This simulator allows computation of both Class I and Class II diffusive reflectance-based OCT images. It was implemented using a tetrahedron-based mesh that allows modeling of an arbitrary-shaped medium with a desired accuracy. This mesh also reduces the computation cost required to obtain the intersection of a photon path with its enclosing tetrahedron. We used Monte Carlo importance sampling to significantly increase the probability of a photon reaching the optical detector, thereby reducing simulation time.

Our simulation results were verified by comparing them to results from two previously validated OCT simulators for multilayered media. We also presented simulation results for OCT imaging of a sphere inside a background slab, which would not have been possible with earlier simulators. We also discussed the resolution and accuracy of our simulator and suggested different ways to reduce computation time. Our simulator could be used to study important OCT phenomena and to design novel OCT systems with improved performance. As future work, we plan to implement our simulator on GPUs using the CUDA programming language, as well as develop a similar simulator for swept-source OCT systems.