Research Papers: Imaging

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

[+] Author Affiliations
Siavash Malektaji

University of Manitoba, Department of Electrical and Computer Engineering, 75A Chancellor’s Circle, Winnipeg, Manitoba R3T 5V6, Canada

Ivan T. Lima, Jr.

North Dakota State University, Department of Electrical and Computer Engineering, 1411 Centennial Boulevard, Fargo, North Dakota 58108-6050

Sherif S. Sherif

University of Manitoba, Department of Electrical and Computer Engineering, 75A Chancellor’s Circle, Winnipeg, Manitoba R3T 5V6, Canada

J. Biomed. Opt. 19(4), 046001 (Apr 02, 2014). doi:10.1117/1.JBO.19.4.046001
History: Received November 28, 2013; Revised February 25, 2014; Accepted March 3, 2014
Text Size: A A A

Open Access Open Access

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.

Figures in this Article

Optical coherence tomography (OCT) is a rapidly growing noninvasive imaging technique with an increasing number of biomedical applications.16 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 simulators712 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.1921 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 Vitkin23 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

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 be p^, its direction be u^, and its mean free-path length be s. Therefore, any point p^ on the packet path satisfies p^'=p^+u^t, where t[0,s]. Let n^.x^+d=0 be the equation of any plane of the enclosing tetrahedron, where x^ 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 Display Formula

t=n^.p^+dn^.u^.(1)

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.

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 dmax and acceptance angle θmax. To calculate the Class I diffuse reflectance, we define a spatial–temporal indicator function as Display Formula

