2.1.

Theoretical DH-PSF

In a wavefront coded microscope system, a phase mask with phase $\varphi (f_x,f_y)$ is placed at the back focal plane of the imaging lens which modifies the appearance of the traditional or clear aperture (CA) PSF (Fig. 2). The DH phase mask design shown by the wrapped phase values in Fig. 2(a) (leftmost image), is used to generate the DH-PSFs that encodes depth information in the spatial rotation of its main lobes [Fig. 1(c)] thereby characterizing a system with new properties. The DH phase mask was developed from the idea of wave propagation with rotating intensity distribution suggested by Piestun et al.²³ which results in a high efficiency rotating PSF,²⁴ named later as the DH-PSF because it appears as a DH along the $Z$-axis of the microscope [Fig. 2(a), rightmost image].

Fig. 2

Comparison of (a) DH-PSF and (b) CA-PSF simulated images shown in $XY$, $XZ$, and $YZ$ central cut-views from left to right, respectively. The maximum intensity projection (MIP) of the $YZ$ cut-view makes the double-helix appearance more evident. Images are displayed on a $256\times 256$ pixel grid using a different range of intensities in each case to allow visualization of structural details. Lens: $63\times /1.4\mathrm{NA}$, oil immersion. Scale bar: 3.2 μm.

Just as in the case of the CA-PSF, the DH-PSF varies with both the imaging depth $z=z_i$ (at which the microscope is focused) and the location of the light point source $z=z_o$:

2.2.

Image Formation in the Presence of SA

The intensity in the 3-D image $g({\mathbf{x}}_i)$ formed by a DH-based system in the presence of depth-variant imaging conditions can be written as

Using this strata-based linear interpolation, the image of the $m$’th stratum can be obtained by

3.1.

DH-PSF Computation

The Gibson and Lanni model³ in the COSMOS PSF module was used to compute depth variant complex-valued coherent clear aperture PSFs (CA-PSFs) for various amounts of SA (due to point sources at different depths under the coverslip either within water or glycerol) for oil-immersion lenses with various NA (1.1 to 1.3) and various wavelengths (514, 488, and 633 nm). The imaging conditions used in each case were chosen to match the imaging conditions of the experimentally determined DH-PSFs. In the case of the experimental DH-PSFs from a fabricated mask with diameter $D=2.7\mathrm{mm}$, the effective NA of the system was found to be 1.28 computed using the following formula:

Using a Matlab script, the 2-D Fourier transform of each $XY$-plane in a 3-D complex-valued CA-PSF was computed to produce the generalized pupil function. The DH phase mask design was then applied to form a modified generalized pupil function and finally, the square magnitude of its inverse Fourier transform was computed to obtain a 2-D DH-PSF with SA [Eq. (1)]. This process was repeated for all the planes in the 3-D CA-PSF resulting in a 3-D DH-PSF.

3.2.

DH-PSF Depth Variability

To quantify the depth variability of the DH-PSF due to SA, a normalized mean square error (NMSE) that quantifies the difference between depth-variant 3-D DH-PSFs (due to point sources at different depths) was used:

3.3.

Experimentally Determined DH-PSF

In order to characterize the 3-D DH-PSF and its depth-variability due to SA, experimentally determined DH-PSFs were acquired from two systems (configured on an optical table at Dr. R. Piestun’s laboratory, University of Colorado, Boulder) which incorporated the DH phase mask in the back focal plane of the objective lens using a 4F system and either: (1) a fabricated phase mask;²⁶ or (2) a spatial-light modulator implementation of the mask.²⁷

In the first system with the fabricated mask, a magnification of $100\times$, a $\mathrm{NA}=1.49$, a wavelength of 488 nm, and an oil-immersion lens were used. DH-PSFs were imaged at two different depths from a sample with quantum dots, 100 nm in diameter (QD-050/515, Invitrogen Inc.), fixed to the coverslip and to the glass slide with a 90% glycerol solution ($\mathrm{RI}=1.47$) in between them. The experimental data consisted of images acquired at 100 different axial positions with 8 images at every axial location and a $z$-step size of 50 nm. The images of the quantum dots fixed to the glass slide below the glycerol solution were measured at a depth of 13.5 μm below the coverslip inferred from the movement of the piezo-electric stage.

In the second system (with the SLM implementation of the mask), DH-PSF images were acquired from a slide containing quantum dots (QD-525, Invitrogen Inc.) mounted in a 90% glycerol solution ($\mathrm{RI}=1.47$) using an emission wavelength of 514 nm, an oil immersion lens ($\mathrm{RI}=1.512$) with a theoretical $\mathrm{NA}=1.3$, and an effective magnification of $91.2\times$. The experimental data consisted of images acquired at 86 different axial positions with 8 images at every axial location and $z$-step size of 50 nm.

