Particle-fixed Monte Carlo model for optical coherence tomography Guanglei Xiong , Ping Xue, Jigang Wu, Qin Miao, Rui Wang and Liang Ji Dept. of Physics and Key Laboratory of Atomic and Molecular Nano Sciences, MOE, Dept. of Automation and National Laboratory for Information Science and Technology, Tsinghua University, Beijing, 100084, P.R. China [email protected]

Abstract: Optical Coherence Tomography (OCT) is a new technique mainly used in biomedical imaging. We present here a Particle-Fixed Monte Carlo (PFMC) simulation for OCT signal. In the PFMC model, the scattering particles of the sample are assumed to be temporarily fixed and randomly distributed in the simulation of the backscattered light. An efficient partitioning scheme is proposed to speed up this simulation process. The new model explains the exponential decay signal at the interfaces of different media layers observed in OCT experimental measurements. © 2005 Optical Society of America OCIS codes: (170.4500) Optical coherence tomography; (290.7050) Turbid media

References and links 1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991). 2. Z. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical Doppler tomographic image of fluid flow velocity in highly scattering media,” Opt. Lett. 22, 64–66 (1997). 3. G. J. Tearney, S. A. Boppart, B. E. Bouma, M. E. Brezinski, N. J. Weissman, J. F. Southern, and J. G. Fujimoto, “Scanning single-mode fiber optic catheter-endoscope for optical coherence tomography,” Opt. Lett. 21, 543– 545 (1996). 4. G. J. Tearney, B. E. Bouma, and J. G. Fujimoto, “High speed phase- and group-delay scanning with a gratingbased phase control delay line,” Opt. Lett. 22, 1811–1813 (1997). 5. Y. Pan, R. Birngruber, J. Rosperich, and R. Engelhardt, “Low-coherence optical tomography in turbid tissue: theoretical analysis,” Appl. Opt. 34, 6564–6574 (1995). 6. J. M. Schmitt, and A. Knuttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14, 1231–1242 (1997). 7. D. J. Smithies, T. Lindmo, Z. Chen, J. S. Nelson, and T. E. Milner, “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol. 43, 3025–3044 (1998). 8. G. Yao, and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44, 2307–2320 (1999). 9. A. Tycho, T. M. Jørgensen, H. T. Yura, and P. E. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt. 41, 6676–6691 (2002). 10. P. E. Andersen, L. Thrane, H. T. Yura, A. Tycho, T. M. Jørgensen, and M. H. Frosz, “Advanced modelling of optical coherence tomography systems,” Phys. Med. Biol. 49, 1307–1327 (2004). 11. Q. Lu, X. Gan, M. Gu, and Q. Luo, “Monte Carlo modeling of optical coherence tomography imaging through turbid media,” Appl. Opt. 43, 1628–1637 (2004). 12. Y. Chen, P. Xue, T. Yuan, W. Chen, and D. Chen, “Simulation of light scattering in optical coherence tomography,” Acta Optica Sinica 19, 486–490 (1999). 13. B. C. Karamata, P. Lambelet, M. Leutenegger, M. Laubscher, S. Bourquin, and T. Lasser, “A semi-analytical model accounting for multiple scattering in optical coherence tomography,” in Coherence Domain Optical Methods and Optical Coherence Tomography in Biomedicine IX , V. V. Tuchin, J. A. Izatt and J. G. Fujimoto, eds. , Proc. SPIE 5690, 386–396 (2005).

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2182

14. M. J. C. V. Gemert, S. L. Jacques, H. J. C. M. Sterenborg, and W. M. Star, “Skin optics,” IEEE Trans. Biomed. Eng. 36, 1146–1154 (1989). 15. L. Wang, S. L. Jacques, and L. Zheng, “MCML — Monte Carlo modeling of light transport in multi-layered tissues,” Computer Methods and Programs in Biomedicine 47, 131–146 (1995). 16. H. J. V. Staveren, C. J. M. Moes, J. V. Marle, S. A. Prahl, and M. J. C. V. Gemert, “Light scattering in Intralipid10% in the wavelength range of 400–1100 nm,” Appl. Opt. 30, 4507–4514 (1991). 17. J. L. Bentley, “Multidimensional binary search trees used for associative searching,” Communications of the ACM 18, 509–517 (1975).

1.

Introduction

