## 1.

## Introduction

Parallelized confocal fluorescence systems are becoming more prevalent in the laboratory and in the clinic because they operate at very high speeds, enabling video rate imaging and real-time visualization of biological processes.^{1, 2} Confocal microscopes use a small aperture to reject-out-of focus light, allowing imaging of thin sections within thick samples. Standard confocal microscopes employ a single pinhole aperture that must be spatially scanned to collect a 2-D or 3-D image. To reduce image acquisition times, parallelized confocal systems use an array of pinholes or a slit aperture to simultaneously collect multiple image points, reducing acquisition times in proportion to the number of simultaneous detection points used. The drawback of parallelized systems is cross talk between the individual apertures. In highly scattering media, such as tissue, the cross talk can be large, resulting in a significant reduction in system performance.

Although the theoretical optical sectioning properties of fluorescence confocal systems for various aperture configurations have been studied,^{3, 4} the actual performance in tissue is typically degraded due to scattering. Previous studies have characterized confocal systems with pinhole apertures in turbid media.^{5, 6} However, little work has been done to characterize slit apertures and pinhole array apertures in turbid media. In this paper, we characterize the performance of parallelized confocal fluorescence systems imaging in turbid media. For comparison, we also characterize single pinhole aperture performance under the same conditions.

Existing real-time parallelized confocal systems use a variety of aperture configurations that can be broken down into three general types. Figure 1a illustrates a single pinhole aperture and Figs. 1b, 1c, 1d illustrate the general categories of parallelized apertures: slit, Nipkow, and a linear array of pinholes. To collect a 2-D image and cover an equivalent field of view (illustrated by the box in Fig. 1), each of the four apertures must employ a scanning system. A pinhole aperture requires scanning in two axes. A slit composed of many adjacent detection points requires only one scan axis. Nipkow apertures are large disks with many pinholes that are rapidly rotated to cover the imaging field. Last, a linear array of pinholes (hereafter referred to as a linear array) requires scanning in two axes like a single pinhole, but the range of motion in one axis can be reduced. Each of the parallelized apertures has benefits and limitations in terms of its light efficiency, optical sectioning, scattering cross talk, scanning time, and instrumentation complexity.

A Monte Carlo model was implemented that simulates the confocal system shown in Fig. 2 . The system consists of a laser source and optical elements that uniformly illuminate the confocal aperture. The aperture can be either a pinhole, slit, Nipkow disk, or linear array. The aperture and illumination beam are imaged into the tissue via the objective lens. Fluorescence signal is collected by the objective lens and imaged back onto the confocal aperture. Light passing through the aperture is brought back to focus. A dichroic beamsplitter directs the emitted fluorescence signal to a detector. In the case of a pinhole aperture, a single detector is used; in the case of a parallelized aperture, an array of detectors is used.

In a nonscattering medium, a confocal imaging system has the ability to reject a significant amount of light generated away from the plane of focus. Signal generated at the imaging plane comes to focus at the confocal aperture and passes through to the detector. Light coming from a point out of focus produces a defocused beam at the aperture because it is not conjugate to this plane. Thus, a significant proportion of the beam’s energy will be rejected by the aperture, resulting in reduced signal from out-of-focus planes. This is what gives the confocal microscope its *optical cross-sectioning* property. Ideally, for a point lying in the plane of focus, all light generated over the numerical aperture (NA) of the objective should pass through the aperture to maximize the collected signal. However, due to diffraction, an infinitely large aperture would be required to collect all the light. Moreover, to reject all light from out-of-focus planes, an infinitely small pinhole would be required. In practice, a pinhole size somewhat larger than the width of the Airy diffraction pattern will provide a good signal-to-noise ratio and reasonable rejection of background signal.^{3}

In a scattering medium, some of the illumination light is scattered out of the beam path, increasing the illuminated volume. When induced fluorescence from scattered illumination couples back into the collection path, the axial and lateral resolution degrade.

Tissue scattering, nonradiative absorption, and fluorophore concentration all effect the system performance. The excitation beam can scatter before reaching the imaging plane, resulting in nonuniform illumination with tissue closer to the surface receiving more excitation energy. This nonuniform illumination leads to an increased probability of fluorescence signal generation near the tissue surface. Tissue nonradiative absorption and fluorophore concentration amplify the surface bias effect. These effects reduce the energy available for conversion to fluorescence signal at the focal plane. Another event that reduces system performance is scattering of the fluorescence emission. Fluorescence signal generated at out-of-focus planes can scatter into the collection path, reducing the axial and lateral performance.

In the following sections, we describe the Monte Carlo model that was implemented to simulate parallelized fluorescence confocal systems imaging in turbid media and then present a comparative performance analysis of the four confocal aperture configurations.

## 2.

## Monte Carlo Model of a Fluorescence Confocal System

A direct Monte Carlo^{7} simulation of individual photons traveling through turbid media was implemented. The model is based on previous work by Wang,^{8} Wilson and Adam,^{9} and Prahl
^{10} Although variance reduction schemes^{11, 12} and multiple independent run methods^{13, 14} exist to increase simulation speed, the direct method was chosen because it makes fewer assumptions and is less susceptible to binning artifacts. With modern computers and careful attention to efficient implementation, the direct model was sufficiently fast to generate the required number of photons for accurate modeling.

Conceptually, the direct Monte Carlo model of a fluorescence confocal system imaging in turbid media starts with the creation of an excitation photon in the confocal aperture. The photon’s position is randomly generated using the intensity distribution in the confocal aperture. Its direction is chosen to be within the optical system’s numerical aperture. Then the photon is propagated through the system optics until it reaches the surface of the tissue, where it refracts and enters the tissue. In tissue, the photon is repeatedly propagated and scattered until it (1) is absorbed and emitted as excited fluorescence energy, (2) is nonradiatively absorbed, or (3) exits the surface of the tissue. When a photon is reemitted as excited fluorescence, the scattering process continues until the photon is either nonradiatively absorbed or exits the tissue surface. In the case where the photon is absorbed nonradiatively, the photon is terminated. If the photon exits the surface, it refracts and is then traced back through the optical system. If an excited photon makes it into the NA of the collection system and passes through the confocal aperture, then it is recorded as collected energy. All other photons that exit the tissue surface are rejected. The process is repeated until a sufficient number of photons are collected to enable analysis of the axial and lateral system response.