3.4.

Simulated Images of a 3-D Object for Particle Localization

A computer-generated object was chosen to investigate particle localization in the presence of SA. The object comprised of three beads, with a 50 nm diameter, located at depths 0, 1.5, and 3.0 μm, respectively, within water [Fig. 3(a)]. 3-D simulated images were computed using the COSMOS Tools module for a DH system without SA [Fig. 3(b) and 3(c)] and with SA [Fig. 3(d)]. The simulated images without SA were computed using a single nonaberrant DH-PSF. To simulate depth-variant (DV) imaging due to SA, the strata-based forward-image model [Eq. (3)] was used with multiple spherically aberrant DH-PSFs defined at depths 0, 0.8, 1.6, 2.4, and 3.2 μm.

Fig. 3

Localization of multiple point sources using a single 2-D image from a DH-based system. $XZ$ central cut-view: (a) image of a three bead object with beads located at depths: 0, 1.5 and 3.0 μm; (b) image from the forward model of a DH-PSF-based system with a 3 bead structure. $XY$ central cut-view from the 3-D image of this object simulated with a DH-based system computed: (c) without SA; and (d) with SA. Dashed lines in (b) and (c) show pairs of DH-PSF lobes associated with each bead. The rotation angle associated with each bead is computed and mapped to the calibration curve (e), where the rotation angle versus axial distance (defocus) from a DH-PSF without any SA is plotted. The circles in (e) indicate the locations of the three beads. Scale bars: 1.6 μm in (a); 3.2 μm in (b), and 0.8 μm in (c) and (d). Lens: $63\times /1.4\mathrm{NA}$ oil immersion. Wavelength: 633 nm.

Particle localization in a nonaberrant system is also shown in Fig. 3(e) where the depths of the three beads in the simple object are obtained from the calibration curve through the rotation angles determined from the simulated image of the object. In the case without SA, the correct depths for the three beads (0, 1.5, and 3.0 μm, correspondingly) were determined. In the presence of SA, the observed change in the rotation angles that can be determined by the lines indicated in Fig. 3(c) and Fig. 3(d) can lead to errors. Particle localization in the presence of SA is discussed in detail in Sec. 4.2.1.

3.5.

Calibration Curve

The angle subtended by the primary lobes of the DH-PSF was computed using an angle estimation algorithm²⁸ (implemented in a Matlab) and it was plotted against axial distance to generate what is known as the calibration curve [Fig. 3(e)]. The calibration curve can be used to obtain the exact depth of point sources within an object by matching rotation angle to axial distance. Historically, the 0 deg angle associated with the best focus has been mapped to the 0 μm reference plane. However, in 3-D samples where different structures appear in focus at different axial positions, it may be more convenient to associate the 0 μm reference plane to the $+90\mathrm{deg}$ angle, so that other axial locations can be calculated relative to this point. We have used this convention throughout the paper.

Several calibration curves were generated for different lenses under various imaging conditions. A small pixel size of 50 nm was used for the DH-PSF computation because this affects the accuracy of the calibration curve as observed in our studies. Better axial resolution gives more points for plotting the calibration curve which improves the accuracy of the calibration curve.

4.1.

Effect of SA on the DH-PSF

As it is known, with increasing SA, an asymmetry in the PSF along the $z$-axis as well as an apparent shift in the peak intensity of the CA-PSF are observed [Fig. 4(a)], resulting in the difference between the apparent and the actual location of a point light source. In the DH-PSF case, this axial shift is evident in the $Z$ location of the focal plane where the two lobes are aligned along the horizontal axis, i.e., the 0 deg angle. Although this focal plane is located at $z=0.0\mu \mathrm{m}$ when SA is not present [Fig. 1(c)] it appears at approximately $z=0.8\mu \mathrm{m}$ in the presence of SA induced due to the point source being at a depth of 50 μm below the coverslip in water [Fig. 1(d)]. This depth variability of the 3-D DH-PSFs is clearly evident in the comparison of $XZ$ central cut-view images from DH-PSFs computed at different depths [Fig. 4(c)].

4.1.1.

Depth variability

