## 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

## 1. Introduction

Optical Coherence Tomography (OCT) [1] is a new optical technique which permits cross-sectional 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 path-length-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 simulation 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.

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 step-lengths 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 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]

$$=\u3008\left[{\int}_{-\infty}^{\infty}{E\prime}_{s}(t,{L}_{s})d{L}_{s}+{E}_{r}\left(t+\tau \right)\right]{\left[{\int}_{-\infty}^{\infty}{E\prime}_{s}(t,{L}_{s})d{L}_{s}+{E}_{r}\left(t+\tau \right)\right]}^{*}\u3009$$

where *E _{r}* and

*E*are the optical fields coming back from the reference and sample arm respectively;

_{s}*L*and

_{r}*L*are the corresponding round-trip path-lengths of the reference and sample arms;

_{s}*E′*is defined as the path-length-resolved field density given by

_{s}*E′*=

_{s}*∂E*(

*t,L*)/

_{s}*∂L*. The angular brackets denote a long-time average over

_{s}*t*and

*τ*is the time delay corresponding to the round-trip path-length difference Δ

*L=L*between two arms. Because the integration term that represents the interference effects is independent of the reference scan characterized by the delay

_{r}-L_{s}*τ*, Eq. 1 can be simplified as

where *I _{s}*=〈∫∞-∞

*E′*(

_{s}*t,L*)

_{s}*dL*∫∞-∞

_{s}*E′**(

_{s}*t,L*)

_{s}*dL*〉 and

_{s}*I*=〈

_{r}*E*(

_{r}*t+τ*)

*E**(

_{r}*t+τ*)〉.

*R*(

*L*)=[

_{s}*dI*(

_{s}*L*)/

_{s}*dL*]/

_{s}*I*is the path-length-resolved diffuse reflectance of the sample.

_{s}*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

$$={\left[R\left({L}_{r}\right)\right]}^{1\u20442}\otimes g\left({L}_{r}\right)$$

where ⊗ denotes convolution. Therefore, the OCT signal *I*(*τ*) is the convolution of the path-length-resolved reflectance of the sample [*R*(*L _{r}*)]

^{1/2}and the self-coherence function of the light source

*g*(

*L*). 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

_{r}*G*(

*k*), which is the Fourier transform of

*g*(

*L*) equivalent to the power spectrum of the light source. For a Gaussian profile,

_{r}*G*(

*k*) can be written as

where *k*
_{0}=2*π/λ*
_{0} and λ_{0} is the center wavelength, *L _{c}* is the coherence length. For a common SLD source with λ

_{0}=850

*nm*and

*L*=17µm, it can be easily shown that

_{c}*k*

_{0}≈ 7.4

*µm*

^{-1}and Δ

*k*=4ln2/

*L*≈0.16

_{c}*µm*

^{-1}.

*G*(

*k*) is actually a narrow band-pass filter with a high center frequency

*k*

_{0}and a bandwidth Δk. Therefore, it is sensitive only to the high-frequency component of [

*R*(

*L*)]

_{s}^{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*(*L _{s}*)]1/2 and then convoluting it with the assumed self-coherence function

*g*(

*L*). 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.

_{r}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:

where *µ _{t}* is the total attenuation coefficient, which is the sum of the absorption coefficient

*µ*and the scattering coefficient

_{a}*µ*.

_{s}Once the photon has taken a step, its weight attenuates due to absorption, and the attenuation is equal to

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

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-length-resolved photon number is then proportional to *R*(*L _{s}*). So we can easily obtain [

*R*(

*L*)]

_{s}^{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*(

*L*) regardless of the number of scattering events. This assumption is valid because in OCT only a single, time-varying speckle is monitored, which by definition is spatially coherent [13]. In other words, the OCT detection scheme ensures the validation of above assumptions.

_{s}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 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.

## 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 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 *N _{a}/N_{p}*, where

*N*and

_{a}*N*are the particles in the medium and in a single partition respectively. This value is usually on the scale of 10

_{p}^{6}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

*T*, where

_{p}/T_{c}*T*and

_{p}*T*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 10

_{c}^{2}in our simulation. The typical computing time is 3~4 hours per 10

^{8}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 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:

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 10^{9} because the light power at sample arm is about 1*mw*). The results are compared after every 10^{8} photons to make sure the irregular behavior of the path-length-resolved reflectance stabilized.

## 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*(*L _{s}*)]

^{1/2}. The OCT signal can be obtained by the convolution of [

*R*(

*L*)]

_{s}^{1/2}with the self-coherence function

*g*(

*L*). This OCT 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 step-lengths of the photons when moving in the scattering medium.

_{r}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*=5*mm*, 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×10^{6}/*mm*
^{3}, albedo *µ _{s}/µ_{t}*=0.98, anisotropy

*g*=0.7, the interaction distance of particle

*r*=0.001

*mm*and refractive index

*n*=1.1. The simulating mfp of the photon is 0.069mm, which corresponds to

*µ*=14.5

_{t}*mm*

^{-1}. The path-length grid size to accumulate the path-length of each photon is 0.0002

*mm*.

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*(*L _{r}*), 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 850*nm* 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 5*mm/s* and the modulation frequency is about 12*kHz*. The lens of the sample arm has a working distance of 22*mm* and outer diameter 7.2*mm*. 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.0*mm*
^{-1} and 3.4*mm*
^{-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.

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 *n*
_{1}=*n*
_{3}=1.1,*n*
_{2}=1.4, where subscripts 1–3 denote layer index. The thicknesses are *d*
_{1}=*d*
_{2}=0.1*mm*,*d*
_{3}=1.5*mm*. Other parameters are: the particle density *ρ*=5.0×10^{6}/*mm*
^{3}, albedo *µ _{s}/µ_{t}*=0.98, anisotropy

*g*=0.7 and the interaction distance of particle

*r*=0.001

*mm*. It can be seen that in each interface of two layers, there is a corresponding peak in the OCT signal. The peak locations 0.23

*mm*and 0.51

*mm*in the Fig. 10(a) agree well with the theoretical round-trip path-lengths of each interface:

*L*

_{1}=2

*n*

_{1}

*d*

_{1}=0.22

*mm,L*

_{2}=2(

*n*

_{1}

*d*

_{1}+

*n*

_{2}

*d*

_{2})=0.50

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

## 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. Fast-varying 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 LihongWang at Texas A&MUniversity 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.

## 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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**4. **G. J. Tearney, B. E. Bouma, and J. G. Fujimoto, “High speed phase- and group-delay scanning with a grating-based phase control delay line,” Opt. Lett. **22**, 1811–1813 (1997). [CrossRef]

**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). [CrossRef] [PubMed]

**6. **J. M. Schmitt and A. Knuttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A **14**, 1231–1242 (1997). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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. SPIE5690, 386–396 (2005).

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**16. **H. J. V. Staveren, C. J. M. Moes, J. V. Marle, S. A. Prahl, and M. J. C. V. Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400–1100 nm,” Appl. Opt. **30**, 4507–4514 (1991). [CrossRef] [PubMed]

**17. **J. L. Bentley, “Multidimensional binary search trees used for associative searching,” Communications of the ACM **18**, 509–517 (1975). [CrossRef]