## 2.1.

### Model Principles

The implemented Monte Carlo model makes six assumptions. First, fluorescence emission is assumed to be isotropic.^{6, 15} Second, we assume that scattering dominates over diffraction effects so that diffraction can be ignored. Previous studies have shown that this is valid for a confocal system in turbid media.^{16} Third, since the diffraction-limited properties of a confocal system are known^{4} and we are interested in studying only the effects introduced by scattering, we model the optical system as ideal with uniform aperture illumination and uniform energy across the NA. Fourth, the detection system is ideal; all collected photons are detected. Fifth, we assume that the tissue absorption and scattering properties do not dramatically change between the excitation and emission wavelengths. Sixth, we assume that the fluorophore quantum conversion efficiency is one and that the fluorophore is uniformly distributed on a macroscopic scale in the media.

Because our model assumes an ideal optical system, propagation of the initial photon from the confocal aperture through the optics to the tissue can be done efficiently. Since the confocal aperture is conjugate to the tissue focal plane, we can avoid skew ray tracing and use the magnification between these two planes to move the photon to the imaging plane. Of course, in turbid media, there is no direct mapping between these planes, but we can envision that the photon follows a path that satisfies this mapping until it reaches the surface of the tissue. Thus, an efficient approach is to “generate” the excitation photon in the ideal image of the confocal aperture in the tissue focal plane and geometrically back-project the photon to the surface of the tissue. At the surface, the photon’s trajectory and position are the same as they would have been if the photon had been propagated via skew ray tracing through the optical system. Using this same logic, when an fluorescence photon exits the tissue, it is geometrically back-projected to the tissue focal plane, where it can be analyzed to determine whether it falls within the image of the confocal aperture and the collection NA of the optical system. In other words, the position and angle of this back-projected photon will determine whether the photon exciting the tissue will be collected by the optical system.

By using this method of generating and collecting the photons, we allow our analysis to be completely carried out in the tissue space. Presenting results in the tissue space is more intuitive and directly useful for characterizing how the system performs. Moreover, as long as the confocal aperture and lens NA can be described in tissue space, the Monte Carlo model can be described independently of the optical system and its magnification.

Figure 3 shows an example of a photon that is absorbed and reemitted as a fluorescence photon and then collected by the system. First the excitation photon’s spatial position and trajectory are generated randomly in the aperture image at ${\mathit{r}}_{\mathbf{0}}$ . Then the photon is back-projected to where it encounters the tissue surface at ${\mathit{r}}_{\mathbf{0}}^{\prime}$ . At this location, it is launched back into the turbid medium and scattered until it is absorbed by the fluorophore and fluoresces at ${\mathit{r}}_{\mathit{s}}$ . The fluorescent photon’s new direction is randomly selected from an isotropic distribution and then scattered until it exits the tissue at ${\mathit{r}}_{\mathit{d}}^{\prime}$ . To determine if the photon is collected by the system, it is back-projected to ${\mathit{r}}_{\mathit{d}}$ . Since ${\mathit{r}}_{\mathit{d}}$ falls within the aperture’s opening and the photon’s angle is within the maximum collection angle ${\varphi}_{\mathrm{max}}$ as determined by the NA, the photon is recorded as collected signal. Ideally, in a nonscattering media with an infinitely small pinhole, ${\mathit{r}}_{\mathit{s}}$ will equal ${\mathit{r}}_{\mathit{d}}$ .

To account for the effects of detector sampling, the photon’s final position is recorded as the center of the $i$ ’th pinhole or detector element that it falls within ${\mathit{r}}_{\mathit{c}\mathit{i}}$ . Figures 4a and 4b depict how ${\mathit{r}}_{\mathit{c}\mathit{i}}$ is determined for pinhole-based apertures. We assume that each pinhole has one detector behind it, resulting in ${\mathit{r}}_{\mathit{d}}$ being recorded as the center of the $i$ ’th pinhole ${\mathit{r}}_{\mathit{c}\mathit{i}}$ that the photon is collected in. Since a slit aperture has many detectors [Fig. 4c], ${\mathit{r}}_{\mathit{c}\mathit{i}}$ is the center of the detector element that contains ${\mathit{r}}_{\mathit{d}}$ . The example shown in Fig. 3, assumes that the aperture is a single pinhole. Therefore, all detected photons are recorded with the collected position of ${\mathit{r}}_{\mathit{c}1}$ .

There is an error in the collected signal as defined by the vector $\mathbf{\epsilon}={\mathit{r}}_{\mathit{c}\mathit{i}}-{\mathit{r}}_{\mathit{s}}$ . The distribution of ${\epsilon}_{z}$ characterizes the axial sensitivity of the confocal system, and the bivariate distribution of ${\epsilon}_{x}$ and ${\epsilon}_{y}$ characterizes the lateral blur induced by scattering. Once ${\mathit{r}}_{\mathit{d}}$ is shifted to the central point in the detection element receiving the photon, the bivariate distribution of ${\epsilon}_{x}$ and ${\epsilon}_{y}$ also accounts for the detector sampling.

## 2.2.

### Implemented Model

Figure 5 depicts the 10 steps in the implemented Monte Carlo model. The model was implemented in C using the GNU Scientific Library.^{17} To generate simulation-quality random numbers, a combined multiple recursive generator^{18} with a period of
${2}^{185}$
(about
${10}^{56}$
) was used. To generate results in a timely manner, the code was run on a
$40\text{-}\mathrm{GHz}$
distributed Xgrid cluster. In the following text, each step is discussed in detail:

**Step 1**. By mapping the uniformly illuminated confocal aperture to the tissue focal plane
$(z=0)$
, the initial spatial position of the photon
${\mathit{r}}_{\mathbf{0}}=({x}_{0},{y}_{0},{z}_{0}=0)$
can be randomly generated.