Fig. 4 summarizes results that quantify the depth variability computed with Eq. (5) [Fig. 4(a)] from 3-D PSFs for a CA and a DH system at different point source depths, i.e., with a varying amount of SA [Fig. 4(b) and 4(c)]. As evident in Fig. 4(a), the NMSE for the DH-PSF changes less rapidly than the NMSE for the CA-PSFs as the axial location (depth) of the point light source in water increases, causing the SA to increase (the slight drop in the NMSE computed for the DH-PSF at depth 100 μm is an artifact due to the spread and truncation of the DH-PSF). The NMSE also suggests that for the imaging conditions simulated, the DH-PSF does not change for depths larger than $\sim 50\mu \mathrm{m}$. This result could be exploited in COSM that accounts for depth-variant imaging of thick samples using multiple depth-variant PSFs. As mentioned earlier, a reduction in the number of PSFs used in the computations reduces the time required to compute the fluorescence concentration in the sample.

Fig. 4

PSF depth variability due to SA quantified by a NMSE metric: (a) the 3-D NMSE computed for the DH-PSF and the CA-PSF. $XZ$ central cut-view images from 3-D PSFs of a 1.3 NA, oil lens with an effective magnification of 91.2X for: (b) a CA; and (c) an aperture with a DH mask. PSFs were computed for point light sources located at different depths (0 to 100 μm) below the coverslip in water. Scale bar: 6.4 μm.

4.1.2.

Calibration curve

Fig. 5 summarizes the impact of both the numerical aperture (NA) of the lens and SA on the calibration curve generated for a DH-mask-based system by comparing calibration curves of DH-PSFs without SA (i.e., due to a point source at 0 μm depth under the coverslip) for lenses with a different NA [Fig. 5(a)] and calibration curves of the DH-PSFs of a single lens (with $\mathrm{NA}=1.2$) with increasing amounts of SA [Fig. 5(b)]. In order to show the change in the range of the DH-PSF rotation (i.e., the distance covered along $Z$ as the two primary lobes rotate from $-90$ to 90 deg) as the NA changes, we aligned the 3 curves shown in Fig. 5(a) at the best focus location (i.e., at the 0 deg angle). The 90 deg angle was mapped to the 0 μm axial position only for the calibration curve for $\mathrm{NA}=1.1$ in order to observe the maximum difference in the curves as the rotation angle changes in either direction from 0 deg. As it is evident in Fig. 5(a), a reduction in the NA of the lens introduces a symmetric increase in the range of rotation on either side of the 0 deg angle, i.e., the best focus. In fact, the range of rotation is proportional to $\sim \lambda /\mathrm{NA}$ where $\lambda$ is the wavelength of light.

Fig. 5

Effect of NA and SA on the calibration curve of theoretically determined DH-PSFs with an effective magnification of $91.2X$; oil-immersion lens. Curves for DH-PSFs with: (a) different NAs and zero SA; and (b) $\mathrm{NA}=1.2$ and different amounts of SA due to changing the location of the point source, i.e., the depth $d$ within the sample ($d=0\mu \mathrm{m}$ corresponds to the zero-SA case). Wavelength: 514 nm. RI of sample: 1.33.

However, symmetry around the 0 deg angle position is lost in the calibration curves of DH-PSFs in the presence of SA shown for the system with $\mathrm{NA}=1.2$ [Fig. 5(b)]. SA causes the range of rotation to increase only on one side of the best focus [denoted by the 0 deg. angle in Fig. 5(b)] resulting in an asymmetric calibration curve with respect to the 0 deg angle. As observed in Fig. 5(b), when the axial position moves away from the best focus towards the coverslip, the calibration curves with different amounts of SA show reasonable agreement; however, when the axial position moves deeper into the sample, the calibration curves show a significant disagreement. This is manifested by the curves converging at $+90\mathrm{deg}$ and completely diverging at $-90\mathrm{deg}$. The change in the calibration curve due to SA impacts traditional particle localization methods that are based on the calibration curve.

4.2.

Particle Localization in the Presence of SA

The accuracy of article localization based on the DH-PSF calibration curve is affected by the accuracy of the calibration curve. Figure 6(a) shows an apparent shift of the system’s best focus due to the presence of SA. For a nonaberrant system, the 0-deg angle point on the calibration curve represents the best focus of the system while for a system with SA, it represents the axial location of the circle of least confusion. In Fig. 6(b), the difference (or shift) in the axial distance evident in Fig. 6(a), is plotted versus the true axial location of the emitter for different rotation angles (which correspond to different amounts of defocus represented by the rotation angle).

Fig. 6