Optical Coherence Tomography (OCT) [1] is a new optical technique which permits crosssectional tomographic imaging in highly turbid media. OCT is analogous to ultrasound imaging except that it uses low coherence light instead of sound waves. Compared with conventional X-ray computed tomography, magnetic resonance imaging, and ultrasound imaging, OCT is non-invasive, safe, portable and of high resolution. Recently, OCT technology has experienced a rapid development. By combining with the Doppler principle [2], catheter-endoscope technology [3] and rapid-scanning optical delay line [4], the applications of OCT have been widely expanded. In contrast, relatively small number of papers on theoretical analysis have been reported though there still remain many fundamental problems to be solved. In the past decade, Monte Carlo (MC) technique has been widely used in OCT simulations [5, 6, 7, 8, 9, 11] and proven to be very helpful for better understanding of OCT imaging process as well as the improvement of OCT related instrumentation and data processing algorithms. Using this method, Pan et al. [5] established the relationship between the path-length-resolved reflectance and the OCT signal using linear system theory. Schmitt and Knuttel [6] described an OCT model based on Huygens-Fresnel diffraction optics. They split the OCT signal into the summation of single backscattered light (coherent) and multiple scattered light (partially coherent). The effect of multiple scattering on the formation of speckle patterns and the degradation of image contrast were demonstrated. Chen et al. [12] studied the distribution of single and multiple backscattered light around the light source and simulated the OCT images formed by different orders of backscattered light. Tycho et al. [9] presented a MC method for modeling OCT measurements of a diffusely reflecting discontinuity embedded in a scattering medium. Andersen et al. [10] reviewed a MC model for calculating the OCT signal based on extended Huygens-Fresnel principle. In these models [9, 10], the extended Huygens-Fresnel principle was used to show the analytical results consistent with other models. Lu et al. [11] described a MC technique with Mie theory and wavelet for simulating OCT imaging through homogeneous turbid media. Among all the above analyses, the conventional Monte Carlo model is used, which employs a continuous exponential distribution to sample the photon’s step-length in the tissue, for they assume the scattering particles locate randomly in the media and every photon “sees” distinct realization of the media. In these models [5, 8], as a result of the averaging, the typical pathlength-resolved reflectance is a continuous curve (see Fig. 1). From the viewpoint of optical information processing [5], the optical transfer function of OCT acts like a narrow bandpass filter typically centered at the high-frequency range and hence slow-varying continuum distribution of path-length does not contribute to the OCT signal. That is to say, due to the coherent gate of OCT, only the fast-varying components contribute to the OCT signal. So when simulating a single-layer phantom with conventional Monte Carlo model, one can only get the OCT signal at the surface (see Fig. 2(a)), where the path-length-resolved reflectance varies dramatically, and there will be no signal in the depth of the scattering medium as the path-length-resolved reflectance decays smoothly. In real experiments, however, an exponential tail penetrating into the sample depth is detected (see Fig. 2(b)). Obviously, there is a discrepancy between the sim#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2183

ulation and the real signal. As a matter of fact, the conventional model assumes the scattering particles to be randomly distributed in the turbid media and thus the sampled path-length of the signal photon is also random, i.e., random distribution of the scattering particle simply results in randomness of the path-length of the signal photon. For the sake of simplicity and less computing time the conventional MC model samples the path-length alone regardless of the location of the scattering particle. The path-length of the signal photon from the same scatter/particle may differ from one sampling to another, which equally means that the locations of the scattering particles vary during the sampling process and consequently the sampling is the average of different distributions of scattering particle. This actually leads to continuum distribution of the scattering particles and the photon path-length. But for real case, the scattering particles are discretely distributed. So a discrete scattering model may potentially account for the exponential tail in depth of medium observed in OCT measurement more accurately than continuum scattering models. 4

x 10 14

Path−length−resovled reflectance

12

10

8

6

4

2

0

0

0.2

0.4

0.6 0.8 Path−length (mm)

1

1.2

1.4

Fig. 1. The conventional simulation of the path-length-resolved reflectance of a semiinfinite medium: refractive index 1.37, absorption coefficient 0.08mm−1 , scattering coefficient 10.0mm−1 and anisotropy factor 0.7. The curve is almost continuous.

In this paper, a novel Monte Carlo model, which we called Particle-Fixed Monte Carlo (PFMC) model, is presented to explain the phenomenon of exponential tail. In this model, the scattering particles are considered to be fixed with a random distribution, so the step lengths of signal photon in the scattering medium are discretely distributed. We will show that this discrete distribution contributes to the convolution of the path-length-resolved reflectance and the coherent function of light source and thus lead to the exponential tail signal, which has been actually observed in OCT measurement but cannot be interpreted by conventional Monte Carlo models. Our model differs from the existing models in the handling of photon propagation in the medium. The conventional MC simulation does not take the locations of particles into account, it simply represents the particles/sample with some statistical numbers, i.e., the scattering coefficient and absorption coefficient, etc.. The photon then propagates with steplengths sampled out of a distribution determined by these coefficients, which are the statistical optical properties of the medium. To some extent, this simplification can interpret the average properties of the sample in OCT measurement, but it is fundamentally unable to explain the phenomena related to some specific or fixed distributions of scatters. For example, as we know that OCT is sensitive to interface, the interface of small particles should also contribute to OCT signal. However, since the photon path sampling method in conventional MC simulations implicitly assumes that the medium is continuously varying at each scattering event, the OCT #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2184

1

0.9

0.9

0.8

0.8

0.7

0.7 Experimental OCT signal

Conventional simulated OCT signal

1

0.6

0.5

0.4

0.6

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.2

0.4

0.6 0.8 Depth (mm)

