Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation

A Monte Carlo model has been developed for optical coherence tomography (OCT). A geometrical optics implementation of the OCT probe with low-coherence interferometric detection was combined with three-dimensional stochastic Monte Carlo modelling of photon propagation in the homogeneous sample medium. Optical properties of the sample were selected to simulate intralipid and blood, representing moderately ( g = 0.7) and highly ( g = 0.99) anisotropic scattering respectively. For shallow optical depths in simulated intralipid (<3 scattering mean free path (mfp) units), the number of detected backscattered photons followed the extinction-single-backscatter model, and OCT was found to detect only minimally scattered photons. Within this depth range the backscatter positions of detected photons corresponded well with the nominal focus position of the probe. For propagation to deeper positions in intralipid, localization of backscattering was quickly lost due to detection of stray photons, and the number of detected photons remained constant with increasing depth in the non-absorbing medium. For strongly forward-directed scattering in simulated blood, the number of detected photons approached the extinction-single-backscatter model only for very shallow depths (<2 mfp units). However, backscattering positions for detected photons correlated well with the nominal focus position of the probe even for optical depths greater than 40 mfp units.


Introduction
Optical coherence tomography (OCT) is a point-scanning method which permits imaging in highly turbid media and thus yields high-resolution three-dimensional backscatter information from depths down to 1-2 mm in biological tissues (reviewed by Fercher 1996). Analogous to ultrasound, the technique can provide information on local flow velocities in scattering media (Gusmeroli and Martinelli 1991). The term optical Doppler tomography (ODT) has been introduced to describe the combination of OCT and Doppler shift detection (Wang et al 1995a, b).
Signal attenuation and loss of resolution due to multiple scattering severely limit the depth penetration of OCT, and do not allow imaging deep into strongly scattering biological tissues. Mathematical models have been developed to promote understanding of the OCT imaging process and thereby enable development of better imaging instrumentation and data processing algorithms. The extinction-single-backscatter model represents the OCT imaging process as an incoming beam which is continuously attenuated according to the total extinction coefficient of the propagation medium (Schmitt et al 1993a, b, 1994b, c, Yadlowsky et al 1995a. After a single backscattering event, light is backpropagated towards the detector under continuous attenuation. Multiple scattering effects have been included into such analytical models (Yadlowsky et al 1995b), but the models are still of limited value for comparison with experiment since only simple geometries can be considered.
It is to a large extent still unclear how multiple scattering effects will influence OCT and ODT. Monte Carlo simulation provides a method for implementing realistic models of the experimental situation based on specified geometries, characteristic parameters of source and detector, as well as full stochastic modelling of light propagation in the scattering medium to be probed. Comparison between simulations and experimental measurements can then be used to test, validate and refine the Monte Carlo model or modify elements of the experimental set-up. Based on achieved agreement between model and measurements in a set of parallel simulations and experiments, the model can be used to predict experimental results for other instrument configurations, and thus lead the way to optimal instrument design.
We have developed a Monte Carlo model to study both OCT and ODT. Our motivation was to establish an exploratory model that could simulate essential aspects of the OCT and ODT measurement process. Although certain output functions could be specified before programming the model, we wanted to keep the model so general that new observation functions could be incorporated as our studies unveiled more and more of the phenomenology of OCT and ODT.
Full three-dimensional stochastic simulation of photon propagation in the medium was implemented. The OCT detector was modelled with low-coherence interferometric path length gating and confocal selectivity determined by the antenna theorem for a heterodyne receiver (Siegman 1966, Schmitt et al 1994b. By this theorem the effective size of the OCT detector is limited to twice the diameter of an ideal confocal detector. This represents nearly full confocality, and OCT has been said to perform as a confocal microscope with an additional narrow time gate (∼50 fs) to further increase the selectivity of photon detection (Yadlowsky et al 1995a).
In terms of Monte Carlo simulation, the high detection selectivity of OCT translates into very low yields of detected photons, and therefore inefficient and time-consuming computations. Although computing time might be reduced by transform methods such as similarity scaling (Stern 1985), we refrained from such techniques, so all data from the Monte Carlo simulation were available for direct interpretation. Simulation results included the number of detected photons, individual accumulative path lengths, and Doppler shifts (the latter for an accompanying study of Doppler flow measurements). Also included were records of the spatial locations for all scattering interactions, which were used to calculate a point spread function for backscattering relative to the nominal focus position of the probe.
Optical properties of the media were selected to mimic samples of biological relevance that could be subjected to OCT measurements in the laboratory. Intralipid solutions (2 and 5%) were modelled to represent moderately anisotropic scattering properties (g = 0.7). Furthermore, blood was modelled, not only for its biological relevance, but also because it represents one of the most anisotropic of biological scattering media (g = 0.99)

Materials and methods
2.1. Experimental OCT system 2.1.1. Instrumentation. A schematic diagram of the fibre-optic OCT instrument is shown in figure 1. Only a brief description will be given here since more detailed information has been published elsewhere (Chen et al 1997a, b). Continuous near-infrared light (λ = 850 nm, λ = 25 nm) emitted by the super-luminescent diode (SLD) is coupled into a fibre-optic Michelson interferometer and split into two beams by a 2 × 2 fibre coupler. The optical phase of the SLD light in the probe and the reference arm of the interferometer is modulated at 1600 Hz by stretching the fibre wrapped around piezoelectric cylinders driven by a ramp waveform. Stress birefringence is used to match the polarization of the probe and the reference beams and thus optimize fringe contrast.
A microlens of NA 0.2 terminates the fibre in the probe arm and focuses light into a 5 µm diameter spot at a working distance of 3 mm. As a detector, light from the focus point is imaged onto the 5 µm diameter fibre core, which effectively serves as the pinhole of a confocal detector. The degree of confocality can be expressed by the dimensionless parameter v p = 2π(NA)r p /λ (Schmitt and Ben-Letaief 1996) where r p is the radius of the detector, i.e. the fibre core. For our detector v p = 3.7, whereas 1.0 corresponds to full confocality. Light that is backscattered from the sample recombines in the 2 × 2 fibre coupler with light from the retroreflector in the reference arm. The two beams interfere and give fringes only when their optical pathlength difference is less than the coherence length, L c , of the SLD source light (FWHM = 0.44λ 2 /n λ, which is 13 µm in air). Because the coherence envelope of the SLD source light yields a rapid phase decorrelation of the beams for pathlength differences greater than the coherence length, good depth resolution is achieved.
Two-dimensional depth scan images are formed by sequential lateral scans, each followed by an incremental probe movement in the vertical direction. To maintain a zero optical pathlength difference between the beam focus in the sample and the reference mirror, we use dynamic focus tracking (Chen et al 1997a, b).
Optical interference fringe intensity is measured by a silicon photovoltaic receiver (New Focus 2001, New Focus Inc:, Santa Clara, CA). The signal is digitized with a 16-bit analogue-digital converter and transferred to a computer workstation for processing. The detected interference fringe intensity at each focus position is transformed by a short-time, fast Fourier transform algorithm (STFFT). A tomographic structural image is obtained by calculating the power in the Fourier component corresponding to the phase-modulation frequency for each focus position and displaying the value as a false-colour dB intensity value in the corresponding pixel of the image. (intralipid 20%,Kabi Pharmacia Inc.,Clayton,NC) were made by dilution in aqua dest., and poured into 5 cm diameter Petri dishes for imaging. Optical parameters for 1% intralipid at 850 nm were taken to be µ a = 0, µ s = 2.0 mm −1 and g = 0.7, as average values of those obtained from van Staveren et al (1991) and Flock et al (1992). Both authors state µ a ≈ 0, and their values for µ s were 2.3 mm −1 and 1.7 mm −1 , respectively; and 0.61 and 0.78 for g.

Samples. Different concentrations of intralipid
A suspension of glutaraldehyde-fixed canine erythrocytes in density-matched CsCl (42%) was used to represent human blood. This model was developed for our ODT studies (Lindmo et al 1998), where a stable, non-settling suspension is required to mimic the optical properties of human blood. Optical parameters for blood at 850 nm were adopted from Jacques (1996) and Yaroslavsky et al (1996). Both authors give values of µ a = 0.75 mm −1 , but report values for µ s of 300 mm −1 and 65 mm −1 respectively, and g-values of 0.996 and 0.99. We chose to use µ a = 0.75 mm −1 , µ s = 150 mm −1 and g = 0.99.

Monte Carlo simulation of OCT
2.2.1. Probe geometry. The geometry of the combined emission and detection probe is shown in figure 2. Since we also measure Doppler shifts (Lindmo et al 1998), the probe is pointed at an incidence angle α of 10 to 20 • . The probe has a numerical aperture of NA = sin θ, and the marginal rays BP and CP intersect the surface at incidence angles β = α + θ and λ = α − θ, respectively. Without refraction at the surface of the sample medium, light from the probe would be focused at point P at depth z. Due to an index of refraction n > 1 (n = 1.37 in our simulations), the in-air focus depth z will be shifted down to a depth in the medium z = k z z. The proportionality constant k z is determined from the geometry in figure 2, by expressing the length of AB in quantities relating to air as well as to the medium Primed variables refer to values after refraction into the medium, and angles in air are related to those in the medium by Snell's law, for example sin β = n sin β . Refraction of light at the surface also introduces a slight shift in the y-position: The probe is positioned so as to have its focus in the medium at P , at depth z in the xz-plane. The probe is segmented into a number of light-emitting elements of equal solid-angle. This is done by subdividing the spherical surface S e into one polar element surrounded by m r concentric rings, each divided into m a azimuthal sectors, resulting in a total of 1 + m r m a elements made to subtend equal solid angles. A pencil beam is launched along the central axis of each element, thus being directed towards the in-air focus point P . The number of photons in each pencil beam is usually chosen to correspond to a Gaussian intensity profile across the surface S e . Upon entry into the medium, each pencil beam refracts and changes its direction towards P .
The physical probe has a fixed dimension (i.e. P H = 3 mm), but for the Monte Carlo simulation, the apparent probe surfaces S e and S d , for emission and detection respectively, are defined as the spherical sector surface, centred at P with aperture angle θ , which will touch the medium surface at B. In this way the probe surface S will always have a size comparable in dimension to the focus depth z, and will therefore be easy to visualize on Figure 2. Schematic view of the combined emission and detection probe of the OCT/ODT instrument as modelled for the Monte Carlo simulations. A focusing lens H H 1 transforms the diverging beam from the fibre tip at P 1 into a beam converging to the point P in air. The spherical wavefront surface on the fibre tip side is S 1 , on the focus side S e or S d , depending on whether emission or detection is considered. The probe is aimed at an angle of incidence α and has numerical aperture angle θ(NA = sin θ). The core of the optical fibre subtends an angle of 2δ as seen from a point on the S 1 wavefront. Refraction at the surface of the medium causes the in-air focus point P to be shifted to P at depth z . the same scale as the ensemble of photon paths. Since r p scales proportionally to P 1 H 1 , the value of δ remains constant even if the apparent probe dimension is changed with depth z to keep point B on the surface of the sample medium. The diameter of the apparent probe is If the emitter and the detector probes are coincident and of equal numerical aperture, the two surfaces S e and S d will be coincident, as in figure 2. The probe lens H H 1 transforms the diverging spherical wavefront S d of backscattered photons into the converging wavefront S 1 . For any photon to be detected, it must couple into the fibre tip which subtends an angle 2δ ≈ 2r p /P 1 H 1 as seen from any point on S 1 , where 2r p is the diameter of the fibre core. The angle δ, which we call the confocality angle, will vary from r p /P 1 H 1 for axial rays to r p cos θ/P 1 H 1 for marginal rays. For our purposes, with typical values of θ = 15 • , we use the approximation cos θ = 1.0, thus neglecting the slight variation in δ with position on the surface S 1 .
Based on a selected emitter-probe configuration, photon histories are generated by Monte Carlo simulation and saved for later analysis using several different detector-probe configurations. However, to reduce storage requirements, the intended orientation of the detector probe is specified before the simulation (always being coincident with the emitter probe for the present work). The other detection parameters, such as NA, L c , and δ are set at suitably selected maximum values for the simulation in order to allow subsequent analysis using the generated photon histories with different values of detection parameters within the originally specified range.

Photon histories generator.
A selected number of photons are generated for each probe position, aimed at depth z at lateral position x, corresponding to the position (x, 0, z ) in the medium. The simulation may contain one or two outer program loops respectively, that generate an axial depth scan by a series of equal increments of the z -variable, or a full two-dimensional xz -scan made from a set of axial z -scans at different x-positions.
Each individual photon launched from the probe is propagated into and within the medium, while accounting for its accumulated pathlength. This parameter also includes the initial path segment in air from S e to the surface of the medium, divided by n to convert it into medium-equivalent distance units. Within the medium, the length of each path segment is determined by stochastically sampling the pathlength distribution according to the method of Buslenko et al (1966) as described previously (Smithies and Butler 1995). At each interaction within the medium, the scattering angle ϕ is determined by stochastically sampling the Henyey-Greenstein phase function (Henyey andGreenstein 1941, Cashwell andEverett 1959). For each interaction, the position of the scattering event is recorded, and the accumulated pathlength is incremented by the length of the latest path segment. For the generation of photon histories, maximum values for detection parameters such as source coherence length in the medium (L c ), numerical aperture (NA), and confocality angle (δ), are pre-set by the user. For each depth position the minimum and maximum pathlength, or coherence gate, for backscattered photons that may be detected, is determined as Propagation of any photon is terminated and its history discarded if its accumulated path length becomes longer than p max less the shortest distance from its present position back to the surface S d . Because backscattering events in highly anisotropic scattering media are very infrequent, increasing the yield of detectable photons just by increasing the number of input photons is not efficient. To take full computational advantage of a backscattering event, backscattered photons can be split into a number of new photons of correspondingly smaller weight, and individually propagated further in the simulation (Jenkins et al 1988). We chose to split photons only once, with a factor of 250 at the first scattering event which gives a photon direction deviating less than 60 • from the normal direction back to the surface.
Backscattered photons that emerge through the surface of the medium are refracted and propagated to the surface S d . If the photon intersects the surface within the maximum numerical aperture, NA max , detection is possible. The incidence angle of the photon with the surface S d is computed to test if the value is within the maximum confocality angle δ max . Path histories are written to a disk file for each photon that satisfies the set of maximum values for the detection parameters NA, L c and δ. Each history includes (x, y, z )-positions of all scattering events, spatial and angular coordinates of the photon's intersection with the detector surface S d , and the accumulated pathlength calculated from leaving S e and incremented all the way into the medium until reaching S d as a backscattered photon.

Photon histories processor.
The data processing program is implemented with the AVS image analysis software (Advanced Visualization System, Waltham, MA), and takes the form of a number of specially developed modules that integrate into the entire software package. The user sets module variables interactively, with on-screen controls, and can use standard AVS modules for further data processing and graphics presentation (Upson et al 1989).
In this report, results will be presented that were generated with the following modules written for the OCT Monte Carlo simulation. The data reader reads photon histories from disk file for a specified subset, or all, of the focus positions simulated in a specified run of the photon histories generator. The photon detector determines which of the read photons are detected, i.e. satisfy the detection criteria for NA, L c and δ. The photon counter registers the number of detected photons for each (x, z ) position in a depth scan. The scattering events distribution generates, for a specified focus position (x, z ), a histogram showing how many photons were detected as a function of the number of scattering events per photon. The path lines generator provides a graphical display of the ensemble of individual path segments for the detected photons from one (x, z ) focus position of the input probe. The paths are drawn from the emitter surface S e , throughout the medium and back to the detector surface S d . This module utilizes other standard AVS graphics modules that allow scaling and rotation of the three-dimensional graphics output. The point spread function depicts, for all detected photons from a certain focus position (x, z ), all the interaction points that the photons have had in the medium. Alternatively, the module displays only the interaction point at maximum depth for individual photons. These are the points of backscattering, and this output is defined as the point spread function. The geometry can be viewed as an xz -image with the data projected along the y-axis, or as a yz -image with the data projected along the x-axis. For each of these two images, projected one-dimensional distributions can be generated corresponding to the depth (zPSF) and lateral point spread functions (xPSF or yPSF respectively).

Implementation.
Simulations and subsequent analyses were run on a SUN Ultra 170E workstation (Sun Microsystems Inc., Mountain View, CA). Special-purpose and standard graphics software were utilized during testing and debugging of the simulation program in order to verify performance. The splitting technique was tested by comparing parallel runs with and without splitting. Good agreement was found as long as the detected split photons represented many different input photons, i.e. many different photon histories prior to splitting. Routinely, the number of independent input photon histories represented among the detected photons, was at least 30% of the number of detected photons.
Run times for a typical optical depth of 10 scattering mean free paths in blood was about 1 h CPU time for tracking 10 9 photons at one focus position. Data files requiring up to 100 Mb of storage space were generated, allowing typical simulations like depth scans to contain a series of 20 to 40 different focus positions.

Attenuation of signal
3.1.1. Experimental results. Attenuation of backscattered light from homogeneous samples was determined by recording two-dimensional depth scans from samples of intralipid or reconstituted blood contained in small Petri dishes. The images consisted of scans 400 µm wide in the x-direction taken at 10 µm intervals of depth in the z-direction. Each image was collapsed into a one-dimensional intensity depth profile by averaging over the x-dimension.
Backscatter intensity distributions for 2 and 5% intralipid are shown in figure 3A, plotted as a function of optical depth, µ s z , i.e. depth measured in number of scattering mean free path units (1 mfp = 1/µ s ). The ordinate is expressed in dB, as S[dB] = 10 log 10 (G(z )) where G(z ) is the square of the short-time fast Fourier transform (STFFT) of the detected temporal interference fringe intensity. The peak to the left in each curve represents the sharp increase in backscatter detected as photons originate from the scattering medium. After the peak, both curves demonstrate a straight segment of sharply decreasing signal with optical depth, followed by a less steep decrease for greater depths. Plotted as a function of optical depth, the results for 2% intralipid which were recorded down to a depth of 2 mm (i.e. optical depth of 4 mm −1 × 2 mm = 8), superimpose quite accurately on the curve for 5% intralipid, demonstrating that depth dependence is governed by the increase in optical depth. The straight broken line in the plot shows attenuation according to the extinction function exp(−2µ s z / cos α ). The factor 2/ cos α takes into consideration that the photon path in and out at an incidence angle of α is twice the distance to the vertical depth z . The agreement between the theoretical curve and the initial part of the experimental curve suggests that for shallow optical depths the initial slope of the OCT intensity distribution for intralipid may be adequately represented by the extinction model.
For optical depths greater than about 4, corresponding to a round-trip path of 8 mean free paths, the experimental curves for intralipid show considerably higher signal than predicted by the extinction model. Figure 3B shows corresponding data for undiluted and 4× diluted canine blood. In these curves the segments of logarithmic attenuation are less pronounced, but still indicate some agreement with the slope represented by the extinction model (straight broken line). The peak of the curve for undiluted blood (full curve) is broadened, probably due to limited depth resolution. (One unit of optical depth for blood is 1/µ s = 7 µm.) 3.1.2. Simulations. Monte Carlo simulations were performed for homogeneous samples of intralipid and reconstituted canine blood, using the same sample concentrations and scanning depths as for the experimental investigation. An input probe configuration of 15 • incidence angle, NA = 0.2 and a Gaussian beam intensity profile were chosen, and photon histories were generated for various sets of sample parameters.
The Monte Carlo simulation result that may most directly be compared with experimental backscattering intensity data is output from the photon counter module, which shows the number of detected photons as a function of depth in the medium. Such curves are shown in figure 4 as a function of optical depth (µ s z ). As in the experimental case, simulated data for two different media concentrations generate superimposable curves when plotted as a function of optical depth. This is understandable because the geometry of the simulation problem scales with the scattering mean free path. A simulation for a more dilute medium only increases the metric distances of the geometry by a factor equal to the ratio of scattering mean free paths. The absolute yield of photons is, however, proportional Figure 3. Experimental data from OCT imaging of intralipid (A) and reconstituted canine blood (B) at an angle of incidence of 10 • (=7.3 • in the medium). Signal intensity was averaged over the x-direction in OCT images and plotted as a function of optical depth (µ s z ). Values of the scattering coefficients were assumed to be 10 mm −1 and 4 mm −1 for 5% and 2% intralipid respectively (A); 150 mm −1 and 37.5 mm −1 for undiluted and 4× diluted blood respectively (B). The full curves represent the higher concentrations of intralipid and blood, respectively, and the dotted curves represent the lower concentrations. Pronounced peaks at significant depth in two of the curves are artefacts, probably due to spurious reflections in the probe optics. For each curve the measure of optical depth should have started with zero at the position of the peak representing the air/medium interface, but the curves were shifted along the x-axis to make the initial segments of logarithmic decrease superimpose. This representation facilitates a comparison with the extinction model indicated by the straight, broken lines depicting the function S = K − 10 × 2µ s z log 10 e/ cos(7.3 • ).
to the source coherence length since this parameter serves as a gate on the interval of photon pathlengths to accept for detection. Complete equivalence of the simulation for two different medium dilutions therefore requires equal source coherence lengths expressed in the respective scattering mean free paths. Recorded photon counts for more dilute samples were therefore multiplied by the ratio of scattering mean free paths for photons in the more dilute to the denser medium in order to generate superimposable curves such as shown in figure 4. For example, a factor of 2.5 was used for the 2% curve of intralipid to plot it with the 5% curve, and a factor of 4 was used for 4× diluted blood versus undiluted blood.
The number of detected photons was found to be approximately proportional to the product of the square of the detector numerical aperture and the source coherence length (results not shown). A change in either of these parameters gave a shift along the ordinate axis of data such as displayed in figure 4. To increase the number of detected photons, the plotted curves therefore represent the maximum values of NA = 0.4 and L c = 50 µm.
The curves for intralipid ( figure 4A) show an initial, nearly logarithmic decrease, followed by nearly constant levels of detected photons for increasing depth. These data demonstrate that increased confocality obtained by decreasing the value of δ increases the steepness and length of the initial logarithmic part of the curve. The slope asymptotically approaches the slope of the extinction function exp(−2µ s z / cos α ) as δ is reduced, and approximates this function quite well down to an optical depth of 3 for δ = 0.5 • (squares). No simulation result had a steeper slope than the extinction function, but values of δ smaller than 0.5 • could not be reliably tested due to insufficient numbers of detected photons. At an optical depth of 3 to 4 mean free paths a transition occurs where the number of detected photons levels off into plateau values independent of optical depth. The plateau values were found to be approximately proportional to the square of δ.
The photon counts for blood ( Figure 4B) do not show the same clear distinction between an initial logarithmic decrease and a plateau region. The broken straight line in Figure 4B shows the extinction function exp(−2µ s z / cos α ), using a value of µ s = 150 mm −1 for blood. Values of the confocality angle lower than 0.3 • were not tested due to an insufficient number of detected photons. The initial slope of the curve may approach the extinction curve as δ is further decreased.
To investigate the transition in photon counts from the logarithmic decrease to the plateau region, outputs from the scattering events distribution module were studied for optical depths characterizing the two regions. Figure 5 shows distributions of the number of scattering events per detected photon at the probing depth of 300 µm in 5% intralipid (A) and 50 µm in 4× diluted blood (B), corresponding to optical depths of 3 and 1.9, respectively. Figure 5A shows that as the confocality is increased by decreasing δ while in the logarithmic region for intralipid, the number of detected photons decreases and the scatter distribution is systematically skewed towards fewer scattering events per detected photon. Thus, the mean number of scattering events per photon decreases from 4.2 for δ = 2 • to 1.9 for δ = 0.5 • . The same effect can be observed in blood at shallow optical depth (figure 5B), but here the effect is less pronounced. The mean number of scattering events is only reduced from 4.3 for δ = 2 • to 3.6 for δ = 0.5 • . At an optical depth of 7 in intralipid or 7.5 in blood, the decrease of δ caused a reduction in the number of detected photons without any significant change in the average number of scattering events per photon. The latter parameter only decreased from 15.3 to 14.5 for intralipid and from 15.6 to 14.6 for blood when δ was changed from 2 • to 0.5 • (results not shown).
To illustrate further the difference between the logarithmic and plateau regions observed in figure 4A, figure 6 shows output from the path lines generator module. The graphs in the two upper panels represent the same detection parameter values as used for figure 5.  Figure 6A shows a well-localized focus of backscattering at an optical depth of 3 in intralipid, apart from one stray photon. This is the photon that suffered nine scattering events (see triangles, figure 5A). Note that eight of the remaining sixteen photons in figure 6A experienced two scattering events, whereas the other eight were single scattering events (see triangles, figure 5A). No apparent deviation from the straight single-backscattering paths can be observed in figure 6A for the twice-scattered photons. This suggests that the interaction experienced in addition to the backscattering event for each of these photons must have been a highly forward-directed scattering event. In fact, after examination of the appropriate numerical output, we found that the average of cos ϕ for the eight forwarddirected scattering events was 0.95, where ϕ is the scattering angle at each interaction. This is in contrast to the value of g = cos ϕ = 0.7 used as input parameter for the simulation of scattering in intralipid. Figure 6C demonstrates that the situation is totally different in the plateau region at greater optical depths. Only stray photons are detected when focusing at a nominal optical depth of 7 in intralipid.
In blood, at the shallow optical depth of 1.9 (figure 6B), all photons show well-focused paths. For the optical depth of 7.5 ( figure 6D) the paths are still quite well-focused, apparently due to the highly forward-directed scattering in blood.
The selective detection of minimally scattered photons from shallow optical depths in intralipid, corresponding to the logarithmic region in figure 4A, is more systematically documented in figure 7. The actual number of scattering events per detected photon is expressed as a ratio relative to the expected number of scattering events given by the round trip pathlength of the photon measured in number of mean free path units. For values of optical depth greater than about 5, the ratios are approximately constant around 1.0, corresponding to the expected number of one scattering event per mean free pathlength. In the interval of optical depths from 1 to 4 the values are less than 1, and these minimum ratios decrease with decreasing values of δ. Ratios higher than 1 are found as the optical depth approaches zero. This is due to the fact that all detected photons have necessarily experienced at least one scattering event whereas their expected round-trip path for µ s z < 0.5 will be shorter than one mean free pathlength. The results in figures 5, 6 and 7 demonstrate that for a high degree of confocality (small value of δ) and shallow optical depths, detection is highly selective in that only minimally scattered photons are registered. Compared with average scattering histories, minimally scattered photons have fewer scattering events per optical pathlength travelled, are more forward-scattered or have only short path segments that deviate significantly in direction from a straight, two-segment path corresponding to a single (back-)scattering event at the nominal focus depth of the probe. In this case of minimal scattering, the signal intensity or number of detected backscattered photons can be asymptotically represented by the extinction model. This model is a much better approximation to the simulation data for intralipid than it is for blood. The approximation only applies for shallow optical depths, which for blood translates into distances too short to leave the model much room for applicability (an optical depth of 3 corresponds to a metric depth of 20 µm in blood).

Signal localization
The results in Figure 6 suggest that minimally scattered photons will represent a welllocalized signal, and therefore that good spatial resolution will be obtained in the minimal scattering regime (i.e. optical depth <4 for intralipid).
Localization of signal was systematically investigated by means of the point spread function module which generates one-dimensional distributions showing the x-, y-or z- localization of the maximum depth interaction point for all photons detected for a specified focus position and a set of detection parameters. Figure 8 shows examples of depth point spread functions (zPSF) at shallow optical depths, i.e. at the depth of 300 µm in 5% intralipid (figure 8A, optical depth of 3.0) and 50 µm in 4× diluted blood (figure 8B, optical depth of 1.9). In such diagrams, the numerical aperture of the detector only influenced the yield of detected photons, and the maximum value of NA = 0.4 was therefore used to increase the number of detected photons. Figure 8A illustrates how the width of the point spread function is decreased by decreasing the confocality angle (circles versus squares), thus introducing the selective detection of minimally scattered photons in intralipid. This applies also for the lateral point spread functions yPSF and xPSF (see later results, figure 10). A selective effect on the width of the depth point spread function zPSF is, however, caused by changing the source coherence length. Such a change in general also influences the yield of detected photons. Both effects are illustrated in figure 8A by the different widths and integrals of the two graphs for δ = 0.5 • corresponding to L c = 50 µm (squares, broken curve) and L c = 14 µm (triangles, full curve). In the latter case the detected photons are confined to one 5 µm wide pixel at 295 µm depth. A similar effect is seen in figure 8B for blood, where a shorter source coherence length gives nearly proportionally smaller width of the zPSF (L c = 14 µm, triangles, full curve; versus L c = 50 µm, squares, broken curve; both for δ = 0.5 • ). By comparison with intralipid (figure 8A), there is little stray scattering in blood, and more direction-selective detection obtained by decreasing the confocality parameter therefore does not decrease the width of the zPSF in figure 8B (δ = 2 • , circles; versus δ = 0.5 • , squares; both for L c = 50 µm).
From statistical parameters calculated for the point spread functions it was generally found that in the minimal scattering regime the full width of the zPSF was approximately equal to L c /2, for example a width ≈25 µm for L c = 50 µm and δ = 0.5 • in figure 8A (squares). This corresponds to a pathlength difference of L c between the two round-trip photon paths defining the width of the zPSF, as expected for OCT signal detection. The full width of a Gaussian-like distribution is well represented by four standard deviations (i.e. ±2σ ). In agreement with this, the values of 4σ for the zPSF were found to be comparable to L c /2 in the minimal scattering regime. At a depth of 700 µm in 5% intralipid (optical depth 7.0), no decrease in width of the point spread functions was found with increasing confocality or decreasing source coherence length. Directionality was, however, so well preserved in blood that decreasing L c still reduced the width of zPSF at an optical depth of 7 (results not shown).
Systematic investigations of mean values and standard deviations of the point spread functions for the z-, y-and x-directions are shown in figures 9 and 10.
In figure 9 the mean position values of the PSFs are expressed in terms of mean free path units and as a function of optical depth. These values therefore give the distance (in mfp units) to the mean z, y or x backscatter position for photons detected with the probe aimed at optical focus depth µ s z .Values for two different intralipid concentrations give the same result when plotted in terms of mfp units (figure 9A open versus full symbols). For intralipid the actual depth position probed (circles, figure 9A) closely follows the optical focus depth µ s z within the minimal scattering regime (µ s z < 3).
Deeper into the medium, the depth probed falls significantly short of the nominal focus position. Furthermore, a systematic, positive off-set value for the mean y-position is observed (squares, figure 9A). The latter effect is due to the probe being tilted towards the positive y-axis (α = 15 • ) and will be more pronounced at larger incidence angles. As expected from symmetry considerations, the mean x-position shows no systematic deviation from 0. The more detailed plot for shallow optical depths inserted in figure 9A indicates that position accuracy is marginally improved by decreasing δ from 2 • to 0.5 • .
In blood (figure 9B) a proportional increase in detected z-position with increasing optical depth is observed, and the detected photons are laterally well localized to the z-axis (i.e. xmean = 0, y-mean = 0). This relationship is well preserved to optical depths over 35, apparently due to highly forward-directed scattering in blood.
Standard deviations of the PSFs are shown in figure 10, expressed in dimensionless mfp units. Within the minimal scattering regime in intralipid (µ s z < 3, figure 10A), the standard  deviations are about equal in the z, y and x-directions, and seem to increase with the square of optical depth. At greater depths the x-and y-standard deviations increase approximately linearly with depth, corresponding to relative values (coefficients of variance, CV) of about 20%. The increase in standard deviation of the zPSF is less. This result is reasonable since the pathlength gate limits the range of acceptable z-positions more effectively than the range of x-and y-positions.
A more detailed plot of the values for shallow optical depths in intralipid is shown in the insert in figure 10A, where also data for higher confocality are included. This plot shows that with increasing confocality the standard deviations remain small over a larger range of optical depths, corresponding to an extension of the minimal scattering regime.
Comparison of the inserts in figures 9A and 10A show that precision in intralipid (i.e. low standard deviation) is lost just before the accuracy in mean position values starts to deteriorate as the probe is focused deeper than the extension of the minimal scattering regime.
Results for blood are shown in figure 10B. As for intralipid, the zPSF has a smaller standard deviation than the x-and yPSF due to pathlength gating. For greater optical depths the relative standard deviation (i.e. CV) of x-and yPSF is about 10%, which is about half of the value found for intralipid. The data in figures 9B and 10B are shown for a source coherence length of L c = 14 µm and confocality angle of δ = 2 • . It was found that a low value of L c resulted in a sharper zPSF even at great optical depths in blood, whereas no significant differences were observed between confocality parameters of 2 and 0.5 • . A value of δ = 2 • was therefore chosen to obtain an increased number of detected photons. Figure 10 shows that depth resolution in blood remains remarkably good down to large optical depths, for example with a relative standard deviation of 4% at an optical depth of 35. In contrast, depth resolution in intralipid deteriorates sharply, with relative standard deviations exceeding 10% already at an optical depth of 3. To further investigate the eventual loss of depth resolution in blood with increasing depth, simulations were performed for depths of 400 and 500 µm in undiluted blood. Relative zPSF standard deviations of 10 and 12% were found for the two depths, which correspond to optical depths of 60 and 75 respectively (results not shown). If a relative standard deviation of zPSF of 10% is selected to define the critical depth where localization of signal is lost, this depth is about 3 for intralipid and 60 for blood. It should, however, be noted that these values correspond to reduced optical depths (µ s z = (1 − g)µ s z ) of 0.9 and 0.6 for intralipid and blood respectively, which are much more comparable values.

Discussion
The present investigation assumes a homogeneous scattering medium characterized only by its absorption and scattering coefficients in addition to the anisotropy parameter which expresses the degree of forward-directed scattering. The stochastic Monte Carlo model for photon propagation within the medium was combined with a geometrical optics model of the OCT probe, and interferometric detection, which was implemented according to the antenna theorem of a heterodyne receiver (Schmitt et al 1994b, Siegman 1966. The interferometric gating on pathlength and the high degree of confocality inherent in the heterodyne receiver constitute a highly selective detector in terms of accepting only photons backscattered from the focus region of the probe. This high selectivity translates into very low photon yields in the Monte Carlo simulation. To improve photon yield, the numerical aperture of the simulated detector was usually set to 0.4 instead of 0.2 which is the numerical aperture of the experimental instrument probe. This increased the yield of detected photons by about a factor of 4 without causing any specific effect on other characteristics of the simulation result. Nevertheless, for simulations at considerable depth, the number of detected photons was low, and we could not implement the full confocality of our physical OCT detector (δ = 0.06 • ) due to the excessive computer run times that would be required. However, our results clearly show how the data are influenced by the various parameters, and in some cases also suggest the asymptotic behaviour for limiting values of the instrument parameters. The stochastic variability due to low photon counts is readily observable in most graphs presented in this paper, and the validity of the data can in most cases be assessed by considering the average over neighbouring data points and the apparent precision represented by the stochastic deviation between neighbouring data points.
The advantage of the present exploratory model approach was that new observation functions could be implemented as our simulation revealed new aspects of the OCT phenomenon. For further studies, more specialized Monte Carlo models could be developed that could be optimally designed for answering specific questions with high statistical precision. For such tasks, only photon histories directly relevant for the specified output function would be generated. Methods such as interaction forcing and next-to-last collision analysis (Aronson et al 1991) could be utilized for improving the statistics of results based on rare events, such as detecting photons in the direction of a point-collimated receiver.
We chose intralipid and blood as model media to study since the present work is part of a larger study to investigate flow of blood in a surrounding turbid medium, where intralipid will be used to represent the latter (Lindmo et al 1998). Intralipid and blood represent very different scattering properties (g = 0.7 versus g = 0.99).
Our simulation results show that for a moderately anisotropic medium such as intralipid, the highly selective confocal-like detection inherent to the OCT probe causes only minimally scattered photons to be detected from depths down to several scattering mean free paths (figure 4A). In the minimal scattering regime, the intensity depth distribution approximates a simple extinction model according to exp(−2µ s z / cos α ). Minimally scattered photons have fewer than average scattering events, have more forward-directed scattering than average, and have only short path segments that may deviate significantly in direction from nominal single (back-)scattering paths (figures 5A, 6A, 7). Minimally scattered photons represent backscattering which is well localized to the nominal focus position of the probe, resulting in narrow point spread functions with accurate mean position values (figures 8A, 9A, 10A).
However, beyond the minimal scattering regime, photon counts no longer follow the extinction model, but rather present a constant level of detected stray photons which shows no further decrease with increasing optical depth in a non-absorbing medium (figure 4A). In the stray scattering regime the backscatter positions of detected photons no longer relate to the nominal focus position of the probe ( figure 6C). Even if the detected photons by virtue of the detection procedure must have path lengths corresponding to the nominal focus depth, they are systematically backscattered from shallower depths due to erratic path trajectories (figure 9A). The transition between the minimal scattering and stray scattering regimes is characterized by a steep increase in standard deviations of the point spread functions (at an optical depth of 2 in figure 10A) leading to a subsequent gradual loss of accuracy in position mean values (at optical depth 3 in figure 9A).
The highly forward-directed scattering in blood combined with the difficulty in studying very high confocality because of low photon yield, make it difficult to identify a regime of minimal scattering in blood. However, the results in figures 4B and 5B suggest that a minimal scattering regime with the same qualitative characteristics also exists for blood at shallow depth (<3 mfp). This depth, however, corresponds to such short metric distances (20 µm, round-trip 40 µm) that the phenomenon will be of little practical importance.
The highly forward-directed scattering in blood leads to preservation of localized backscattering from the nominal focus position even for greater optical depths ( figure 6D). Since directionality is preserved in spite of multiple scattering, the depth resolution is improved by a short source coherence length, and accurate position localization may be obtained for greater optical depths ( figure 9B). Depth resolution, as expressed by the standard deviation of the depth point spread function, is maintained at the value determined by the source coherence length down to optical depths of about 10. For greater depths, the relative standard deviation (CV) of the zPSF is about 2.5%, and eventually increases to 10% and more beyond an optical depth of 60 (results not shown). In intralipid 10% relative depth standard deviation is already found at an optical depth of 3. However, the reduced optical depths corresponding to 10% depth resolution are 0.9 and 0.6 units of reduced scattering mean free pathlengths for intralipid and blood respectively. The latter values are quite comparable, and suggest that loss of depth resolution occurs at a threshold value of reduced optical depth rather than optical depth, and that depth resolution is essentially lost in signals detected from depths corresponding to 1 reduced scattering mean free pathlength, i.e. 1/(1 − g)µ s , corresponding to 1.7 mm in 1% intralipid and 0.67 mm in blood.
Our experimental data ( figure 3) give support to the relevance of the simulation results. For intralipid (figure 3A), the agreement with the extinction model is good for the initial part of the curve. For greater optical depth, the experimental curve deviates significantly from the extinction model but does not show a plateau as demonstrated by the simulation data ( figure 4A). Other effects, such as deviation from perfect focus tracking for increasingly greater scanning depths, or non-negligible absorption, may contribute to the slope of the final segment of the experimental intensity depth profile in intralipid. For blood the initial decrease in signal is very steep. The data for 4× diluted blood seem to agree better with the extinction model than the data for undiluted blood.
The extinction model which sees the OCT signal process as a continuous extinction of the photon beam going into the medium, followed by a single (back-)scattering event and continuous extinction of the beam on the way out, has been taken as a starting point for several investigations of OCT (Schmitt et al 1993a, b, 1994a, Yadlowsky et al 1995a. Some studies relate to confocal microscopy, but are still relevant for OCT since by the antenna theorem for a heterodyne receiver (Siegman 1966), OCT is functionally a confocal microscope with a narrow (50 fs) time gate or equivalent gate on photon pathlength (Schmitt et al 1994b, Yadlowsky et al 1995a. In agreement with our results in figure 4A, others have found that the initial slope of signal intensity depth profiles is given by the extinction model, and that the single scatter contribution is dominant down to optical depths of about 3 to 5, whereas for larger depths a less steep slope is observed due to multiple scattering (Pan et al 1997, Schmitt et al 1994a, Yadlowsky et al 1995a. In view of our results (figures 4A, 5A) it is clear that in addition to the perhaps dominant single scattering events, the minimal scattering regime includes photons that have experienced more than one scattering event. The additional scattering events are very much forward-directed and/or characterized by very short path segments deviating significantly from the nominal single backscattering path. The intensity will get asymptotically closer to the extinction function as confocality is increased, and at the same time the minimal scattering regime will be extended to greater depth.
An extension of the extinction-single-backscatter model has been formulated, which includes an additional term due to multiple scattering (Yadlowsky et al 1995b). The contribution from multiple scattering was modelled from statistical optics and shows a depth dependence of the form exp[−µ s (1−g)z], i.e. governed by the reduced scattering coefficient µ s = µ s (1 − g). In our simulations for intralipid using g = 0.7, the reduced scattering coefficient would be µ s = 0.3µ s . The constant level of photon counts demonstrated by the curve in figure 4A for greater optical depths is, however, not in agreement with a decrease corresponding to 30% of the slope of the extinction curve shown in the same figure.
Several investigators have studied axial localization of signal in confocal microscopy and OCT by Monte Carlo simulation. They found that the coherence or time gate inherent to OCT leads to much better depth resolution than confocality alone (Izatt et al 1994, Schmitt et al 1994b. On the other hand, time gating alone gives very poor depth resolution since it permits detection of photons multiply scattered close to the surface (Schmitt et al 1994b). Such effects are demonstrated also by our data in figure 4A for low confocality (δ = 2 • ). Yadlowsky et al (1995b) showed that blurring of the lateral point spread function, which would imply an increase in the standard deviation of the PSF, was more pronounced for more anisotropic scattering (g = 0.96 versus g = 0.90) at the same optical depth (8 mfp units). This result seems to be in disagreement with our findings (for example figures 10A and 10B) that demonstrate better localization of signal both laterally and axially for more forward-directed photon scattering (i.e. higher g-value).
An expression for the standard deviation in axial position of a multiply forward-scattered photon due to Ishimaru (1978) has been used in an analysis by Schmitt and Knüttel (1997). If this expression is applied to our simulation of scattering in blood (g = 0.99, θ rms = 8.1 • = 0.141 rad), we obtain for the axial standard deviation in mfp units µ s S z = µ s z cos α θ rms 2 2 = 5 × 10 −3 (µ s z ) 2 .
If plotted in figure 10B, this expression would yield values quite similar to the standard deviations of xPSF and yPSF, which are considerably larger than the standard deviations of zPSF (circles). Evidently, the coherence gate and the directional selectivity of the OCT detector towards detecting photons that have undergone less than average scattering yield axial depth resolution better than the average uncertainty in axial position of all photons that have propagated the same pathlength as nominally required for the focus position of the probe. The present simulation study only considers homogeneous scattering properties of the medium. It does not consider inhomogeneities in refractive index or other medium properties and no spatial correlations in the medium. Such properties will influence OCT detection in real media. One such effect was pointed out by Yadlowsky et al (1995a) who concluded that the multiple forward scattering from anisotropic scatterers produces a net probe field that is stronger than predicted by the beam extinction theory. The excess probe field has less effective spatial coherence than the unscattered field and therefore leads to higher effective backscattering for structures smaller than the transverse correlation length of the probe field. The multiple scattered light can therefore lead to a depth-dependent shift in the relative backscatter contributions from large and small structures.