Impact of SA on the DH-PSF calibration curve. (a) Calibration curves for DH-PSFs with different amount of SA due to a point light source located at different depths within water. (b) The shift in the emitter’s axial location versus the true axial location (i.e., depth within water) of the emitter for different rotation angles (in the range of $-90$ to 90 deg) for the axial location of the emitters changing from 0 to 100 μm in steps of 10 μm The inset shows a view of the observed shift for small changes in the axial locations of the emitter in the range of 0 to 3.2 μm using 5 locations from 0.0 to 3.2 μm, in steps of 0.8 μm Lens: $63\times /1.4\mathrm{NA}$, oil immersion. Wavelength: 633 nm.

As shown in Fig. 6(b), the shift in the calibration curve changes with defocus and it is nonlinearly related with the true axial location of the emitter. For increasing SA, the apparent best focus of the emitter changes as we go deeper into the sample and thus accurate localization is not possible without the use of the calibration curve specific to the depth where the point source is located. It is interesting to note that the rate of change in the apparent location of the emitter changes through the entire range of rotation [Fig. 6(b)].

In summary, Fig. 6 shows that in the presence of SA, an error in localization occurs due to the apparent shift in the true location of a point source. In a thick 3-D sample, with multiple point sources at different depths throughout the sample, a single calibration cannot be used without incurring a particle localization error. In order to overcome this, all the calibration curves shown in Fig. 6(a) must be used. Here we propose the use of only two calibration curves: the traditional calibration curve from a nonaberrant DH-PSF and a second curve which provides a map of the apparent depth to the true depth. Such a calibration curve can be constructed from the curves in Fig. 6(a) where the depths in the legend of the figure represent the true axial location of the emitter and the apparent depth that corresponds to each of these is obtained by adding the observed shift of the “best focus” to the true depth. For small amounts of SA, this new curve is linear with a slope given by the RI ratio of the sample and the immersion medium of the lens:

This relation is valid only for small amounts of SA as previously reported in a confocal microscopy study.²⁹ We have confirmed this relation with a particle localization simulation described in the next section using the 3 bead object (discussed in Sec. 3.4).

4.2.1.

Particle localization in a three-bead object

For particle localization in the presence of SA using the angle estimation and calibration curve method, a single $XY$ section image taken from a 3-D DH-system image was used [Fig. 7(a)]. The angle subtended by the two lobes, associated with each particle in the object, was estimated. As expected, these angles mapped to the traditional calibration curve of a nonaberrant DH-PSF [Fig. 7(b)] results in incorrect depths for the three beads (0, 1.7, and 3.3 μm, correspondingly). This change in precision is due to the fact that the calibration curve used does not account for the SA present in the image. This error is compensated using a second calibration curve which maps the apparent depth to the true depth [Fig. 7(c)], constructed using the axial location shift observed in 0 deg angle position of the calibrations curves with SA in Fig. 6(a) as explained above.

Fig. 7

Localization of particles in the presence of SA. (a) $XY$ view of the forward DH-PSF image of a 3 bead object in the presence of SA. Proposed calibration curves for a system in the presence of SA: (b) rotation angle versus axial distance (defocus) from a DH-PSF without any SA (the circles indicate the apparent locations of the three beads); and (c) true axial location of emitter versus apparent depth obtained by plotting the shift observed in the best focus location from Fig. 6(a). Scale bar: 0.8 μm. Lens: $63\times /1.4\mathrm{NA}$ oil immersion. Wavelength: 633 nm.

As the amount of SA in this study is relatively small, the apparent depth is the sum of the true depth and the shift resulting in a linear curve with a slope equal to the RI ratio in Eq. (6). However, this relation does not hold for large amounts of SA as evident by the variability of the shift in the axial position of an emitter with increasing amount of SA [Fig. 6(b)].

4.3.

Comparison of Theoretical and Experimental DH-PSFs

A comparison between experimental and theoretical DH-PSFs is summarized in this section. The calibration curve for the experimental data is prepared using the same method as the one used for the simulated data described in Sec. 3.5.

The first comparison study was conducted to assess the impact of SA on experimental DH-PSFs imaged from quantum dots placed on the coverslip and glass slide, respectively, as described in Sec. 3.3. Calibration curves computed from experimental DH-PSFs were compared to calibration curves from simulated DH-PSFs computed at depths 0 and 13.5 μm using the experimental imaging conditions. Four isolated quantum dots were chosen from each data set [Fig. 8(a) and 8(b)] in order to compute calibration curves for the comparison study [Fig. 8(c)]. The observed shift at the 0 deg angle between the calibration curves from the coverslip and glass slide quantum dots is determined to be 250 nm. Shift in the axial location of the beads is predicted by Eq. (6) and reinforced by the length and shift in the 0 deg angle location of the simulated DH-PSF calibration curves. The similarity in the trends observed in the simulated and experimental calibration curves of Fig. 8 is promising and validates our model.