1

1.2

1.4

0

0.2

(a)

0.4

0.6 0.8 Depth (mm)

1

1.2

1.4

(b)

Fig. 2. The OCT signal obtained by the conventional models and by experiment. (a) The conventional simulation of the OCT signal: coherence length 17µ m and center wavelength 850nm. The signal only occurs at the surface; (b) The OCT signal of Intralipid 10% obtained by experiment. An exponential tail penetrating into the media is observed. There is a obvious discrepancy between the conventional model and experiment.

signal contributed by the interface of scatters are averaged and washed out. This is why the conventional MC simulation fails to interpret the exponential tail observed in OCT measurement. To better represent the medium, we need to take the distribution of scatters into account and study its intrinsic contribution to the OCT signal. We suggest a new sampling method of the photon’s step-length, which is determined by the locations of local scatters rather than random variable sampled out of a distribution governed by the optical properties of the medium. To get the step-length, one needs to know the location of the scatter interacting with the photon. So we first realize the distribution of scatters and then actually calculate the step-lengths based on the interaction between the scatters and the photon. We will demonstrate that PFMC is better to simulate the OCT signal in depth of the medium. However, it also should be noted that the path-length-resolved reflectance obtained by conventional MC simulations may still be very useful to describe the average OCT signal. 2.

Theory

The OCT system is a low coherence interferometer as shown in Fig. 3. The light intensity collected by the detector is given by [5] Id (τ ) = [Es (t) + Er (t + τ )][Es (t) + Er (t + τ )]∗   ∞

= [

−∞



 ∞

Es (t, Ls )dLs + Er (t + τ )][

−∞



Es (t, Ls )dLs + Er (t + τ )]∗ 

(1)

where Er and Es are the optical fields coming back from the reference and sample arm respectively; Lr and Ls are the corresponding round-trip path-lengths of the reference and sample   arms; Es is defined as the path-length-resolved field density given by Es = ∂ E(t, Ls )/∂ Ls . The angular brackets denote a long-time average over t and τ is the time delay corresponding to the round-trip path-length difference ∆L = Lr − Ls between two arms. Because the integration term that represents the interference effects is independent of the reference scan characterized by the

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2185

electronics

PC

detector

DC SLD

beam splitter

reference mirror

sample

Fig. 3. Schematic of OCT experimental setup: SLD, super-luminescent diode; DC, dispersion compensator; PC, personal computer. The reference mirror is moving at a constant speed to produce interference modulation, which is demodulated by the electronics and digitized via A/D converter to the computer for processing and OCT imaging display.

delay τ , Eq. 1 can be simplified as Id (Lr ) = Is + Ir + 2(Ir Is )1/2 



 ∞ −∞

[R(Ls )]1/2 g(Lr − Ls )dLs

(2)

∞ ∞ Es (t, Ls )dLs −∞ Es∗ (t, Ls )dLs  and Ir = Er (t + τ )Er∗ (t + τ ). R(Ls ) = where Is =  −∞ [dIs (Ls )/dLs ]/Is is the path-length-resolved diffuse reflectance of the sample. g(τ ) = E(t)E ∗ (t + τ )/I is the self-coherence function of the light source, where I = E(t)E ∗ (t) is the incident irradiance of the light source. From the last term of Eq. 2, we can get the OCT signal 



I(Lr ) ∝

 ∞ −∞

[R(Ls )]1/2 g(Lr − Ls )dLs

= [R(Lr )]

1/2

(3)

⊗ g(Lr )

where ⊗ denotes convolution. Therefore, the OCT signal I(τ ) is the convolution of the pathlength-resolved reflectance of the sample [R(Lr )]1/2 and the self-coherence function of the light source g(Lr ). From the viewpoint of linear system theory, the OCT signal can be obtained after the path-length-resolved reflectance passes through a filter with a optical transfer function as G(k), which is the Fourier transform of g(Lr ) equivalent to the power spectrum of the light source. For a Gaussian profile, G(k) can be written as G(k) ∝ exp[−(k − k0 )2 (Lc /4)2 ]

(4)

where k0 = 2π /λ0 and λ0 is the center wavelength, Lc is the coherence length. For a common SLD source with λ0 = 850nm and Lc = 17µ m, it can be easily shown that k0 ≈ 7.4µ m−1 and #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2186

∆k = 4 ln 2/Lc ≈ 0.16µ m−1 . G(k) is actually a narrow band-pass filter with a high center frequency k0 and a bandwidth ∆k. Therefore, it is sensitive only to the high-frequency component of [R(Ls )]1/2 within a path-length difference around λ0 . 3.

The PFMC model

From the above theory, the OCT signal can be expected by simulating the [R(Ls )]1/2 and then convoluting it with the assumed self-coherence function g(Lr ). Monte Carlo method can be used to simulate the photon transport in the sample and to obtain the path-length-resolved photon number which contributes to the OCT signal. A photon’s history consists of moving, absorption, and scattering. Conventional MC models use the continuous exponential distribution to sample the step-length d of photon moving: p(d) =

