Translator Disclaimer
1 July 2009 Monte Carlo characterization of parallelized fluorescence confocal systems imaging in turbid media
Author Affiliations +
We characterize and compare the axial and lateral performance of fluorescence confocal systems imaging in turbid media. The aperture configurations studied are a single pinhole, a slit, a Nipkow disk, and a linear array of pinholes. Systems with parallelized apertures are used clinically because they enable high-speed and real-time imaging. Understanding how they perform in highly scattering tissue is important. A Monte Carlo model was developed to characterize parallelized system performance in a scattering media representative of human tissues. The results indicate that a slit aperture has degraded performance, both laterally and axially. In contrast, the analysis reveals that multipinhole apertures such as a Nipkow disk or a linear pinhole array can achieve performance nearly equivalent to a single pinhole aperture. The optimal aperture spacing for the multipinhole apertures was determined for a specific tissue model. In addition to comparing aperture configurations, the effects of tissue nonradiative absorption, scattering anisotropy, and fluorophore concentration on lateral and axial performance of confocal systems were studied.



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.

Fig. 1

Four types of confocal apertures. The dashed arrows show the scanning patterns required to cover equivalent fields of view.


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.

Fig. 2

Simplified optical layout of a fluorescence confocal system. Cones illustrate nonscattering system illumination and collection beam paths. The dashed line depicts an excitation photon that scatters outside the typical beam path, fluoresces, and couples back into the confocal system.


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.


Monte Carlo Model of a Fluorescence Confocal System

A direct Monte Carlo7 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 schemes11, 12 and multiple independent run methods13, 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.


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 known4 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 r0 . Then the photon is back-projected to where it encounters the tissue surface at r0 . At this location, it is launched back into the turbid medium and scattered until it is absorbed by the fluorophore and fluoresces at rs . The fluorescent photon’s new direction is randomly selected from an isotropic distribution and then scattered until it exits the tissue at rd . To determine if the photon is collected by the system, it is back-projected to rd . Since rd falls within the aperture’s opening and the photon’s angle is within the maximum collection angle ϕmax as determined by the NA, the photon is recorded as collected signal. Ideally, in a nonscattering media with an infinitely small pinhole, rs will equal rd .

Fig. 3

Example of a photon that scatters, fluoresces, scatters, and is collected by the confocal system with a single pinhole aperture.


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 rci . Figures 4a and 4b depict how rci is determined for pinhole-based apertures. We assume that each pinhole has one detector behind it, resulting in rd being recorded as the center of the i ’th pinhole rci that the photon is collected in. Since a slit aperture has many detectors [Fig. 4c], rci is the center of the detector element that contains rd . 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 rc1 .

Fig. 4

Illustration showing the final collected position rci of a detected photon incident in an aperture at rd . For a single pinhole aperture (a), all detected photons are recorded with a collected position at the center of the aperture rc1 . If a detected photon is inside the fifth pinhole of a linear array aperture (b), its final collected position is rc5 . If a detected photon is inside the fourth detector of a slit aperture (c), its final collected position is rc4 .


There is an error in the collected signal as defined by the vector ε=rcirs . The distribution of εz characterizes the axial sensitivity of the confocal system, and the bivariate distribution of εx and εy characterizes the lateral blur induced by scattering. Once rd is shifted to the central point in the detection element receiving the photon, the bivariate distribution of εx and εy also accounts for the detector sampling.


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 generator18 with a period of 2185 (about 1056 ) was used. To generate results in a timely manner, the code was run on a 40-GHz distributed Xgrid cluster. In the following text, each step is discussed in detail:

Fig. 5

Flow diagram for the Monte Carlo simulation of photons through a fluorescence confocal system imaging in turbid media.


Step 1. By mapping the uniformly illuminated confocal aperture to the tissue focal plane (z=0) , the initial spatial position of the photon r0=(x0,y0,z0=0) can be randomly generated.

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

Eq. 1

Then F(x) is inverted to yield the inverse cumulative probability density function F1(x) . Last, the random variable X can be generated as

Eq. 2

Note that in the following sections, each instance of ξ 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 r0=(x0,y0,z0=0) coordinates for the photon in the aperture image can be generated for a single pinhole with radius a using