I1(z,i)={1,lc>|Δsi2zmax|,ri<dmax,θz,i<θmax,|Δsi2z|<lc0,otherwise,(2)
where lc is the coherence length of the source, ri is the distance of the i’th-reflected photon packet from the probe, Δsi is the optical path of i’th photon packet, θz,i is the angle of the photon packet direction with the z-axis, and zmax is the maximum depth reached by the photon packet. Similarly, we can define I2 for Class II diffuse reflectance as Display Formula
I2(z,i)={l,lc<|Δsi2zmax|,ri<dmax,θz,i<θmax,|Δsi2z|<lc0,otherwise.(3)

We calculate Class I diffusive reflectance, R1(z), and Class II diffusive reflectance, R2(z), from a particular depth, z, as the weighted mean of these indicator functions Display Formula

R1,2(z)=1Ni=1NI1,2(z,i)L(i)W(i).(4)

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

σ1,22(z)=1N(N1)i=1N[I1,2(z,i)L(i)W(i)R1,2(z)]2,(5)
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 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, Display Formula

fHG(cosθs)=1g22(1+g22gcosθs)3/2,(6)
where g is the anisotropy factor, and θs is the longitudinal angle between the photon packet propagation direction u^=(ux,uy,uz) prior to the scattering and its new scattered direction u^. After sampling cosθs from the above distribution, we rotate the scattering direction u'^ 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 (uz>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: Display Formula

fB(cosθs)={(11aa2+1)1a(1a)(1+a22acosθB)3/2,cosθB<[0,1]0,otherwise,(7)
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, Display Formula

L(cosθB)=fHG(cosθs)fB(cosθB)=1g22a(1a)(11aa2+1)(1+a22acosθB1+g22gcosθs)3/2.(8)

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: Display Formula

L(cosθB)=fHG(cosθs)pfB(cosθB)+(1p)fHG(cosθs).(9)

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=1L will be launched from the previous first scattering event position with a direction sampled from the unbiased Henyey–Greenstein phase function.8

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 lc=0.015mm. 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 107 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 Wang9 and Lima et al.8 Since the simulator by Yao and Wang9 does not use the advanced importance sampling method in Sec. 4, 109 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. 79. In our novel tetrahedron-based simulator, this multilayered object was divided into 9600 tetrahedrons that resulted in 2205 vertices.

Table Grahic Jump Location
Table 1Optical parameters of the multilayered object used to validate our new tetrahedron-based OCT simulator.

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.

Graphic Jump LocationF1 :

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 Wang9 and Lima et al.,8 respectively. The red line represents the Class I signal of our novel tetrahedron-based OCT simulator.

Graphic Jump LocationF2 :

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 Wang9 and Lima et al.8 respectively. The red line represents the Class II signal of our novel tetrahedron-based OCT simulator.

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×3mm lateral dimensions and 1 mm axial dimension. The optical properties of this nonlayered object are shown in Table 2.

Graphic Jump LocationF3 :

Spatial structure of our nonlayered object.

Table Grahic Jump Location
Table 2Optical parameters of the nonlayered object used in our new tetrahedron-based OCT simulator.

We used the mesh generator NETGEN31 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.15mm. 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.

Graphic Jump LocationF4 :

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.

Graphic Jump LocationF5 :

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

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.15mm. 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.

Graphic Jump LocationF6 :

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 Δsi and has reached a maximum depth, zmax, satisfies the conditions Display Formula

lc>|Δsi2zmax|,(10)
or Display Formula
lc<|Δsi2zmax|,(11)
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): Display Formula
Δsilc2zΔsi+lc2.(12)

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.

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 CI1,2(z), is defined as Display Formula

CI1,2(z)[R1,2(z)σ1,2(z),R1,2(z)+σ1,2(z)].(13)

Graphic Jump LocationF7 :

(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 107 and 105 photons and without importance sampling using 107 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,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 Display Formula

SNR1,2(z)=R1,2(z)σ1,2(z).(14)

In Fig. 7(b), we show signal-to-computational-noise-ratios, SNR1(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., Display Formula

SNR1,2(z)N.(15)

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 107 and 105 photons. The third A-scan was obtained using 107 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 media3234 and of OCT in multilayered media35 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.

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.

Drexler  W., Fujimoto  J. G., Optical Coherence Tomography. ,  Springer ,  Berlin, New York  (2008).
Popescu  D. P. et al., “Optical coherence tomography: fundamental principles, instrumental designs and biomedical applications,” Biophys. Rev.. 3, (3 ), 155 –169 (2011). 1793-0480 CrossRef
Welzel  J., “Optical coherence tomography in dermatology: a review,” Skin Res. Technol.. 7, (1 ), 1 –9 (2001). 0909-752X CrossRef
Pierce  M. C. et al., “Advances in optical coherence tomography imaging for dermatology,” J. Invest. Dermatol.. 123, (3 ), 458 –463 (2004). 0022-202X CrossRef
Drexler  W. et al., “Ultrahighresolution ophthalmic optical coherence tomography,” Nat Med.. 7, (4 ), 502 –507 (2001). 1078-8956 CrossRef
Fercher  A. F. et al., “In vivo optical coherence tomography,” Am. J. Ophthalmol.. 116, (1 ), 113 –114 (1993). 0002-9394 
Lima  I., Kalra  A., Sherif  S., “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express. 2, (5 ), 1069 –1081 (2011). 2156-7085 CrossRef
Lima  I. et al., “Fast calculation of multipath diffusive reflectance in optical coherence tomography,” Biomed. Opt. Express. 3, (4 ), 692 –700 (2012). 2156-7085 CrossRef
Yao  G., Wang  L., “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol.. 44, (9 ), 2307 –2320 (1999). 0031-9155 CrossRef
Tycho  A. et al., “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt.. 41, (31 ), 6676 –6691 (2002). 0003-6935 CrossRef
Wang  R. K., “Signal degradation by multiple scattering in optical coherence tomography of dense tissue: a Monte Carlo study towards optical clearing of biotissues,” Phys. Med. Biol.. 47, (13 ), 2281 –2299 (2002). 0031-9155 CrossRef
Kirillin  M. Yu. et al., “Monte Carlo simulation of optical clearing of paper in optical coherence tomography,” Quantum Electron.. 36, (2 ), 174 –180 (2006). 1063-7818 CrossRef
Wang  L., Jacques  S. L., Zheng  L., “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed.. 47, (2 ), 131 –146 (1995). 0169-2607 CrossRef
Kirillin  M. Yu. et al., “Visualization of paper structure by optical coherence tomography: Monte Carlo simulations and experimental study,” J. Europ. Opt. Soc. Rap. Public.. 2, , 07031  (2007). 1990-2573 CrossRef
Shi  B. et al., “Monte Carlo modeling of human tooth optical coherence tomography imaging,” J. Opt.. 15, (7 ), 075304  (2013). 0150-536X CrossRef
Kirillin  M. et al., “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express. 18, (21 ), 21714 –21724 (2010). 1094-4087 CrossRef
Meglinski  I. et al., “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett.. 33, (14 ), 1581 –1583 (2008). 0146-9592 CrossRef
Pfefer  T. et al., “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron.. 2, (4 ), 934 –942 (1996). 1077-260X CrossRef
Binzoni  T. et al., “Light transport in tissue by 3D Monte Carlo: influence of boundary voxelization,” Comput. Methods Programs Biomed.. 89, (1 ), 14 –23 (2008). 0169-2607 CrossRef
Boas  D. et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express. 10, (3 ), 159 –170 (2002). 1094-4087 CrossRef
Li  T., “MCVM: Monte Carlo modeling of photon migration in voxelized media,” J. Innov. Opt. Health Sci.. 3, (2 ), 91 –102 (2010). 1793-7205 CrossRef
Li  H. et al., “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method,” Acad. Radiol.. 11, (9 ), 1029 –1038 (2004). 1076-6332 CrossRef
Cote  D., Vitkin  I., “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express. 13, (1 ), 148 –163 (2005). 1094-4087 CrossRef
Margallo-Balbs  E., French  P. J., “Shape based Monte Carlo code for light transport in complex heterogeneous tissues,” Opt. Express. 15, (21 ), 14086 –14098 (2007). 1094-4087 CrossRef
Fang  Q., “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express. 1, (1 ), 165 –175 (2010). 2156-7085 CrossRef
Shen  H., Wang  G., “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol.. 55, (4 ), 947 –962 (2010). 0031-9155 CrossRef
Zhu  C., Liu  Q., “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt.. 18, (5 ), 050902  (2013). 1083-3668 CrossRef
Yadlowsky  M. J., Schmitt  J. M., Bonner  R. F., “Multiple scattering in optical coherence microscopy,” Appl. Opt.. 34, (25 ), 5699 –5707 (1995). 0003-6935 CrossRef
Kirillin  M. Y., Priezzhev  A. V., Myllylä  R. A., “Role of multiple scattering in formation of OCT skin images,” Quantum Electron.. 38, (6 ), 570 –575 (2008). 1063-7818 CrossRef
The Gnu Project, “Gnu scientific library,” http://www.gnu.org/s/gsl/ (15  November 2013).
Schoberl  J., “NETGEN an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Visual. Sci.. 1, (1 ), 41 –52 (1997). 1432-9360 CrossRef
Alerstam  E., Svensson  T., Andersson-Engels  S., “Parallel computing with graphics processing units for high-speed Monte-Carlo simulation of photon migration,” J. Biomed. Opt.. 13, (6 ), 060504  (2008). 1083-3668 CrossRef
Fang  Q., Boas  D. A., “Monte-Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express. 17, (22 ), 20178 –20190 (2009). 1094-4087 CrossRef
Leung  T. S., Powell  S., “Fast Monte Carlo simulations of ultrasound-modulated light using a graphics processing unit,” J. Biomed. Opt.. 15, (5 ), 055007  (2010). 1083-3668 CrossRef
Kalra  A., Lima  I.  Jr., Sherif  S., “Almost instantaneous Monte Carlo calculation of optical coherence tomography signal using graphic processing unit,” in  Proceedings of IEEE Photonics Conference (IPC) 2013 ,  Bellevue, Washington , paper MA2.3 (2013).
Lo  W. C. Y. et al., “Hardware acceleration of a Monte Carlo simulation for photodynamic therapy treatment planning,” J. Biomed. Opt.. 14, (1 ), 014019  (2009). 1083-3668 CrossRef

Siavash Malektaji received his BSc degree in computer engineering from K. N. Toosi University of Technology, Tehran, Iran, in 2011. He is currently pursuing his MSc degree in electrical and computer engineering at the University of Manitoba, Winnipeg, Canada. His research interests include optical coherence tomography (OCT) and Monte Carlo methods.

Ivan T. Lima Jr. received his PhD degree in electrical engineering from the University of Maryland, Baltimore County, USA, in 2003. He is an associate professor with tenure in the Department of Electrical and Computer Engineering at North Dakota State University. His research interests include optical communications and biomedical engineering. He is a member of SPIE, and he is also a senior member of the IEEE Photonics Society and the IEEE EMBS.

Sherif S. Sherif received his PhD degree in optics from the University of Colorado-Boulder, Boulder, Colorado, and an MS degree in digital image processing from the University of Wisconsin-Madison, Madison, Wisconsin. After postdoctoral training at University of Oxford and Imperial College London, he was a lecturer at the University of Kent, UK. He is currently an associate professor of electrical engineering at the University of Manitoba, Winnipeg, Canada. He has over 75 scientific publications, 3 patents, and he was the recipient of 3 teaching awards.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation

Siavash Malektaji ; Ivan T. Lima, Jr. and Sherif S. Sherif
"Monte Carlo simulation of optical coherence tomography for turbid media with arbitrary spatial distributions", J. Biomed. Opt. 19(4), 046001 (Apr 02, 2014). ; http://dx.doi.org/10.1117/1.JBO.19.4.046001


Figures

Graphic Jump LocationF1 :

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 Wang9 and Lima et al.,8 respectively. The red line represents the Class I signal of our novel tetrahedron-based OCT simulator.

Graphic Jump LocationF2 :

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 Wang9 and Lima et al.8 respectively. The red line represents the Class II signal of our novel tetrahedron-based OCT simulator.

Graphic Jump LocationF3 :

Spatial structure of our nonlayered object.

Graphic Jump LocationF4 :

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.

Graphic Jump LocationF5 :

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

Graphic Jump LocationF6 :

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

Graphic Jump LocationF7 :

(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 107 and 105 photons and without importance sampling using 107 photons.

Tables

Table Grahic Jump Location
Table 1Optical parameters of the multilayered object used to validate our new tetrahedron-based OCT simulator.
Table Grahic Jump Location
Table 2Optical parameters of the nonlayered object used in our new tetrahedron-based OCT simulator.

References

Drexler  W., Fujimoto  J. G., Optical Coherence Tomography. ,  Springer ,  Berlin, New York  (2008).
Popescu  D. P. et al., “Optical coherence tomography: fundamental principles, instrumental designs and biomedical applications,” Biophys. Rev.. 3, (3 ), 155 –169 (2011). 1793-0480 CrossRef
Welzel  J., “Optical coherence tomography in dermatology: a review,” Skin Res. Technol.. 7, (1 ), 1 –9 (2001). 0909-752X CrossRef
Pierce  M. C. et al., “Advances in optical coherence tomography imaging for dermatology,” J. Invest. Dermatol.. 123, (3 ), 458 –463 (2004). 0022-202X CrossRef
Drexler  W. et al., “Ultrahighresolution ophthalmic optical coherence tomography,” Nat Med.. 7, (4 ), 502 –507 (2001). 1078-8956 CrossRef
Fercher  A. F. et al., “In vivo optical coherence tomography,” Am. J. Ophthalmol.. 116, (1 ), 113 –114 (1993). 0002-9394 
Lima  I., Kalra  A., Sherif  S., “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express. 2, (5 ), 1069 –1081 (2011). 2156-7085 CrossRef
Lima  I. et al., “Fast calculation of multipath diffusive reflectance in optical coherence tomography,” Biomed. Opt. Express. 3, (4 ), 692 –700 (2012). 2156-7085 CrossRef
Yao  G., Wang  L., “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol.. 44, (9 ), 2307 –2320 (1999). 0031-9155 CrossRef
Tycho  A. et al., “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt.. 41, (31 ), 6676 –6691 (2002). 0003-6935 CrossRef
Wang  R. K., “Signal degradation by multiple scattering in optical coherence tomography of dense tissue: a Monte Carlo study towards optical clearing of biotissues,” Phys. Med. Biol.. 47, (13 ), 2281 –2299 (2002). 0031-9155 CrossRef
Kirillin  M. Yu. et al., “Monte Carlo simulation of optical clearing of paper in optical coherence tomography,” Quantum Electron.. 36, (2 ), 174 –180 (2006). 1063-7818 CrossRef
Wang  L., Jacques  S. L., Zheng  L., “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed.. 47, (2 ), 131 –146 (1995). 0169-2607 CrossRef
Kirillin  M. Yu. et al., “Visualization of paper structure by optical coherence tomography: Monte Carlo simulations and experimental study,” J. Europ. Opt. Soc. Rap. Public.. 2, , 07031  (2007). 1990-2573 CrossRef
Shi  B. et al., “Monte Carlo modeling of human tooth optical coherence tomography imaging,” J. Opt.. 15, (7 ), 075304  (2013). 0150-536X CrossRef
Kirillin  M. et al., “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express. 18, (21 ), 21714 –21724 (2010). 1094-4087 CrossRef
Meglinski  I. et al., “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett.. 33, (14 ), 1581 –1583 (2008). 0146-9592 CrossRef
Pfefer  T. et al., “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron.. 2, (4 ), 934 –942 (1996). 1077-260X CrossRef
Binzoni  T. et al., “Light transport in tissue by 3D Monte Carlo: influence of boundary voxelization,” Comput. Methods Programs Biomed.. 89, (1 ), 14 –23 (2008). 0169-2607 CrossRef
Boas  D. et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express. 10, (3 ), 159 –170 (2002). 1094-4087 CrossRef
Li  T., “MCVM: Monte Carlo modeling of photon migration in voxelized media,” J. Innov. Opt. Health Sci.. 3, (2 ), 91 –102 (2010). 1793-7205 CrossRef
Li  H. et al., “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method,” Acad. Radiol.. 11, (9 ), 1029 –1038 (2004). 1076-6332 CrossRef
Cote  D., Vitkin  I., “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express. 13, (1 ), 148 –163 (2005). 1094-4087 CrossRef
Margallo-Balbs  E., French  P. J., “Shape based Monte Carlo code for light transport in complex heterogeneous tissues,” Opt. Express. 15, (21 ), 14086 –14098 (2007). 1094-4087 CrossRef
Fang  Q., “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express. 1, (1 ), 165 –175 (2010). 2156-7085 CrossRef
Shen  H., Wang  G., “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol.. 55, (4 ), 947 –962 (2010). 0031-9155 CrossRef
Zhu  C., Liu  Q., “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt.. 18, (5 ), 050902  (2013). 1083-3668 CrossRef
Yadlowsky  M. J., Schmitt  J. M., Bonner  R. F., “Multiple scattering in optical coherence microscopy,” Appl. Opt.. 34, (25 ), 5699 –5707 (1995). 0003-6935 CrossRef
Kirillin  M. Y., Priezzhev  A. V., Myllylä  R. A., “Role of multiple scattering in formation of OCT skin images,” Quantum Electron.. 38, (6 ), 570 –575 (2008). 1063-7818 CrossRef
The Gnu Project, “Gnu scientific library,” http://www.gnu.org/s/gsl/ (15  November 2013).
Schoberl  J., “NETGEN an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Visual. Sci.. 1, (1 ), 41 –52 (1997). 1432-9360 CrossRef
Alerstam  E., Svensson  T., Andersson-Engels  S., “Parallel computing with graphics processing units for high-speed Monte-Carlo simulation of photon migration,” J. Biomed. Opt.. 13, (6 ), 060504  (2008). 1083-3668 CrossRef
Fang  Q., Boas  D. A., “Monte-Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express. 17, (22 ), 20178 –20190 (2009). 1094-4087 CrossRef
Leung  T. S., Powell  S., “Fast Monte Carlo simulations of ultrasound-modulated light using a graphics processing unit,” J. Biomed. Opt.. 15, (5 ), 055007  (2010). 1083-3668 CrossRef
Kalra  A., Lima  I.  Jr., Sherif  S., “Almost instantaneous Monte Carlo calculation of optical coherence tomography signal using graphic processing unit,” in  Proceedings of IEEE Photonics Conference (IPC) 2013 ,  Bellevue, Washington , paper MA2.3 (2013).
Lo  W. C. Y. et al., “Hardware acceleration of a Monte Carlo simulation for photodynamic therapy treatment planning,” J. Biomed. Opt.. 14, (1 ), 014019  (2009). 1083-3668 CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

PubMed Articles
Driven Boson Sampling. Phys Rev Lett 2017;118(2):020502.
Advertisement
  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.