In general, a random variable $X$ with an arbitrary probability density function $f\left(x\right)$ can be generated using a random variable $\xi $ uniformly distributed between 0 and 1. (Ref. 19). To generate $X$ , first the cumulative probability density function $F\left(x\right)$ is computed using

Then $F\left(x\right)$ is inverted to yield the inverse cumulative probability density function ${F}^{-1}\left(x\right)$ . Last, the random variable $X$ can be generated asNote that in the following sections, each instance of $\xi $ represents a new uncorrelated sampling of a uniformly distributed (0 to 1) random variable.Since we assume that the aperture is uniformly illuminated, the initial ${\mathit{r}}_{\mathbf{0}}=({x}_{0},{y}_{0},{z}_{0}=0)$ coordinates for the photon in the aperture image can be generated for a single pinhole with radius $a$ using

## 3.

Once the initial photon position is generated, its initial trajectory must also be randomly generated. Let the optical axis be along
$z$
. Assuming that the intensity of the beam does not vary with angle, the trajectory of the photon can be in any orientation in the
$x\text{-}y$
azimuthal plane and maximally depart from the
$z$
axis by the zenith angle
${\varphi}_{\mathrm{max}}={\mathrm{sin}}^{-1}(\mathrm{NA}\u2215n)$
, where NA and
$n$
are the numerical aperture and index of refraction in the turbid media. The initial trajectory can be described with the random variables
$\theta $
and
$\varphi $
, which are the zenith angle, and azimuthal angle, respectively. The azimuthal angle
$\theta $
has a uniform distribution that ranges from 0 to
$2\pi $
. The zenith angle requires a uniform sampling of
$\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\varphi $
for *phi* ranging from 0 to
${\varphi}_{\mathrm{max}}$
in order to achieve uniform sampling within the NA.

**Step 2**. This step consists of back-projecting the photon from
${\mathit{r}}_{\mathbf{0}}$
to
${\mathit{r}}_{\mathbf{0}}=({x}^{\prime},{y}^{\prime},{z}^{\prime}=-d)$
at the surface of the turbid media so that the scattering process can begin.

In general, to propagate a photon, the new coordinate can be related to the old coordinate by

## 4.

## 5.

**Step 3**. With the photon position and angle initialized, the propagation and scattering process begins. The photon can be scattered through the semi-infinite volume of the tissue until it is absorbed or exits the surface. The relevant tissue properties are
${\mu}_{s}$
$\left[{\mathrm{cm}}^{-1}\right]$
,
${\mu}_{a}$
$\left[{\mathrm{cm}}^{-1}\right]$
, and
$g$
[unit-less], which are the scattering coefficient, the nonradiative absorption coefficient, and the anisotropy factor describing the scattering angle distribution. The fluorophore properties are described by the fluorescence coefficient
${\mu}_{f}=\Gamma \cdot C$
$\left[{\mathrm{cm}}^{-1}\right]$
, which can be expressed in terms of the fluorescence extinction coefficient
$\Gamma $
$[{\mathrm{cm}}^{-1}\cdot \mathrm{L}\u2215\mathrm{mol}]$
(assuming a quantum efficiency of one) and the fluorophore concentration
$C$
$[\mathrm{mol}\u2215\mathrm{L}]$
.

The mean free path a photon travels before it scatters is $1\u2215{\mu}_{s}$ . The distribution of path lengths $\Delta s$ that photons travel before scattering is described by the Beer-Lambert law probability density function;

## Eq. 7

$f\left(\Delta s\right)={\mu}_{s}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\Delta s\cdot {\mu}_{s}).$^{20}

## Eq. 9

$f\left(\Delta \varphi \right)=\frac{1-{g}^{2}}{{(1+{g}^{2}-2g\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\Delta \varphi )}^{3\u22152}}.$## Eq. 10

$\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\Delta \varphi =\frac{1}{2g}[1+{g}^{2}-{\left(\frac{1-{g}^{2}}{1-g+2g\xi}\right)}^{2}].$^{21}have shown experimentally that the Henyey-Greenstein phase function provides a good characterization of single scattering events in tissue. In the limit of $g=0$ in Eq. 10, scattering is isotropic, and $\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\Delta \varphi $ is uniformly distributed between $-1$ and 1.

Given the scattered angle $(\Delta \theta ,\Delta \varphi )$ and the propagation direction $(\alpha ,\beta ,\gamma )$ of the incoming photon, the new direction $({\alpha}^{\prime},{\beta}^{\prime},{\gamma}^{\prime})$ is given by

## 11.

The propagation and scattering process continues until the photon either exits the surface, is absorbed and reemitted as fluorescence, or is nonradiatively absorbed. When an incident photon reaches the tissue surface, it reflects back into the tissue with a probability $R$ , where $R$ is the reflection coefficient. The test condition for reflection is $\xi \u2a7dR$ . The reflection coefficient $R$ is the average of the Fresnel reflection coefficients for the $s$ and $p$ polarization states:

## Eq. 12

$R=\frac{1}{2}{\left\{\frac{{n}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}-{n}_{2}{[1-{\left(\frac{{n}_{1}}{{n}_{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}\right)}^{2}]}^{1\u22152}}{{n}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}+{n}_{2}\sqrt{1-{\left(\frac{{n}_{1}}{{n}_{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}\right)}^{2}}}\right\}}^{2}+\frac{1}{2}{\left\{\frac{{n}_{1}{[1-{\left(\frac{{n}_{1}}{{n}_{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}\right)}^{2}]}^{1\u22152}-{n}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}}{{n}_{1}{[1-{\left(\frac{{n}_{1}}{{n}_{2}}\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}\right)}^{2}]}^{1\u22152}+{n}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}{\varphi}_{i}}\right\}}^{2},$To determine where the excitation photon terminates in the tissue, the path length to fluorescence, ${l}_{f}$ , and the path length to absorption, ${l}_{a}$ , are determined in the same manner that $\Delta s$ was determined using the random generator