where a is a random radial distance from the aperture center, and θ is a random angle. To uniformly generate photons in x and y , the random radial position a must have an increasing probability away from the center. In the case of a linear pinhole array or a Nipkow aperture, first a pinhole in the aperture is randomly selected with equal probability for all pinholes. Then equation set (3) is used to find the random offset position for that pinhole. For the slit aperture, x0 and y0 are uniformly distributed within the aperture width and height.

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-y azimuthal plane and maximally depart from the z axis by the zenith angle ϕmax=sin1(NAn) , 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 θ and ϕ , which are the zenith angle, and azimuthal angle, respectively. The azimuthal angle θ has a uniform distribution that ranges from 0 to 2π . The zenith angle requires a uniform sampling of cosϕ for phi ranging from 0 to ϕmax in order to achieve uniform sampling within the NA.

Step 2. This step consists of back-projecting the photon from r0 to r0=(x,y,z=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


where direction cosines α , β , and γ are given by


and the propagation distance is given by

Eq. 6

In this case, we are propagating from z=0 to the surface lying a distance d above, and therefore Δz=d . Since scattering will begin in the tissue after this point, we need not consider refraction at the tissue surface.

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 μs [cm1] , μa [cm1] , 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 μf=ΓC [cm1] , which can be expressed in terms of the fluorescence extinction coefficient Γ [cm1Lmol] (assuming a quantum efficiency of one) and the fluorophore concentration C [molL] .

The mean free path a photon travels before it scatters is 1μs . The distribution of path lengths Δs that photons travel before scattering is described by the Beer-Lambert law probability density function;

Eq. 7

To randomly generate an individual path length Δs , Eqs. 1, 2 are used to define the generator function for Δs ,

Eq. 8

Once Δs is generated, the photon is propagated using Eq. 4 to the position at which the scattering event occurs. The scattering event causes the photon to change its trajectory. The azimuthal direction change is described by the uniformly distributed random variable Δθ ranging from 0 to 2π. The change in zenith angle Δϕ is modeled using the Henyey-Greenstein phase function20

Eq. 9

The random generator computed from Eqs. 1, 2 for cosΔϕ is

Eq. 10

Jacques 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 cosΔϕ is uniformly distributed between 1 and 1.

Given the scattered angle (Δθ,Δϕ) and the propagation direction (α,β,γ) of the incoming photon, the new direction (α,β,γ) is given by



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 ξR . The reflection coefficient R is the average of the Fresnel reflection coefficients for the s and p polarization states:

Eq. 12

where ϕi=cos1γ is the angle of incidence, n1 is the tissue index, and n2 is the index outside of the tissue. For ϕi greater than the critical angle, the photons always reflect back into the tissue. If the photon exits the surface, it is lost, since the system rejects nonfluorescence signal.

To determine where the excitation photon terminates in the tissue, the path length to fluorescence, lf , and the path length to absorption, la , are determined in the same manner that Δs was determined using the random generator

Eq. 13


Eq. 14

Note that lf and la are determined only once per photon. If la<lf , 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 lf<la , the photon will be absorbed and reemitted as fluorescence when its total path length L is equal to the fluorescence path length lf . The total path length L is the accumulated path length for the number of scattering events

Eq. 15


Photons with lf<la travel until L=lf , where they are converted to excited fluorescence photons in the next step. If the next propagation step size Δs causes the total path length to exceed lf , the propagation terminates partway through that path at the point where the total path length equals lf .

Step 4: When the photon has traveled a distance L=lf , it is absorbed and emitted as fluorescent energy. The position where fluorescence occurs represents the location of the signal and is recorded as rs . Since we assume that fluorescent emission photons are isotropically generated from the fluorophore, a new photon direction (Δθ,Δϕ) 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 Δs . The propagation continues until the photon exits the tissue surface or the photon is absorbed when L=la using a newly generated la .

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 rd on the surface to the point rd 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 rci 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 μaμ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 Nc are collected by the system.


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 μaμs with matched boundary conditions (n=n) . 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 results24 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 Saloma6 for a fluorescence confocal system with a pinhole aperture. Since Blanca and Saloma did not report μ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 Nc 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 Nc 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 Nc 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-μm -square region. The pinhole aperture was 1.5μm in diameter. The slit aperture was 1.5μm by 450μm , with the long dimension oriented along y . Inside the slit, three hundred 1.5-μm -square detectors recorded signal in parallel. The linear array consisted of a line of pinhole apertures 1.5μm in diameter spaced 30μm center-to-center, spanning 450μm (15 pinholes total). For the Nipkow aperture, a 2-D array of sixty 1.5μm apertures were spaced 60μ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μ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 is3

Eq. 16

For a confocal system with an NA of 0.5, the optimal pinhole diameter would be 0.633μm , assuming λ is 632.8nm . The pinhole diameter in the modeled system is larger because better signal efficiency is achieved with radii greater than that given by Eq. 16.

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 μaμs , g , and μfus 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.8nm with values μs(1g)=12cm1 , μf=40cm1 , and μa=0.4cm1 (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 μs=80cm1 . We also assume that the average fluorophore concentration is uniform throughout the tissue with μf=0.5μs . Our analysis used a nominal imaging depth of d=62.5μ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μs . Similarly, we can generalize the axial error from focus εz as the average number of mean free paths from focus εzμ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.


Characterizing Aperture Imaging Performance

Videos 1 shows the 3-D distributions28 of the collected fluorescence photon signal positions rs for the four aperture configurations imaging in simulated tissue at an imaging depth of dμs=0.5 mean free paths (62.5μm) . In each of the figures, the tissue surface lies in the x-y plane at z=62.5μ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.

Video 1

Movies of the 3-D distribution of collected photon signal positions rs for each aperture configuration in simulated tissue. Tissue surface is at z=62.5μm (grid plane); focus is at z=0 (QuickTime 5 MB). .


Figure 6 shows the projection29 of the collected fluorescence photon signal rs onto the x-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.

Fig. 6

Lateral x-y signal distribution for rs for each aperture configuration in esophagus tissue. Note log color map scale.


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 ε=rcirs , as shown in Videos 2 . In an ideal confocal system, |ε| should be zero for all collected fluorescent photons. As the distribution of ε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.

Video 2

Data from Fig. 6 shown in terms of error in collected position ε for each aperture configuration in esophagus tissue (QuickTime 3 MB). .


Compared 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 ε on the axial planes εzεx and εzεy and the lateral plane εxεy for the four aperture configurations. Since the majority of the energy is localized near ε=0 , the plots are shown with a log color scale.

Fig. 7

Bivariate projections of the data shown in Videos 2 in the axial planes εzεx (left column) and εzεy (middle column) and the lateral plane εxεy (right column) for each aperture configuration. Tissue surface is at z=62.5μm , focus is at z=0 ; dμs=0.5 .


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 ε=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 εzεy plane. The interesting V-pattern in the εzε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 ε=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 εzμs (unit-less) away from focus. The nominal focus is at εzμs=0 , and the tissue surface is at εzμ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.

Fig. 8

Collected photon signal density broken down into ballistic and scattered densities for each confocal aperture. Focus is at εzμs=0 , and the tissue surface is at εzμs=0.5 .


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μ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 Nc 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μs=1 . The linear array and Nipkow maintain at least 90% ballistic signal down to dμs=0.5 but drop off thereafter. The drop off after dμs=0.5 is a result of the optimization for imaging down to dμs=0.5 . The figure highlights that the slit aperture has the worst performance, limited to 71% at best for dμs=0.25 .

Fig. 9

Comparison of the ratio of the number of collected ballistic photons to the total number of collected photons as a function of imaging depth dμs for each of the apertures is shown in (a). The standard deviation of 10 runs for each point in (a) is shown as an estimate of the Monte Carlo model error in (b).


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 ε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 ε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 εzμ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.

Fig. 10

Axial density at three imaging depths for each of the aperture configurations.


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 ε and the Airy IQR. To find the Airy IQR, the intensity of the Airy pattern (Fraunhofer diffraction pattern)

Eq. 17

can be integrated to give the total power contained in a radius ρ

Eq. 18

J0 and J1 are the zeroth and first-order Bessel functions, and k=2πλ is the wave number. Solving for the radius that contains the middle 50% of the power, we find that the Airy IQR is twice this radius, or 0.535λNA . At λ=632.8nm and an NA of 0.5, the Airy IQR is 0.662μm .

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μm and 8.72μm , respectively. The central lobe of the axial response can be reasonably approximated with a Gaussian profile. Since the standard deviation σ for a Gaussian function in terms of the FWHM is

Eq. 19

and the IQR is ±0.674σ , the diffraction-limited axial IQR for the modeled pinhole and slit systems are 3.52μm and 4.99μm , respectively.

Figure 11a presents lateral and axial IQR values as a function of imaging depths dμs . To characterize the lateral distribution of ε , we compute the IQR of the signed radial response ερ=(εx2+εy2)12sign(εx)sign(εy) . 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 Nc is sufficiently large to reliably estimate the IQR.

Fig. 11

Lateral and axial IQR performance (in μm) as a function of imaging depth dμs for each aperture is shown in (a). The standard deviation of 10 runs for each point in (a) is shown as an estimate of the Monte Carlo model error in (b).


Figure 11 shows how the pinhole aperture has a stable axial performance of about 4μm and lateral performance of about 1μm down to dμs=1 . The linear array and Nipkow aperture maintain performance comparable to the pinhole aperture down to their optimized depth dμs=0.5 (optimization discussed in Sec. 3.2). The slit aperture has substantially degraded performance. All apertures fail for dμ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.


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μs must be specified. Since the axial and lateral performance degrades as the imaging depth is increased, the pinhole spacing δ should be optimized to obtain the minimum acceptable performance at the maximum required imaging depth. To stay consistent with our previous discussion, δ 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 δ at dμs=0.5 and dμs=1 , respectively. As δ 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 δ increases. For the target imaging depth of dμs=0.5 , the axial and lateral performance do not substantially improve beyond δ=30μm . Thus, for the system modeled, an optimized linear array would have 15 pinholes spaced 30microns apart in a 450-micron line.

Fig. 12

Linear array aperture pinhole spacing δ optimization. No substantial improvement for δ> 30μm at target depth dμs=0.5 .


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 δ 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 δ greater than 60μm Therefore, for the system modeled, a Nipkow disk with 60 pinholes spaced 60μm apart in a 450-μm -square area is nearly optimal.

Fig. 13

Nipkow aperture pinhole spacing δ optimization. No substantial improvement for δ> 60μm at target depth dμs=0.5 .



Characterizing Effects of Tissue Properties

Changes in μs , μa , and g can have substantial effects on the axial and lateral performance. Since we present our results in terms of εzμs and dμs , they are inherently generalized for μs . To determine how μs , μa , and g affect the system performance, we analyzed a range of μaμ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 μaμ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 μaμs . The effect on the axial sensitivity is subtle. As μaμs increases, the focus peak slightly increases because the increased absorption causes a decrease in the tails for εzμ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 μaμs .

Fig. 14

Effect on the lateral and axial performance with changing μa for the pinhole and slit apertures.


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 μaμs , the results show that the axial sensitivity density is not very sensitive to changes in g .

Fig. 15

Effect on the lateral and axial performance with changing g for the pinhole and slit apertures.



Characterizing the Effect of Fluorophore Concentration

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

Fig. 16

Effect on the lateral and axial performance with changing μf for the pinhole and slit apertures.



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μs=0.5 with an NA of 0.5, a 1.5μm by 450μm micron slit has almost an order of magnitude worse axial and lateral performance compared to a 1.5-μ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 μa , scattering anisotropy g factor, and fluorescence absorption coefficient μ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 μf to increase the signal to noise ratio and decreasing μ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.


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.



A. A. Tanbakuchi, A. R. Rouse, J. A. Udovich, K. D. Hatch, and A. F. Gmitro, “Clinical confocal microlaparoscope for real-time in vivo optical biopsies,” J. Biomed. Opt., 14 044030 (2009). 1083-3668 Google Scholar


R. Mustonen, M. McDonald, S. Srivannaboon, A. Tan, M. Doubrava, and C. Kim, “Normal human corneal cell populations evaluated by in vivo scanning slit confocal microscopy,” Cornea, 17 485 –494 (1998). 0277-3740 Google Scholar


T. R. Corle and G. S. Kino, Confocal Scanning Optical Microscopy and Related Imaging Systems, Academic Press, San Diego (1996). Google Scholar


Confocal Microscopy, Academic Press, London (1990). Google Scholar


M. Magnor, P. Dorn, and W. Rudolph, “Simulation of confocal microscopy through scattering media with and without time gating,” J. Opt. Soc. Am. B, 18 1695 –1700 (2001). 0740-3224 Google Scholar


C. Blanca and C. Saloma, “Monte Carlo analysis of two-photon fluorescence imaging through a scattering medium,” Appl. Opt., 37 8092 –8102 (1998). 0003-6935 Google Scholar


N. Metropolis and S. Ulam, “The Monte Carlo method,” J. Am. Stat. Assoc., 44 (247), 335 –341 (1949). 0162-1459 Google Scholar


L. Wang, S. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed., 47 131 –146 (1995). 0169-2607 Google Scholar


B. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys., 10 824 –830 (1983). 0094-2405 Google Scholar


S. Prahl, M. Keijzer, S. Jacques, and A. Welch, “A Monte Carlo model of light propagation in tissue,” 102 –111 (1989). Google Scholar


A. Liebert, H. Wabnitz, N. Zolek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express, 16 13188 –13202 (2008). 1094-4087 Google Scholar


J. Schmitt and K. Ben-Letaiet, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A, 13 952 –961 (1996). 0740-3232 Google Scholar


R. Crilly, W. Cheong, B. Wilson, and J. Spears, “Forward-adjoint fluorescence model: Monte Carlo integration and experimental validation,” Appl. Opt., 36 6513 –6519 (1997). 0003-6935 Google Scholar


J. Swartling, A. Pifferi, A. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A, 20 714 –727 (2003). 0740-3232 Google Scholar


A. Welch, C. Gardner, R. Richards-Kortum, E. Chan, G. Criswell, J. Pfefer, and S. Warren, “Propagation of fluorescent light,” Lasers Surg. Med., 21 166 –178 (1997).<166::AID-LSM8>3.0.CO;2-O 0196-8092 Google Scholar


J. Schmitt, A. Knuettel, and M. Yadlowsky, J. Opt. Soc. Am. A, 11 2226 –2235 (1994). 0740-3232 Google Scholar


M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, 2Network Theory Ltd, Bristol, UK (2006). Google Scholar


P. L’Ecuyer, “Combined multiple recursive random number generators,” Oper. Res., 44 816 –822 (1996). 0030-364X Google Scholar


B. R. Frieden, Probability, Statistical Optics, and Data Testing: A Problem-Solving Approach, 10 3rd ed.Springer, Berlin (2001). Google Scholar


L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Ann. Astrophys., 3 117 –137 (1940). 0365-0499 Google Scholar


S. Jacques, C. Alter, and S. Prahl, “Angular dependence of HeNe laser light scattering by human dermis,” Lasers Life Sci., 1 309 –333 (1987). 0886-0467 Google Scholar


A. Welch and M. van Gemert, Optical-Thermal Response of Laser-Irradiated Tissue, Plenum Press, New York (1995). Google Scholar


H. van de Hulst, Multiple Light Scattering: Tables, Formulas, and Applications, Academic Press, New York (1980). Google Scholar


R. Giovanelli, “Reflection by semi-infinite diffusers,” Opt. Acta, 2 153 –162 (1955). 0030-3909 Google Scholar


R Development Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria (2008). Google Scholar


W. Cheong, S. Prahl, and A. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron., 26 2166 –2185 (1990). 0018-9197 Google Scholar


Biomedical Photonics Handbook, CRC Press, Boca Raton, FL (2003). Google Scholar


D. Adler and D. Murdoch, rgl: 3D Visualization Device System (OpenGL), (2008) Google Scholar


D. Carr, Nicholas Lewin-Koh and M. Maechler, hexbin: Hexagonal Binning Routines, (2008) Google Scholar


B. Wilson and S. Jacques, IEEE J. Quantum Electron., 26 2186 –2199 (1990). 0018-9197 Google Scholar
©(2009) Society of Photo-Optical Instrumentation Engineers (SPIE)
Anthony A. Tanbakuchi, Andrew R. Rouse, and Arthur F. Gmitro "Monte Carlo characterization of parallelized fluorescence confocal systems imaging in turbid media," Journal of Biomedical Optics 14(4), 044024 (1 July 2009).
Published: 1 July 2009

Back to Top