1 exp (−µt d) µt

(5)

where µt is the total attenuation coefficient, which is the sum of the absorption coefficient µa and the scattering coefficient µs . Once the photon has taken a step, its weight attenuates due to absorption, and the attenuation is equal to µa (6) ∆W = W µt After the photon has been moved and its weight decremented, it is ready to be scattered. The Henyey-Greenstein function [14] is used to get the deflected angle, θ . The probability distribution of cos θ can be described as p(cos θ ) =

1 − g2 2(1 + g2 − 2g cos θ )3/2

(7)

where g is the anisotropy factor. For a spherical particle, the azimuthal angle Ψ is sampled uniformly over the interval 0 to 2π . Assume that the phase shift caused by each scattering event is constant, the path-lengthresolved photon number is then proportional to R(Ls ). So we can easily obtain [R(Ls )]1/2 . Here we implicitly assume that the loss of coherence and depolarization effects due to scattering can be neglected. All photons that match the detection scheme contribute to R(Ls ) regardless of the number of scattering events. This assumption is valid because in OCT only a single, timevarying speckle is monitored, which by definition is spatially coherent [13]. In other words, the OCT detection scheme ensures the validation of above assumptions. The reflectance curve obtained by conventional MC models is varying dramatically only at the interface of two adjacent mediums and rather smooth elsewhere. Therefore, after the convolution, peaks only occur at interfaces in the simulated OCT signal and the exponential tail seen from experimental result does not appear. The PFMC model, which assumes the scattering particles in the sample are fixed during imaging process, can resolve this discrepancy. The physical image is clear: when the scattering particles are fixed, unlike the conventional models which involve sampling a continuous distribution, averaging of vicinal scattering particles cannot mask the fast-varying signal. In addition, the simulation results demonstrate that the path-length-resolved reflectance obtained by PFMC is not continuous especially at small path-length locations. Therefore, after convolution and demodulation, the obtained OCT signal stretches with an exponential tail penetrating to the depth of the scattering medium and matches well with the experiment. In the conventional models, each photon’s path-length is sampled from a continuous distribution. That is to say, every photon “sees” distinct realization of the medium. The dependency of #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2187

individual scattering particles is averaged out from the backscattered reflectance. The PFMC model differs from conventional MC on the sampling of photon step-length in the scattering medium. It proposes the scattering particles are fixed and every photon “sees” the same realization of the medium. We argue that this case is more close to real situations in tissue due to the short period of one A-scan process. Although, strictly speaking, the step-length distribution is also an exponential one as conventional MC, the discrete property of step-length caused by the discrete distribution of scattering particles must be taken into account, when the discrete scale is compared with the longitudinal resolution of the system. The PFMC model first samples the location of scattering particles, then computes the photon’s step-length by considering the interaction between the photon and scattering particles. The interaction distance is determined by the average particle size in the medium and it is related to the attenuation factor and particle density of the sample. The particle distribution pattern considered in our model is illustrated in Fig. 4. The location of particles is randomly sampled once and fixed throughout the simulation process. Every particle has an interaction sphere around it, shown as a dotted circle. It should be noted that although the particles in Fig. 4 is 2-D case, the PFMC model treat them in 3-D in real simulation.



Fig. 4. The particle distribution pattern in the PFMC model. The location of particles is randomly sampled and fixed throughout the simulation process. Every particle has an interaction sphere around it, shown as a dotted circle. The upside red dot denotes the incident photon.

4.

Implementation

For the implementation, when treating the uniform randomly distributed particles, we face the difficulty of huge computation burden because at every step of a single photon, all the particles must be checked to see if one interacts with it. What is more, the more steps the photon takes, the more computation the photon consumes, but the less its weight will contribute to the overall path-length-resolved reflectance. To overcome this dilemma we propose a fast algorithm. The basic idea of the algorithm is simple. Not all the particles need to be searched for whether they interact with the photon, that is, but only the photon’s neighboring ones are necessary. We suggest that the medium is partitioned by cubic lattice (as shown in Fig. 4). Though other complex partitioning schemes may be possible (like oct-tree or kd-tree [17]), they help little #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2188