andNote that ${l}_{f}$ and ${l}_{a}$ are determined only once per photon. If ${l}_{a}<{l}_{f}$ , the photon will be nonradiatively absorbed before it has a chance to be converted to fluorescent energy. Since we are interested only in the distribution of collected fluorescent signal, the photon need not be propagated. If ${l}_{f}<{l}_{a}$ , the photon will be absorbed and reemitted as fluorescence when its total path length $L$ is equal to the fluorescence path length ${l}_{f}$ . The total path length $L$ is the accumulated path length for the number of scattering eventsPhotons with ${l}_{f}<{l}_{a}$ travel until $L={l}_{f}$ , where they are converted to excited fluorescence photons in the next step. If the next propagation step size $\Delta s$ causes the total path length to exceed ${l}_{f}$ , the propagation terminates partway through that path at the point where the total path length equals ${l}_{f}$ .

**Step 4**: When the photon has traveled a distance
$L={l}_{f}$
, it is absorbed and emitted as fluorescent energy. The position where fluorescence occurs represents the location of the signal and is recorded as
${\mathit{r}}_{\mathit{s}}$
. Since we assume that fluorescent emission photons are isotropically generated from the fluorophore, a new photon direction
$(\Delta \theta ,\Delta \varphi )$
is randomly generated with a uniform probability. The photon path length
$L$
is reset to zero.

**Step 5**: Just as in step 3, the fluorescence photon is propagated and scattered through the tissue with random path lengths
$\Delta s$
. The propagation continues until the photon exits the tissue surface or the photon is absorbed when
$L={l}_{a}$
using a newly generated
${l}_{a}$
.

**Step 6**: If the fluorescence photon reaches the surface and exits, the scattering process is terminated.

**Step 7**: At this step, the fluorescence photon is geometrically back-projected from
${\mathit{r}}_{\mathit{d}}^{\prime}$
on the surface to the point
${\mathit{r}}_{\mathit{d}}$
in the focal plane.

**Step 8**: The angle of the photon at the tissue focal plane is compared to the acceptance NA of the optical system to determine whether the photon is collected. If the angle of the photon from the
$z$
axis is too large, it is rejected.

**Step 9**: The spatial position of the photon is checked to see if it falls within the image of the confocal aperture. If it is outside the aperture, the photon is rejected.

**Step 10**: If the photon is not lost via confocal rejection, NA rejection, absorption, or nonfluorescence, it is detected. To account for the effect of multiple detectors, the final signal position is recorded as the central coordinate of the detector element
${\mathit{r}}_{\mathit{c}\mathit{i}}$
that collects the photon (illustrated in Fig. 4).

Since the modeled tissue is semi-infinite, if the tissue absorption is low enough, it is possible for a photon to scatter many times. A photon that scatters a large number of times has an extremely low probability of being collected because it tends to drift far away from the imaging region. To prevent the simulation from following a photon that has an extremely low probability of being collected in the system, we terminate any excitation photon or fluorescent photon that scatters more than 100 times. In systems where ${\mu}_{a}\u2aa1{\mu}_{s}$ , it is not unusual for an individual photon to undergo a large number of scattering events. However, under the conditions studied in this paper, it was extremely rare for a photon with more than 100 scatter events to be collected. The performance metrics we use to characterize the axial and lateral response are resistant to rare events, and therefore terminating these photons does not affect the results.

With the termination of a photon by collection or it being lost, the process is repeated until the desired number of photons ${N}_{c}$ are collected by the system.

## 3.

## Results and Discussion

The scattering model in the Monte Carlo code was checked by measuring the total reflectance and total transmission for finite slabs of various optical thicknesses and ratios of
${\mu}_{a}\u2215{\mu}_{s}$
with matched boundary conditions
$(n={n}^{\prime})$
. The predicted values were validated against solutions to the radiative transport equation and results of van de Hulst.^{22, 23} Total diffuse reflectance and total diffuse transmission were validated against Giovanelli’s 1955 results^{24} for mismatched boundary conditions (tissue
$n>1$
) for a semi-infinite slab bounded by glass slides. Finally, to verify that the fluorescence confocal aspects were correctly implemented, we were able to reproduce the normalized axial intensity functions modeled by Blanca and Saloma^{6} for a fluorescence confocal system with a pinhole aperture. Since Blanca and Saloma did not report
${\mu}_{s}$
, it was not possible to reproduce the same
$z$
scale, but the relative intensity shapes and heights matched.

Our primary interest was to compare the axial and lateral performance of parallelized apertures to a single pinhole aperture. In all simulations, the confocal apertures were uniformly illuminated. We used realistic parameters based on an existing clinical confocal microendoscope for imaging esophagus tissue.^{1} An NA of 0.5 in tissue space was used. The simulations were run until the number of collected photons
${N}_{c}$
was 100,000 when axial distributions were studied and 10,000 when measures of the lateral and axial distribution’s spread were studied. Our results show that the values of
${N}_{c}$
were large enough to limit the error of the spread measurements to two orders of magnitude less than the nominal value. To achieve these values of
${N}_{c}$
for all configurations presented in this paper, some simulations had to generate in excess of 1.5 billion photons. The analysis and visualization of the data were accomplished with R. (Ref. 25).

The four apertures were specified in terms of their size in tissue space for a system designed to image a $450\text{-}\mu \mathrm{m}$ -square region. The pinhole aperture was $1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ in diameter. The slit aperture was $1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ by $450\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , with the long dimension oriented along $y$ . Inside the slit, three hundred $1.5\text{-}\mu \mathrm{m}$ -square detectors recorded signal in parallel. The linear array consisted of a line of pinhole apertures $1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ in diameter spaced $30\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ center-to-center, spanning $450\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ (15 pinholes total). For the Nipkow aperture, a 2-D array of sixty $1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ apertures were spaced $60\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ apart. The linear array and the Nipkow apertures were optimized to maximize the number of apertures while maintaining performance comparable to a single pinhole aperture down to an imaging depth of $62.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . The optimization of the linear array and Nipkow apertures is discussed in Sec. 3.2.

