Monte Carlo simulation of optical coherence tomography for turbid media with arbitrary spatial distributions

Abstract. We developed a Monte Carlo-based simulator of optical coherence tomography (OCT) imaging for turbid media with arbitrary spatial distributions. This simulator allows computation of both Class I diffusive reflectance due to ballistic and quasiballistic scattered photons and Class II diffusive reflectance due to multiple scattered photons. It was implemented using a tetrahedron-based mesh and importance sampling to significantly reduce computational time. Our simulation results were verified by comparing them with results from two previously validated OCT simulators for multilayered media. We present simulation results for OCT imaging of a sphere inside a background slab, which would not have been possible with earlier simulators. We also discuss three important aspects of our simulator: (1) resolution, (2) accuracy, and (3) computation time. Our simulator could be used to study important OCT phenomena and to design OCT systems with improved performance.


Introduction
Optical coherence tomography (OCT) is a rapidly growing noninvasive imaging technique with an increasing number of biomedical applications. [1][2][3][4][5][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][8][9][10][11][12] are based on the Monte Carlo simulation of light transport in multilayered turbid media (MCML) that was developed by Wang et al. 13 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. 14 This approach to modeling nonplanar multilayered media has been used to simulate OCT imaging of paper, 14 human enamel, 15 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. 18 In this method, the medium is modeled as a set of voxels with different optical properties. [19][20][21] Another method to model the shape of turbid media uses standard geometrical building blocks such as ellipsoids, cylinders, and polyhedrons. 22 A more recent approach to modeling turbid media is based on a surface mesh method proposed by Cote and Vitkin 23 and Margallo-Balbas and French. 24 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. 25 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. 26 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. 27

Photon Tracing in Turbid Media with Tetrahedron-Based Mesh
In TIM-OS, each arbitrary-shaped region defined by optical parameters, scattering coefficient μ s , absorption coefficient μ 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 bep, its direction beû, and its mean free-path length be s. Therefore, any pointp 0 on the packet path satisfiesp 0 ¼p þût, where t ∈ ½0; s. Let n:x þ d ¼ 0 be the equation of any plane of the enclosing tetrahedron, wherex is any point on this plane,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.

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 θ max . To calculate the Class I diffuse reflectance, we define a spatial-temporal indicator function as where l c is the coherence length of the source, r i is the distance of the i'th-reflected photon packet from the probe, Δs i is the optical path of i'th photon packet, θ z;i is the angle of the photon packet direction with the z-axis, and z max is the maximum depth reached by the photon packet. Similarly, we can define I 2 for Class II diffuse reflectance as We calculate Class I diffusive reflectance, R 1 ðzÞ, and Class II diffusive reflectance, R 2 ð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: where N is the number of simulated photon packets, and LðiÞ is a likelihood ratio due to an importance sampling biasing scheme defined in the next section.

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, where g is the anisotropy factor, and θ s is the longitudinal angle between the photon packet propagation direction u ¼ ðu x ; u y ; u z Þ prior to the scattering and its new scattered directionû 0 . After sampling cos θ s from the above distribution, we rotate the scattering directionû 0 by an azimuthal angle φ sampled from uniform distribution from 0 to 2π.
The importance sampling technique used in this article is similar to the one described in Lima et al. 8 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,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 θ B : where a is the given bias coefficient in the range [0,1].
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 0 ¼ 1 − L will be launched from the previous first scattering event position with a direction sampled from the unbiased Henyey-Greenstein phase function. 8

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). 30 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 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 9 and Lima et al. 8 Since the simulator by Yao and Wang 9 does not use the advanced importance sampling method in Sec. 4, 10 9 photon packets were used to obtain results with comparable accuracy.

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. 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. 8 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.

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 × 3 mm lateral dimensions and 1 mm axial dimension. The optical properties of this nonlayered object are shown in Table 2. We used the mesh generator NETGEN 31 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 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. Figures 5(a) and 5(b) show the simulated B-scan OCT images from our nonlayered object, each comprised 512 Ascans 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.

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.

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 mm. These B-scans have a lower lateral-simulated resolution than the Bscans, 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.
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 Δs i and has reached a maximum depth, z max , satisfies the conditions l c > jΔs i − 2z max j; (10) or l c < jΔs i − 2z max j; (11) Fig. 4 A depiction of the tetrahedrons representing our nonlayered object. The solid black lines shown in the imaged cross-section (Bscan) depict intersections of these tetrahedrons with the imaging cross-section.  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): 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 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. simulated photons, rather than the coherence length or oversampling ratio, any simulated axial resolution value will have a negligible impact on computation time.

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 CI 1;2 ðzÞ, is defined as CI 1;2 ðzÞ ≡ ½R 1;2 ðzÞ − σ 1;2 ðzÞ; R 1;2 ðzÞ þ σ 1;2 ðzÞ: 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, 9 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 In Fig. 7(b), we show signal-to-computational-noise-ratios, 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., 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 Ascan 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.

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][33][34] and of OCT in multilayered media 35 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. 36 This FPGA-based approach could also be used to reduce computation time of our simulator.

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 arbitraryshaped 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.