due to the assumption of uniform particle distribution. The size of partition is chosen based on two aspects: if it is too large, searching in the partition still cost too much time; if it is too small, many partitions have no particles and thus they are useless. Therefore, appropriate partition size is a compromise between these two extremes (we find several tens of particles in each partition is sufficient.). The particles in each partition are stored at the beginning of the simulation. In the lifetime of photon, we record its current partition and events of entering and leaving partitions. In this way, Potential scattering particles can be only searched in the current partition. When the photon’s trajectory cuts through a particle’s sphere of radius interaction distance, the particle scatters the photon. If more than one particle is on the path of the photon, the nearest particle is chosen to scatter the photon (Note there is a rare case that more than one particles departs from the photon’s trajectory with the same distance smaller than their interaction distance. For this, the interaction can be regarded to cancel each other, and hence the photon is unchanged.). If no such an event, The photon will go into another adjacent partition or emit from the medium. Obviously, the ratio of the computation complexity without and with the partition scheme is Na /Np , where Na and Np are the particles in the medium and in a single partition respectively. This value is usually on the scale of 106 in our simulation. It is the partition scheme that makes the interaction between photon and randomly distributed particles efficient and practical. The ratio of the computation complexity between our algorithm with the partition scheme and the conventional MC simulation algorithm is Tp /Tc , where Tp and Tc are the time cost for the search of the next interacting particle and the sampling of the continuous exponential distribution. This value is usually on the scale of 102 in our simulation. The typical computing time is 3 ∼ 4 hours per 108 photons on an ordinary PC and this value reduces linearly with the number of nodes on cluster computers. We summarize the PFMC model and illustrate the flow chart in Fig. 5. 1. Launch photon : A new photon is initialized. Its path-length is set to zero (L = 0) and the weight to 1 (W = 1). A Gaussian beam focusing at the sample surface is simulated. The location of photon is sampled out of a Gaussian distribution determined by the radius of beam focus, which is 10µ m and the directional cosine set to (0, 0, 1). The current partition is set according to the partitioning scheme. As the waist of the Gaussian beam is on the sample surface, each photon is launched normal to the surface. In fact, the direction of photon is randomized after one scattering in the simulation. 2. Interact? : All particles in the local partition are checked to see if the photon’s trajectory can interact them. The particle is regarded as a small sphere with radius r. It is easy to deduce r = (µs /πρ )1/2 , where ρ is the particle density. We define r is the interaction distance (typically 1µ m). If interaction happens, go to (3); otherwise go to (5). 3. Move, Absorb and Scatter : The photon steps forward to the interacting particle, then its weight is decremented and its new direction is calculated as in Eq. (6) and Eq. (7). 4. Photon die? : A threshold value of the termination weight (e.g., Wth = 0.001) is set. After the photon’s weight falls below it, further propagation of the photon is supposed to yield little information about final result. A technique called roulette [15] of a dieing scheme is utilized to ensure the conservation of energy. If the photon’s weight is high, it is still alive and go to (2). 5. To partition boundary : The photon reaches the partition boundary and its path-length is updated. 6. Hit interface? : There are two situations when a photon hits the partition boundary. One is the photon hits tissue-air or tissue-tissue interface, and the other is just hitting a boundary #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2189

Begin

Launch photon

(1)

N

Interacts? (2)

(9)

(5)

To partition boundary

Update current partition

Y (3)

N

Hit interface? (6)

Move

Y R

Absorb

Transmit (7) or Reflect?

T Scatter

Emit?

(8)

N

Y N

N

Photon die? (4)

Collect?(10)

Y

N

Y Record

(11)

Last photon? (12)

Y End

Fig. 5. The flow chart of PFMC model

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2190

within tissue (Note that a tissue may have several partition layers). 7. Transmit or Reflect? : The internal reflectance r(αi ) is calculated by Fresnel’s formulas in statistical form:   1 sin2 (αi − αt ) tan2 (αi − αt ) + r(αi ) = (8) 2 sin2 (αi + αt ) tan2 (αi + αt ) where angle of incidence αi and the angle of transmission αt determined by Snell’s law. We then generate a random number, ξ . If ξ < r(αi ), the photon is internally reflected, go to (4); otherwise transmitted and go to (8). 8. Emit? : After transmitting, the photon may emit from tissue or still propagate within it, but in another layer when dealing with multi-layer tissue. 9. Update current partition : when the photon travels to the boundary of current partition, there are 26 possible adjacent partitions to jump into according to its direction. The current partition is updated and interaction within the new partition is checked as (2). 10. Collect? : Only a small amount of backscattered light can be detected by OCT. The geometry of detecting optics is employed to record the photon as shown in Fig. 6. 11. Last photon? : When a predefined number of simulation photons is traced, the simulation ends. This predefined number can be determined by the light power and the detection time (we choose 109 because the light power at sample arm is about 1mw). The results are compared after every 108 photons to make sure the irregular behavior of the pathlength-resolved reflectance stabilized.

φ0

2l0

tissue

Fig. 6. The light transport process and detection geometry of the PFMC model. The shaded area represents the simulating tissue. Dots in tissue represent the scattering particles and dashed circle around it marks the interaction region. The escaping photons within the restriction: φ < φ0 , l < l0 are collected and thus contribute to the calculation of path-lengthresolved reflectance according to the experiment setup in Fig. 3. Two example trajectories of photons, which are supposed to be collected, are also shown. The typical values are φ0 = 15◦ and l0 = 10µ m.

5.

Simulation and experiment

The followings are the results of simulating a semi-infinite scattering medium by PFMC model. The simulation gives the path-length-resolved reflectance [R(Ls )]1/2 . The OCT signal can be obtained by the convolution of [R(Ls )]1/2 with the self-coherence function g(Lr ). This OCT #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2191