Fig. 8

Effect of SA on calibration curves from experimental DH-PSFs acquired from quantum dots (QDs) in glycerol at different depths: (a) 13.5 μm; and (b) 0 μm below the coverslip. (c) Calibration curves from 4 QDs selected in each case [shown by boxes in (a) and (b)] are compared to simulated DH-PSFs with SA at a depth of 13.5 μm and without SA. Lens: 100× / 1.25 NA, oil-immersion. Scale bar: 6.4 μm.

A second study was conducted to match a theoretically determined DH-PSF to a PSF measured from quantum dots (described in Sec. 3.3) using simulations with different amounts of SA. Calibration curves and images from both the experimental and several simulated DH-PSFs are shown in Fig. 9. The asymmetrical range of rotation on either side of the “best focus” (indicating the presence of SA) covered by the experimentally determined calibration curve is determined to be $\sim 4.0\mu \mathrm{m}$ from 90 to $-80\mathrm{deg}$ [Fig. 9(a)]. A gradual increment of the amount of SA in the [Fig. 9(c) to 9(e)] computed for the comparison study, resulted in a reasonable match between the calibration curves and the images of the measured DH-PSF [Fig. 9(b)] and a computed DH-PSF [Fig. 9(c)]. This result is quantitatively supported by Fig. 9(a), which shows that the best agreement in the calibration curves is achieved with an aberrant DH-PSF whose SA is equivalent to the amount induced by a point source located in water at a 100-μm depth below the coverslip [Fig. 9(c)]. The asymmetrical range of the theoretical DH-PSF rotation determined from the 100-μm calibration curve is $\sim 3.9\mu \mathrm{m}$ from 90 to $-90\mathrm{deg}$ which is close to the one determined for the experimental case.

Fig. 9

Comparison of experimental and theoretical DH-PSFs: (a) calibration curves; and (b) to (e) $XY$ cut-view images from experimental (b), and theoretically determined DH-PSFs (c) to (e) of a 91.2× / 1.3 NA oil-immersion lens. Theoretical DH-PSFs determined for a light point-source in water at a depth: (c) $d=100\mu \mathrm{m}$; (d) $d=30\mu \mathrm{m}$; and (e) $d=0\mu \mathrm{m}$. Scale bar: 2 μm.

Due to the noise evident in the experimental data [Fig. 9(b)], the calibration curve in Fig. 9(a) does not have the precision of the experimental calibration curves in Fig. 8(b). As expected, the large amount of SA in the system characterized in Fig. 9 (which is approximately 8 times larger than the amount of SA reported in Fig. 8) causes a reduction in the signal to noise ratio.

4.4.

Effect of SA on the Frequency Content of the DH-PSF

The Fourier transform of the DH-PSF was computed to investigate its frequency content. Modulation transfer functions (MTFs) (i.e., magnitude images of computed 2-D optical transfer functions) are shown in Fig. 10 for conventional and DH-based imaging systems without SA [Fig. 10(a)] and with SA [Fig. 10(b)], respectively. The cutoff frequency for the conventional system without SA in Fig. 10 is $f_{x,\max }=2\mathrm{NA}/\lambda =3791.4$ [$1/\mathrm{mm}$].

Fig. 10

Comparison of 2-D MTFs for a $63\times /1.4\mathrm{NA}$, oil-immersion lens with a 633-nm wavelength computed assuming that the light point source is in water at: (a) a 0-μm depth using a CA system (top row) and a DH system (bottom row); (b) a 2.5-μm depth (resulting in SA) using a CA system (top row) and a DH system (bottom row). Profiles along the horizontal axis through the center of each MTF for: (c) the CA system without SA; (d) the DH system without SA; (e) the CA system with SA; and (f) the DH system with SA.

As evident from this figure, the 2-D DH-MTF shows a rotation as a function of defocus (consistent with the rotation that encodes depth information in the DH-PSF) and a redistribution of frequency attenuation compared to the conventional MTFs. Except for the rotation, there is less change with defocus in the 2-D DH-MTF than in the conventional MTF. This is also evident in the profiles from the corresponding MTFs [Fig. 10(c) to 10(f)] that compare the defocused 2-D MTFs of the two systems. The oscillatory frequency attenuation of the DH-MTF [Fig. 10(d) and 10(f)] shows that the DH-MTF becomes almost zero several times along the spatial frequencies axes, $f_x$ and $f_y$, respectively. SA attenuates the MTF frequencies as expected, but for small defocus the attenuation is less for the DH-MTF than for the MTF of the conventional system.