For a pinhole aperture system, an approximate formula for the optimum pinhole radius as determined by the half-power width of the Airy diffraction pattern in tissue space is^{3}

We simulated human esophagus tissue for a range of imaging depths below the tissue surface. Because tissue and fluorophore properties can vary significantly, we analyzed a range of typical tissue and fluorophore parameters ${\mu}_{a}\u2215{\mu}_{s}$ , $g$ , and ${\mu}_{f}\u2215{u}_{s}$ to show how they affect the results. We did not study the effects of the tissue index $n$ since it does not substantially vary, and it affects only the Fresnel reflections at the surface of the tissue.

Unless otherwise specified, we present results for esophagus tissue at a wavelength of $632.8\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ with values ${\mu}_{s}(1-g)=12\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , ${\mu}_{f}=40\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , and ${\mu}_{a}=0.4\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ (Ref. 26). The optical index of refraction modeled was $n=1.37$ (Ref. 27). For a typical tissue value of $g=0.85$ , the scattering coefficient is ${\mu}_{s}=80\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . We also assume that the average fluorophore concentration is uniform throughout the tissue with ${\mu}_{f}=0.5{\mu}_{s}$ . Our analysis used a nominal imaging depth of $d=62.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . However, the results can be generalized for other tissues with different scattering coefficients by expressing the imaging depth in terms of the average number of mean free paths $d\cdot {\mu}_{s}$ . Similarly, we can generalize the axial error from focus ${\epsilon}_{z}$ as the average number of mean free paths from focus ${\epsilon}_{z}\cdot {\mu}_{s}$ .

Our approach in modeling the fluorophore concentration as uniform throughout the tissue is different than previous approaches. Others have modeled a high-frequency fluorescence pattern in the tissue or modeled a plane of fluorescence that is moved through the plane of focus. Modeling a high-frequency fluorescent pattern makes the results dependent on the pattern that was used. A fluorescent plane moving through focus masks the effects of decaying excitation signal due to the absorption of photons at planes above the plane being analyzed. Real tissue will have spatially varying fluorescence due to the cellular structure, but modeling a uniform concentration of fluorophore throughout the tissue approximates the fluorescence as an average fluorescence signal over a volume. This results in simulations that realistically characterize axial sensitivity and properly account for decaying signal with depth due to fluorescence and nonradiative absorption.

## 3.1.

### Characterizing Aperture Imaging Performance