signal is then obtained by quadrature demodulation. To compare with the conventional MC model, the mean free path (mfp) of the photon transport is calculated by averaging the steplengths of the photons when moving in the scattering medium. Figure 7(a) shows the simulation result of a medium whose parameters are designed to fit the real testing medium, Intralipid [16]. The sample thickness is d = 5mm, which is much larger than the penetrating depth of OCT and can be treated as semi-infinite. Other parameters are: the particle density ρ = 4.4 × 106 /mm3 , albedo µs /µt = 0.98, anisotropy g = 0.7, the interaction distance of particle r = 0.001mm and refractive index n = 1.1. The simulating mfp of the photon is 0.069mm, which corresponds to µt = 14.5mm−1 . The path-length grid size to accumulate the path-length of each photon is 0.0002mm. 5

10 3.5 4

10

3

3

Path−length−resovled reflectance

Path−length−resolved reflectance

2.5

10

2

2

1.5 0.2

10

0.3

0.4

0.5

1

10

4

10

3

10

0

10

2

−1

10

10

0

0.2

0.4

0.6

0.8 1 1.2 Path−length (mm)

(a)

1.4

1.6

1.8

2

0

0.2

0.4

0.6

0.8 1 1.2 Path−length (mm)

1.4

1.6

1.8

2

(b)

Fig. 7. Path-length-resolved reflectance simulated by the PFMC model and by conventional model. (a) The simulation of a semi-infinite medium: n = 1.1, d = 5mm, ρ = 4.4 × 106 /mm3 , µs /µt = 0.98, g = 0.7, r = 0.001mm by the PFMC model. The simulating mean free path is 0.069mm, corresponding to µt = 14.5mm−1 . The upper-right inset is the detail curve with path-length between 0.2 and 0.5mm; (b) The simulation of a semiinfinite medium: n = 1.1, d = 5mm, µs /µt = 0.98, g = 0.7, µt = 14.5mm−1 by conventional model.

We can see from the Fig. 7(a) that the curve simulated by PFMC is obviously not smooth, but varies sharply in the medium. These discontinuities come from the proposed discrete nature of particles by the model. Dominant photon paths lead to steps and impulses shown in the Fig. 7(a). Because of the scattering and absorbing nature of the medium, the intensity of the backscattered light attenuates. After convolution with g(Lr ), discontinuity (steps and impulses, see inset of Fig. 7(a)) can contribute to the OCT signal and the exponential tail will appear. For comparison, we also simulate by using the same parameters with the conventional model (see Fig. 7(b)). The curve is rather smooth and will yield no signal within the medium after convolution. The convolution of the curve by PFMC is demonstrated in Fig. 8(a) where the path-length is halved to get the longitudinal position. Figure 8(b) shows the demodulation curve. It can be seen that the exponential tail appears. Figure 9 compares the envelope of simulation and experiment result with Intralipid-10%. In the experiment, we used a standard OCT system (see Fig. 3) based on a linear translating stage. The light from SLD has the central wavelength 850nm and coherence length 17µ m. An axial scan of the reference arm is performed to get the depth-resolved one-dimension OCT signal. The scanning speed is 5mm/s and the modulation frequency is about 12kHz. The lens of the #6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2192

1200 1000

1000

Demodulation signal

Convolution signal

500

0

800

600

400

−500

200 −1000

0

0.1

0.2

0.3

0.4

0.5 Depth (mm)

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Depth (mm)

(a)

(b)

Fig. 8. The convolution and demodulation OCT signal of Fig.(a). (a) The convolution signal between path-length-resolved reflectance and the self-coherence function. It contains an envelope modulated by a carrier frequency; (b) The OCT signal obtained by quadrature demodulation of the convolution signal.

sample arm has a working distance of 22mm and outer diameter 7.2mm. Its focus is located on the air-sample interface. By exponential fitting of the curve with depth larger than 0, we can get the decay coefficients of simulation and experiment results, which are 3.0mm−1 and 3.4mm−1 , respectively. They agree well enough to prove that there is indeed an OCT signal inside the medium and the OCT signal attenuates exponentially. It should be noted that the decay coefficient is different from the attenuation coefficient of the sample. It reveals the decay behavior of the fast-varying components of the path-length-resolved reflectance. Small peaks and impulses are both observed in the simulation and experiment. They are caused by discrete nature of scattering particles. Obviously, if we can precisely know where the particles are located in the experimental tissue, the simulation and experiment will match even better at these peaks. 1 simulation experiment

0.9

Normalized OCT signal

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4 0.6 Depth (mm)

0.8

1

Fig. 9. Comparison of simulation and experiment results. The red dashed curve is the OCT signal of Intralipid-10% and The blue solid curve shows the PFMC simulation OCT signal. The maximum of both curved are normalized to 1 for easy comparison.

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2193

To show the usefulness of the PFMC model on tissue with multiple layers, we also simulate the path-length-resolved reflectance of a three-layer medium. The result is shown in Fig. 10(a). All the medium parameters of the three layers are the same except refractive index. The refractive indices are n1 = n3 = 1.1, n2 = 1.4, where subscripts 1−3 denote layer index. The thicknesses are d1 = d2 = 0.1mm, d3 = 1.5mm. Other parameters are: the particle density ρ = 5.0 × 106 /mm3 , albedo µs /µt = 0.98, anisotropy g = 0.7 and the interaction distance of particle r = 0.001mm. It can be seen that in each interface of two layers, there is a corresponding peak in the OCT signal. The peak locations 0.23mm and 0.51mm in the Fig. 10(a) agree well with the theoretical round-trip path-lengths of each interface: L1 = 2n1 d1 = 0.22mm, L2 = 2(n1 d1 +n2 d2 ) = 0.50mm. Also the intensity of photon number decreases rapidly in the first layer, and there are lots of discontinuities along this trend. They come from the discrete distribution of scattering particles. After convolution with the self-coherence function of optical field, these discontinuities will be magnified in the final OCT signal. Also it can be seen that in the second layer, the curve is still not continuous, which will result in an exponential tail penetrate into the third layer in the OCT signal. However, the curve beyond the last interface does not show large peaks because the photons, which can reach there and return, have experienced so many times of scattering which can mask the discrete nature of scattering particles. The small variation of the curve comes from the small number of backscattered photons. The result shows that the PFMC model can potentially account for the high sensitivity of OCT to sharp variations of medium parameters and explain the exponential tail of the OCT signal observed in the experiment. For comparison, we also simulate the same medium using the conventional model as shown in Fig. 10(b). As the simulation of the single layer, the curve is rather smooth except at the interfaces. It is also shown that as the number of layers increases, the PFMC and conventional MC come to be consistent. The difference between the results by the PFMC and the conventional model can strengthen the idea that the discrete nature of scattering particles in the medium can contribute to the OCT signal, especially near the surface of sample. first interface

first interface discontinuities

4

10

4

3

10

third interface

3

10

third interface

2

10

2

10

0.23mm

0.51mm

0

0.1

0.2

0.3

0.4 0.5 0.6 Path−length (mm)

(a)

0.23mm

1

10

1

10

second interface

second interface Path−lengh−resolved reflectance

Path−length−resolved reflectance

10

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.51mm 0.4

0.5

0.6

0.7

0.8

0.9

1

Path−lengh (mm)

(b)

Fig. 10. Path-length-resolved reflectance simulated by the PFMC model and by conventional model. (a) The simulation of a three layer medium with different refractive indices: n1 = n3 = 1.1, n2 = 1.4, d1 = d2 = 0.1mm, d3 = 1.5mm, ρ = 5.0×106 /mm3 , µs /µt = 0.98, g = 0.7, r = 0.001mm by PFMC; (b) The simulation of the same medium by conventional model (corresponding µt = 15.0mm−1 ). Both of the simulation results are not continuous at the interfaces, but only the result by PFMC shows discontinuities between interfaces.

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2194

6.

Conclusion and discussion

The Particle-Fixed Monte Carlo (PFMC) model is used to simulate OCT signal of scattering medium. Unlike the conventional models, it supposes the scattering particles are fixed in the process of imaging. Therefore, every photon “sees” the same realization of the medium. Fastvarying components due to the discreteness of path-length contribute to the OCT signal. The simulation result can explain the exponential tail not shown in conventional models but observed in experiment. In addition, a three-layer medium is simulated to show the effectiveness of the PFMC model on multiple-layered tissue. For the implementation, the partition scheme of the particles greatly reduce the computation cost of the interaction between photons and randomly distributed particles and makes PFMC practical. We also compare our algorithm with the conventional one. Considering the more accurate of our algorithm over the conventional one, the additional cost caused by the search of the next interacting particle is desirable and acceptable. It must be pointed out that, as the particles in the media may experience small movements around its location during the short imaging process, averaging of the result of several simulations, in which the locations of particles vary a bit, tends to agree with the real case more closely. Besides, if we can precisely know the locations of particles, the small peaks in Fig. 9 should match those of the experiment well. These will be investigated in our future work. Acknowledgments The authors thank Professor Lihong Wang at Texas A&M University for the helpful discussions. We are also grateful to Rowett Research Institute in Aberdeen, Scotland for allowing us to use its cluster and Dr. Tony Travis for the tip of programming on the cluster. We also thank the editor and referees for their valuable comments. This work was partially supported by THSJZ, Ministry of Science and Technology of China under contract No. 001CB510307 and National Science Foundation of China under grant No. 69908004.

#6559 - $15.00 US

(C) 2005 OSA

Received 8 February 2005; revised 13 March 2005; accepted 14 March 2005

21 March 2005 / Vol. 13, No. 6 / OPTICS EXPRESS 2195

Particle-fixed Monte Carlo model for optical coherence ...

Mar 21, 2005 - Y. Pan, R. Birngruber, J. Rosperich, and R. Engelhardt, ..... complex partitioning schemes may be possible (like oct-tree or kd-tree [17]), they help little ..... its cluster and Dr. Tony Travis for the tip of programming on the cluster.