Videos 1 shows the 3-D distributions^{28} of the collected fluorescence photon signal positions
${\mathit{r}}_{\mathit{s}}$
for the four aperture configurations imaging in simulated tissue at an imaging depth of
$d\cdot {\mu}_{s}=0.5$
mean free paths
$\left(62.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$
. In each of the figures, the tissue surface lies in the
$x\text{-}y$
plane at
$z=-62.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, and the plane of focus is at
$z=0$
.
Videos 1(c) and
Videos 1(d) appear more granular than the single pinhole result in
Videos 1(a) because the same number of collected photons is distributed over a larger area.

Figure 6 shows the projection^{29} of the collected fluorescence photon signal
${\mathit{r}}_{\mathit{s}}$
onto the
$x\text{-}y$
plane. This plot identifies how much signal is collected outside the aperture. Although the signal is extremely large within the aperture, the distribution of energy outside the aperture indicates blurring due to scattering. The slit aperture has the most spread out signal, indicating a substantial lateral blur. Both the linear array and the Nipkow aperture have a distribution of photons centered around each aperture that seems comparable to the single pinhole aperture.

It is difficult to directly compare the performance of the aperture configurations using the distribution of fluorescence photon locations shown in Videos 1 and Fig. 6. A better comparison can be made by plotting the error distribution of photon positions $\mathit{\epsilon}={\mathit{r}}_{\mathit{c}\mathit{i}}-{\mathit{r}}_{\mathit{s}}$ , as shown in Videos 2 . In an ideal confocal system, $\left|\mathit{\epsilon}\right|$ should be zero for all collected fluorescent photons. As the distribution of ${\epsilon}_{z}$ broadens, the axial performance of the system decreases. Similarly, as the distribution in the lateral plane broadens, the lateral resolution performance of the system degrades.

10.1117/1.3194131.2Compared to the pinhole aperture in Video 2(a) the slit aperture in Video 2(b) shows significant error, especially near the surface. Video 2(b) also shows an interesting V-shaped region of sensitivity error moving away from the surface. The parallelized pinhole apertures in Video 2(c) and Video 2(d) appear to have distributions nearly identical to the single pinhole except for some subtle signal at intervals of the pinhole spacing.

To more quantitatively compare the aperture performance shown in Videos 2 the data are presented in Fig. 7 as projections of $\mathit{\epsilon}$ on the axial planes ${\epsilon}_{z}\u2013{\epsilon}_{x}$ and ${\epsilon}_{z}\u2013{\epsilon}_{y}$ and the lateral plane ${\epsilon}_{x}\u2013{\epsilon}_{y}$ for the four aperture configurations. Since the majority of the energy is localized near $\mathit{\epsilon}=0$ , the plots are shown with a log color scale.

The axial plots (left and center columns) in Fig. 7 show that the apertures have a strong sensitivity within the NA (cone region) of the optical system. The sensitivity decreases away from $\mathit{\epsilon}=0$ . The pinhole aperture appears to have very little sensitivity outside the NA, whereas the parallelized apertures are more sensitive outside this region. The asymmetry of the slit and linear array apertures is apparent, as they both have increased sensitivity toward the surface in the ${\epsilon}_{z}\u2013{\epsilon}_{y}$ plane. The interesting V-pattern in the ${\epsilon}_{z}\u2013{\epsilon}_{y}$ plot for the slit is due to the cross talk between the detectors along the slit. The repetitive pattern in the parallelized pinhole apertures is a result of cross talk between pinhole apertures. This pattern is not apparent in the slit aperture because the detectors are adjacent to one another.

The lateral error plots (right column) in Fig. 7 illustrate the lateral blur induced by the tissue scattering. The elongation of the slit lateral plot in Fig 7b is due to the cross talk between detector elements in the aperture. The repetitive clusters of signal away from $\mathit{\epsilon}=0$ in the linear array Fig. 7c and the Nipkow Fig. 7d apertures indicates the potential for ghosting in the image.

To understand how the apertures differ in terms of their ability to reject out-of-focus signal, Fig. 8 plots the axial sensitivity for each of the apertures and breaks it down into the ballistic photons (fluorescence photons that do not scatter before collection) and scattered photons. The horizontal axis is given in terms of the number of mean free paths ${\epsilon}_{z}\cdot {\mu}_{s}$ (unit-less) away from focus. The nominal focus is at ${\epsilon}_{z}\cdot {\mu}_{s}=0$ , and the tissue surface is at ${\epsilon}_{z}\cdot {\mu}_{s}=-0.5$ . Each of the three curves are probability density functions describing the probability of collected photons as a function of mean free paths away from focus. For a probability density function, the area under the curve is always one and the vertical axis units are the reciprocal of the horizontal axis units.

Figure 8 shows that the pinhole-based apertures (pinhole, linear array, and Nipkow) have a strong peak sensitivity at focus that rapidly drops off. This represents the desired axial sensitivity in a confocal system. The pinhole-based aperture signals are composed primarily of ballistic photons near the focus. This is indicated by the overlap of the ballistic photon density plot (dashed line) with the density for all photons (solid line). There are small performance differences in the scattered photon density (dotted lines) for the pinhole-based apertures. The parallelized pinhole apertures have more sensitivity moving away from the surface compared to the single pinhole aperture. The slit aperture Fig. 8b has a significant photon density at the tissue surface due to scattered photons.

To illustrate how the confocal system’s sectioning properties change with imaging depth, Fig. 9a shows the ratio of the number of collected ballistic photons to the total number of collected photons as a function of imaging depth (focal position in tissue) $d\cdot {\mu}_{s}$ for each of the apertures. Figure 9b shows the standard deviation of 10 runs for each point in Fig. 9a to estimate the error in the Monte Carlo model. The model errors are typically at least two orders of magnitude less than the estimated values, indicating that ${N}_{c}$ is sufficiently large to reliably estimate the relative ballistic signal. Ideally, the relative ballistic component should be equal to one. The pinhole aperture maintains more than 90% ballistic signal down to $d\cdot {\mu}_{s}=1$ . The linear array and Nipkow maintain at least 90% ballistic signal down to $d\cdot {\mu}_{s}=0.5$ but drop off thereafter. The drop off after $d\cdot {\mu}_{s}=0.5$ is a result of the optimization for imaging down to $d\cdot {\mu}_{s}=0.5$ . The figure highlights that the slit aperture has the worst performance, limited to 71% at best for $d\cdot {\mu}_{s}=0.25$ .

While the relative ballistic signal helps to highlight the difference in the four aperture configurations, it fails to capture the confocal system’s ability to reject defocused light that is ballistic. In general, the spread of the axial response along ${\epsilon}_{z}$ quantifies the axial sensitivity of the confocal system. However, a single measure of axial spread can be misleading since the axial sensitivity function is not a simple unimodal function. To provide a characterization of axial response, we present the axial density function ${\epsilon}_{z}$ for a range of imaging depths, tissue properties, and fluorophore concentrations.

Figure 10 shows how the axial sensitivity for the four apertures varies as a function of departure from focus ${\epsilon}_{z}\cdot {\mu}_{s}$ for three different imaging depths. The pinhole, linear array, and Nipkow apertures have nearly identical performance. The slit aperture clearly has increased sensitivity at the surface.

Quantifying the lateral and axial response in a single number such as FWHM or RMS is problematic because the distributions have long tails. FWHM does not capture the spread in the tails and RMS gives greater weight to extreme values. To quantify the spread of the lateral and axial response, we report the interquartile range (IQR), which represents the range that bounds an area of 0.5 centered about the median value. The IQR is a stable estimate of spread in the presence of extreme values.

In lens design, the lateral response of an abberated system is often related to the diffraction-limited response by comparing the RMS spot diameter to the diffraction limited Airy disk diameter. We can make a similar comparison using the lateral distribution of $\mathit{\epsilon}$ and the Airy IQR. To find the Airy IQR, the intensity of the Airy pattern (Fraunhofer diffraction pattern)

## Eq. 17

$I\left(\rho \right)={I}_{0}{\left[\frac{2{J}_{1}\left(k\rho \mathrm{NA}\right)}{k\rho \mathrm{NA}}\right]}^{2},$## Eq. 18

$P\left(\rho \right)={P}_{0}[1-{J}_{0}^{2}\left(k\rho \mathrm{NA}\right)-{J}_{1}^{2}\left(k\rho \mathrm{NA}\right)].$The diffraction-limited nonscattering axial FWHM responses for pinhole and slit aperture confocal systems have been previously described.^{3} For the pinhole and slit systems modeled here, the axial FWHM is
$6.15\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
and
$8.72\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, respectively. The central lobe of the axial response can be reasonably approximated with a Gaussian profile. Since the standard deviation
$\sigma $
for a Gaussian function in terms of the FWHM is

## Eq. 19

$\sigma =\frac{\mathrm{FWHM}}{{\left(2\phantom{\rule{0.2em}{0ex}}\mathrm{ln}\phantom{\rule{0.2em}{0ex}}2\right)}^{1\u22152}},$Figure 11a presents lateral and axial IQR values as a function of imaging depths $d\cdot {\mu}_{s}$ . To characterize the lateral distribution of $\mathit{\epsilon}$ , we compute the IQR of the signed radial response ${\epsilon}_{\rho}={({\epsilon}_{x}^{2}+{\epsilon}_{y}^{2})}^{1\u22152}\cdot \text{sign}\left({\epsilon}_{x}\right)\cdot \text{sign}\left({\epsilon}_{y}\right)$ . Since the lateral IQR values are at least twice the Airy IQR, the combination of the finite detector size and scattering effects are substantially more than the diffraction-limited performance, which supports our decision to ignore diffraction effects in our model. Similarly, the axial IQR values for the pinhole and slit are also at least as large as the diffraction-limited values. Figure 11b shows the standard deviation of 10 runs for each point in Fig. 11a to estimate the error in the Monte Carlo model. The model errors are typically at least two orders of magnitude less than the estimated values, indicating that ${N}_{c}$ is sufficiently large to reliably estimate the IQR.

Figure 11 shows how the pinhole aperture has a stable axial performance of about $4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ and lateral performance of about $1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ down to $d\cdot {\mu}_{s}=1$ . The linear array and Nipkow aperture maintain performance comparable to the pinhole aperture down to their optimized depth $d\cdot {\mu}_{s}=0.5$ (optimization discussed in Sec. 3.2). The slit aperture has substantially degraded performance. All apertures fail for $d\cdot {\mu}_{s}>2$ , and the point of failure decreases with increasing parallelization (more pinholes).

Figure 11 shows that after a certain depth, the aperture performance appears to stabilize or even improve; however, this is not true. At very deep imaging depths, there is almost no signal being collected at or below the plane of focus. The resultant axial density falls off from the surface with no peak at focus, causing the IQR to shift toward the surface. At these depths, the image is almost completely composed of defocused signal from near the surface.

## 3.2.

### Optimizing Pinhole Spacing

To maximize the speed performance of systems using parallelized apertures, the highest possible aperture density should be used. To determine the maximum possible density that can be used while still maintaining reasonable confocal performance, the maximum imaging depth $d\cdot {\mu}_{s}$ must be specified. Since the axial and lateral performance degrades as the imaging depth is increased, the pinhole spacing $\delta $ should be optimized to obtain the minimum acceptable performance at the maximum required imaging depth. To stay consistent with our previous discussion, $\delta $ is specified in terms of the aperture spacing imaged into tissue space.

Figures 12a and 12b show the axial sensitivities of a linear array for three values of $\delta $ at $d\cdot {\mu}_{s}=0.5$ and $d\cdot {\mu}_{s}=1$ , respectively. As $\delta $ increases, the performance increases. Figures 12c and 12d show how the axial and lateral IQR improve asymptotically toward the limit of the single pinhole aperture as $\delta $ increases. For the target imaging depth of $d\cdot {\mu}_{s}=0.5$ , the axial and lateral performance do not substantially improve beyond $\delta =30\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . Thus, for the system modeled, an optimized linear array would have 15 pinholes spaced $30\phantom{\rule{0.3em}{0ex}}\text{microns}$ apart in a $450\text{-}\text{micron}$ line.

Figure 13 shows the same analysis for optimizing the pinhole spacing in the Nipkow aperture. Since the Nipkow aperture has pinholes in a 2-D arrangement, a larger $\delta $ is required to achieve the same level of performance since cross-talk can occur in two dimensions. Figure 13c and 13d show that the IQR performance does not improve substantially for values of $\delta $ greater than $60\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ Therefore, for the system modeled, a Nipkow disk with 60 pinholes spaced $60\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ apart in a $450\text{-}\mu \mathrm{m}$ -square area is nearly optimal.

## 3.3.

### Characterizing Effects of Tissue Properties

Changes in ${\mu}_{s}$ , ${\mu}_{a}$ , and $g$ can have substantial effects on the axial and lateral performance. Since we present our results in terms of ${\epsilon}_{z}\cdot {\mu}_{s}$ and $d\cdot {\mu}_{s}$ , they are inherently generalized for ${\mu}_{s}$ . To determine how ${\mu}_{s}$ , ${\mu}_{a}$ , and $g$ affect the system performance, we analyzed a range of ${\mu}_{a}\u2215{\mu}_{s}$ and $g$ values that might be encountered when imaging tissues. We present the effects only for the pinhole and slit apertures, since these two represent the extreme cases.

Soft tissues generally have ${\mu}_{a}\u2215{\mu}_{s}<1$ ; therefore, we analyzed this ratio in the range 0 to 1 (Ref. 30). Figure 14 shows how the axial sensitivity and the lateral and axial IQR are effected by changing ${\mu}_{a}\u2215{\mu}_{s}$ . The effect on the axial sensitivity is subtle. As ${\mu}_{a}\u2215{\mu}_{s}$ increases, the focus peak slightly increases because the increased absorption causes a decrease in the tails for ${\epsilon}_{z}\cdot {\mu}_{s}>0$ . There is also a small increase in signal toward the tissue surface, although the effect is appreciable only in the slit aperture, since it is more sensitive to signal from the surface. The relatively flat IQR lines at different imaging depths indicate that both the pinhole aperture and the slit aperture are resistant to the effects of changing ${\mu}_{a}\u2215{\mu}_{s}$ .

We studied
$g$
in the range of 0.8 to 0.95 since this is the typical range for soft tissues.^{30} Figure 15 shows how the axial sensitivity and the lateral and axial IQR are affected by changing
$g$
. As was the case with variations in
${\mu}_{a}\u2215{\mu}_{s}$
, the results show that the axial sensitivity density is not very sensitive to changes in
$g$
.

## 3.4.

### Characterizing the Effect of Fluorophore Concentration

We also studied the effect of changing ${\mu}_{f}\u2215{\mu}_{s}$ in the range of weak fluorophore concentration $({\mu}_{f}\u2215{\mu}_{s}=0.1)$ to a very strong concentration $({\mu}_{f}\u2215{\mu}_{s}=2)$ . Figure 16 shows that there are changes in the axial sensitivity and the lateral and axial IQR values for high ${\mu}_{f}\u2215{\mu}_{s}$ . Increasing ${\mu}_{f}\u2215{\mu}_{s}$ causes the surface sensitivity of the slit to increase. The IQR plots indicate that the pinhole performance is stable down to $d\cdot {\mu}_{s}=0.5$ , but beyond this depth, increasing ${\mu}_{f}\u2215{\mu}_{s}$ degrades the performance.

## 4.

## Summary and Conclusions

A Monte Carlo model was developed and implemented to characterize the axial and lateral performance of parallelized fluorescence confocal systems. The results indicate that although a slit aperture offers the potential for high-speed imaging, its axial and lateral performance are degraded. When imaging at a depth of $d\cdot {\mu}_{s}=0.5$ with an NA of 0.5, a $1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ by $450\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ micron slit has almost an order of magnitude worse axial and lateral performance compared to a $1.5\text{-}\mu \mathrm{m}$ -diam pinhole aperture. Sparse parallelized apertures such as a linear pinhole array and a Nipkow aperture can be optimized to achieve fast imaging with performance comparable to a single pinhole aperture. The results also show that performance for all apertures degrades significantly for depths greater than two mean free paths.

Although the results modeled a specific set of tissue properties and NA, they provide important insight into the benefits and limitations of different confocal pinhole configurations. There is a direct trade-off between axial and lateral performance, maximum imaging depth, and parallelization. In general, if time resolution is not important, a single pinhole aperture provides the best overall performance. Increasing the pinhole density (number of pinholes in the field of view) to speed up data acquisition will limit the axial performance and the maximum imaging depth. As a general rule, one should minimize the degree of parallelization in order to limit the effect of scattering cross talk between pinholes.

The results of this study indicate that the lateral performance of a confocal system imaging in tissue is likely to be limited by scattering and not diffraction effects. This result is useful for the optical designer. For example, it would be wasteful to optimize a system for diffraction limited performance when tissue scattering imposes the overall limit on image quality.

In addition to modeling aperture effects, we also investigated how the tissue’s nonradiative absorption coefficient ${\mu}_{a}$ , scattering anisotropy $g$ factor, and fluorescence absorption coefficient ${\mu}_{f}$ affect system performance. We found that the tissue absorption coefficient and anisotropy have little effect on the system performance for typical ranges encountered in tissue. We found that the best imaging performance is achieved when the fluorescence absorption coefficient is small relative to the scattering coefficient. Since the fluorescence coefficient is directly coupled to the amount of collected signal, there is a trade-off between increasing ${\mu}_{f}$ to increase the signal to noise ratio and decreasing ${\mu}_{f}$ to improve the lateral and axial system performance. Therefore, the system designer and user of a confocal system must determine the appropriate balance between signal-to-noise ratio and the axial and lateral resolution.

## Acknowledgments

The authors would like to thank Urs Utzinger for reviewing this paper and providing helpful feedback. This work was supported by the National Institutes of Health and the Arizona Disease Control Research Commission through the following grants: NIH CA73095, NIH CA115780, and ADCRC 9711.

## References

**,” J. Biomed. Opt., 14 044030 (2009). 1083-3668 Google Scholar**

*Clinical confocal microlaparoscope for real-time**in vivo*optical biopsies**,” Cornea, 17 485 –494 (1998). 0277-3740 Google Scholar**

*Normal human corneal cell populations evaluated by**in vivo*scanning slit confocal microscopy**,” J. Opt. Soc. Am. B, 18 1695 –1700 (2001). https://doi.org/10.1364/JOSAB.18.001695 0740-3224 Google Scholar**

*Simulation of confocal microscopy through scattering media with and without time gating***,” Appl. Opt., 37 8092 –8102 (1998). https://doi.org/10.1364/AO.37.008092 0003-6935 Google Scholar**

*Monte Carlo analysis of two-photon fluorescence imaging through a scattering medium***,” J. Am. Stat. Assoc., 44 (247), 335 –341 (1949). https://doi.org/10.2307/2280232 0162-1459 Google Scholar**

*The Monte Carlo method***,” Comput. Methods Programs Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar**

*MCML-Monte Carlo modeling of light transport in multi-layered tissues***,” Med. Phys., 10 824 –830 (1983). https://doi.org/10.1118/1.595361 0094-2405 Google Scholar**

*A Monte Carlo model for the absorption and flux distributions of light in tissue***,” 102 –111 (1989). Google Scholar**

*A Monte Carlo model of light propagation in tissue***,” Opt. Express, 16 13188 –13202 (2008). https://doi.org/10.1364/OE.16.013188 1094-4087 Google Scholar**

*Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media***,” J. Opt. Soc. Am. A, 13 952 –961 (1996). https://doi.org/10.1364/JOSAA.13.000952 0740-3232 Google Scholar**

*Efficient Monte Carlo simulation of confocal microscopy in biological tissue***,” Appl. Opt., 36 6513 –6519 (1997). https://doi.org/10.1364/AO.36.006513 0003-6935 Google Scholar**

*Forward-adjoint fluorescence model: Monte Carlo integration and experimental validation***,” J. Opt. Soc. Am. A, 20 714 –727 (2003). https://doi.org/10.1364/JOSAA.20.000714 0740-3232 Google Scholar**

*Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues***,” Lasers Surg. Med., 21 166 –178 (1997). https://doi.org/10.1002/(SICI)1096-9101(1997)21:2<166::AID-LSM8>3.0.CO;2-O 0196-8092 Google Scholar**

*Propagation of fluorescent light***,” Oper. Res., 44 816 –822 (1996). https://doi.org/10.1287/opre.44.5.816 0030-364X Google Scholar**

*Combined multiple recursive random number generators***,” Ann. Astrophys., 3 117 –137 (1940). 0365-0499 Google Scholar**

*Diffuse radiation in the galaxy***,” Lasers Life Sci., 1 309 –333 (1987). 0886-0467 Google Scholar**

*Angular dependence of HeNe laser light scattering by human dermis***,” Opt. Acta, 2 153 –162 (1955). https://doi.org/10.1080/713821040 0030-3909 Google Scholar**

*Reflection by semi-infinite diffusers***,” IEEE J. Quantum Electron., 26 2166 –2185 (1990). https://doi.org/10.1109/3.64354 0018-9197 Google Scholar**

*A review of the optical properties of biological tissues*