384KB Sizes 3 Downloads 198 Views

Recommend Documents

Particle-fixed Monte Carlo model for optical ... - OSA Publishing
Mar 21, 2005 - tissues,” Computer Methods and Programs in Biomedicine 47, .... simply represents the particles/sample with some statistical numbers, i.e., the.

Monte Carlo simulations for a model of amphiphiles ...
Available online at www.sciencedirect.com. Physica A 319 (2003) .... many sites. The amphiphiles possess internal degrees of freedom with n different states.

Hamiltonian Monte Carlo for Hierarchical Models
Dec 3, 2013 - eigenvalues, which encode the direction and magnitudes of the local deviation from isotropy. data, latent mean µ set to zero, and a log-normal ...

Monte Carlo Simulation
You are going to use simulation elsewhere in the .... If we use Monte Carlo simulation to price a. European ...... Do not put all of your “business logic” in your GUI.

a monte carlo study
Mar 22, 2005 - We confirm this result using simulated data for a wide range of specifications by ...... Federal Reserve Bank of Kansas City and University of Missouri. ... Clements M.P., Krolzig H.$M. (1998), lA Comparison of the Forecast ...

Sequential Monte Carlo multiple testing
Oct 13, 2011 - can be reproduced through a Galaxy Pages document at: ... Then, in Section 3, we show on both simulated and real data that this method can ...

Statistical Modeling for Monte Carlo Simulation using Hspice - CiteSeerX
To enable Monte Carlo methods, a statistical model is needed. This is a model ..... However, it is difficult to determine the correlation without a lot of statistical data. The best case .... [3] HSPICE Simulation and Analysis User Guide. March 2005.

Using the Direct Simulation Monte Carlo Approach for ...
The viability of using the Direct Simulation Monte Carlo (DSMC) approach to study the blast-impact ... by computing load definition for two model geometries - a box and an 'I' shaped beam. ... On the other hand, particle methods do not make the conti

A Non-Resampling Sequential Monte Carlo Detector for ... - IEEE Xplore
SMC methods are traditionally built on the techniques of sequential importance sampling (SIS) and resampling. In this paper, we apply the SMC methodology.

accelerated monte carlo for kullback-leibler divergence ...
When using ˜Dts a (x) with antithetical variates, errors in the odd-order terms cancel, significantly improving efficiency. 9. VARIATIONAL IMPORTANCE ...

Introduction to Monte Carlo Simulation
Crystal Ball Global Business Unit ... Simulation is the application of models to predict future outcomes ... As an experimenter increases the number of cases to.

Sequential Monte Carlo multiple testing
Oct 13, 2011 - An example of such a local analysis is the study of how the relation ... and then perform a statistical test of a null hypothesis H0 versus. ∗To whom ... resampling risk (Gandy, 2009), and prediction of P-values using. Random ...

A Monte Carlo study of the Bell–Lavis model for water
mean-field solution. ..... However, on a global look, the Husimi cactus solution40 ..... 42 K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical.

ii SPECTROSCOPIC OPTICAL COHERENCE ...
Spectroscopic optical coherence tomography (SOCT) is a recent functional ...... broadband optical Gaussian beam; then for most cases the incident wave can be ..... constructed to alter the intensity of backscattered light from specific locations.

Fundamentals of the Monte Carlo method for neutral ...
Sep 17, 2001 - Fax: 734 763 4540 email: [email protected] cс 1998—2001 Alex F .... 11.3 Convergence of Monte Carlo solutions . . . . . . . . . . . . . . . . . . . . . .

Monte Carlo Null Models for Genomic Data - Project Euclid
To the rescue comes Monte Carlo testing, which may ap- pear deceptively ... order the null models by increasing preservation of the data characteristics, and we ...

Introduction to Monte Carlo Simulation - PDFKUL.COM
Monte Carlo Simulation Steps. • Following are the important steps for Monte Carlo simulation: 1. Deterministic model generation. 2. Input distribution identification. 3. Random number generation. 4. Analysis and decision making ..... perform output

A quasi-Monte Carlo method for computing areas ... - Semantic Scholar
Our method operates directly on the point cloud without any surface ... by counting the number of intersection points between the point cloud and a set of.

Copy of CMG14 monte-carlo methodology for network capacity ...
Quality of Service (QoS) = a generic term describing the performance of a system ... We have: 1. A network a. Topology, including Possible Points of Failure b.

Monte Carlo Simulation for the Structure of Polyolefins ...
unsaturated (usually vinyl) groups can be incorporated by the. CGC into a ... broaden the molecular weight distribution, since chains formed at the second site will .... to published experimental data for a lab-scale synthesis of a dual catalyst ...

Monte Carlo null models for genomic data - Semantic Scholar
Section 3 presents null model preservation hierarchies and signifi- cance orderings. Sections 4–6 ... A Galaxy Pages (Goecks et al., 2010) document allowing for ...