## Abstract

In this paper, wavefront-encoded (WFE) computational optical sectioning microscopy (COSM) using a fabricated square cubic (SQUBIC) phase mask, designed to render the system less sensitive to depth-induced aberration, is investigated. The WFE-COSM system is characterized by a point spread function (PSF) that does not vary as rapidly with imaging depth compared to the conventional system. Thus, in WFE-COSM, image restoration from large volumes can be achieved using computationally efficient space-invariant (SI) algorithms, thereby avoiding the use of depth-variant algorithms. The fabricated SQUBIC phase mask was first evaluated and found to have a 75% fidelity compared to the theoretical design; it was then integrated in a commercial wide-field microscope to implement a WFE-COSM system. Evaluation of the WFE-COSM system is demonstrated with comparative studies of theoretical and experimental PSFs and simulated and measured images of spherical shells located at different depths in a test sample. These comparisons show that PSF and imaging models capture major trends in experimental data with a 99% correlation between forward image intensity distribution in experimental and simulated images of spherical shells. Our experimental SI restoration results demonstrate that the WFE-COSM system achieves more than a twofold performance improvement over the conventional system of up to a 65 μm depth below the coverslip investigated in this study.

© 2017 Optical Society of America

## 1. INTRODUCTION

Three-dimensional (3D) microscopic techniques enable investigation of biological processes and generation of new knowledge. 3D imaging methods suitable for live specimen investigation should be non-invasive and should provide fast data acquisition. Although these two features are achieved in computational optical sectioning microscopy (COSM), which facilitates high-resolution 3D imaging with a traditional wide-field fluorescence microscope [1–3], challenges due to the depth variability of the point spread function (PSF) [4] have resulted in the use of computationally intensive 3D depth-variant (DV) image reconstruction algorithms for quantitative imaging [5–9]. Because DV algorithms use multiple DV PSFs defined over the imaging volume, they require higher computational resources than space-invariant (SI), a.k.a. deconvolution, algorithms [8].

An alternative approach to overcome the challenges associated with DV imaging systems is based on the idea of engineering the PSF to be less sensitive to depth-induced spherical aberration (SA) while still suitable for COSM [10,11], so that computationally efficient SI methods can be used to reconstruct images over a large volume. During the last 5 years we have been investigating the integration of wavefront coding (traditionally used for extended depth of field microscopy [12]) to COSM [11,13–15], a system we call wavefront-encoded COSM (WFE-COSM). The WFE-COSM PSF, based on a squared cubic (SQUBIC) phase mask (PM), has been shown to have reduced variability over different amounts of SA compared to the PSF of the conventional COSM system [13]. Recently, we presented a promising experimental proof-of-concept study, in which a WFE-COSM system was built using a reflective liquid crystal spatial light modulator (LC-SLM) for the implementation of the SQUBIC PM [15].

Although the LC-SLM has provided us with flexibility in testing multiple designs of the PM [11,14] without changing any hardware in the system, the SLM-based WFE implementation suffers from limitations, such as (i) reduction of light because of the polarizers used in the system; (ii) distortion of the pupil due to oblique illumination; and (iii) finite pixel size of the LC-SLM, which restricts the implementation of a PM with high spatial frequencies [14,15]. For this reason, our previous study was limited to a SQUBIC PM with low spatial frequencies (or lower design parameter, e.g., $A=20$), although it is predicted by theory that a higher design parameter makes the imaging system less sensitive to depth-induced SA over a longer imaging depth [13,15]. In this study, a fabricated SQUBIC ($A=50$) PM is used to implement the WFE system. As the fabricated PM is a transmissive optical element, our WFE-COSM system based on this PM does not suffer from the LC-SLM specific limitations. In addition, the use of a customized fabricated PM provides a compact and cost-effective implementation of a permanent and more accessible system for widespread use. Thus, the current study’s contribution is to demonstrate experimentally the applicability of this type of WFE-COSM implementation using a commercial microscope.

In this paper, a fabricated SQUBIC ($A=50$) PM is investigated in order to evaluate (1) the accuracy of the fabrication process for the implementation of our application-specific PM, and (2) the performance of the WFE system based on this fabricated PM. Initial findings from some of the investigation studies in this paper were presented in a conference publication [16]. The fabricated PM is first evaluated by comparing the wrapped phase measured from an interferometric setup to the theoretical phase. Performance of the WFE-COSM imaging system is evaluated by comparing the experimentally acquired WFE PSF with the simulated PSF. The depth invariance of the WFE imaging system is demonstrated by comparing the experimentally acquired images of test objects at different depths. Finally, the engineered system is validated by reconstructing the experimentally acquired images from different depths using the SI expectation maximization algorithm [17], and compared to the images restored from the conventional imaging system. Similar results obtained from the LC-SLM-based implementation of the WFE imaging are also shown for comparison.

The paper is organized as follows: Section 2 briefly describes the SQUBIC PM and models for the WFE PSF and intermediate image intensity; Section 3 describes methods used in this study, which include the PM fabrication process, experimental setup, and image acquisition. Simulation methods and image reconstruction methods are also described in Section 3. Results from this investigation study are shown in Section 4 and finally, the main findings of this work are summarized in Section 5.

## 2. BACKGROUND

This section briefly describes the SQUBIC PM, the WFE PSF model, and the forward imaging model for the WFE-COSM system.

The phase distribution of the SQUBIC PM defined inside the pupil of the imaging system is given by [13,15]

The WFE PSF is defined by the Fourier optics formulation as

where $u$ and $v$ are the spatial frequency coordinates of the spatial coordinates $x$ and $y$, respectively; $H(u,v;z)$ is the 2D Fourier transform of each axial plane of the complex-valued PSF for COSM with a clear pupil; and ${e}^{j\varphi (u,v)}$ is the phase term due to the PM of Eq. (1).Because depth variability in COSM is minimized by the use of the SQUBIC PM in the WFE system (see, for example, Fig. 5 in [15]), the ideal intensity in the intermediate forward image of the WFE-COSM system, $g(\mathbf{x})$, can be approximated by the convolution between the object intensity, $f(\mathbf{x})$, and the WFE PSF:

where $\mathbf{x}=(x,y,z)$. Because experimental intermediate forward images are perturbed by noise (mainly photon noise), the intensity in the WFE system forward image can be better modeled as the realization from a Poisson random process [18]: where $b$ is a constant that accounts for background emission from the sample, and for CCD camera noise (due to background radiation, luminescent radiation from the chip itself, bias or fat zero, and thermal noise [19]). CCD camera noise due to the on-chip amplifiers, which is Gaussian distributed [19], is not accounted by Eq. (4).## 3. METHODS

This section describes the PM fabrication process, the experimental setup and image acquisition process, as well as image simulation and reconstruction methods used in our investigation studies.

#### A. Phase Mask Fabrication

A custom-made SQUIBIC ($A=50$) PM design was fabricated by Fresnel Technologies Inc. (Fort Worth, TX) on a 24.4 mm diameter and 3 mm thick cynic olefin copolymer plate (Zeonex E48R), which has a RI equal to $1.5355\pm 0.005$ at 515 nm wavelength (Fig. 1). The fabricated PM was mounted in a lens holder [Fig. 1(a)] to facilitate its use in the experimental setup. To obtain a high-quality PM, a broadband anti-reflection coating was added to the surface of the mask. The mask was specifically fabricated for a $63\times /1.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{NA}$ oil immersion objective lens on the center of the glass plate [Fig. 1(a)], which has a diameter equal to the pupil of the objective lens (3.39 mm). The sag profile of the mask [Fig. 1(b)] computed from Eq. (1) represents the unwrapped phase, which can be converted directly to the physical height, $t$, of the mask [Fig. 1(c)] via the known relation $\varphi =2\pi (n-1)t/\lambda $, where $\lambda $ is the design wavelength (equal to 515 nm) and $n$ is the RI of the material used (equal to 1.5355). The height profile was fit to a polynomial equation in order to be able to use the data in the fabrication process. The mask was fabricated using the single-point diamond turning process, which is a well-known method of manufacturing geometrical shapes on polymer plastic [20]. The sag profile machined on the surface was measured by a Tylor Hobson profilometer, which has an accuracy of $\pm 0.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$. The physical dimensions of the PM are shown by the schematics displayed in Figs. 1(d) and 1(e).

#### B. Phase Mask Evaluation

To test the accuracy of the fabricated PM, we measured the wrapped phase produced by a modified Mach–Zehnder interferometer, shown in Fig. 2(a). In the experimental interferometer system, the collimated output of a 4.5 mW laser diode module at center wavelength 532 nm (Thorlabs, Inc., Newton, New Jersey) illuminated the PM. The collimated light emerging from the laser was focused to a point source by a $10\times /0.25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{NA}$ objective lens and then the spherical wavefront was collimated by L1 lens (a 100 mm focal length achromatic doublet). Data were recorded using a CCD camera (Zeiss Axiocam MRm) for the reference wave, the PM wave, and the interference pattern between both imaging paths.

The beam intensity for the reference arm, ${I}_{\text{ref}}(x,y)$, and for the PM arm, ${I}_{\mathrm{PM}}(x,y)$, was captured by blocking the alternate beam arm, and the interfering distribution of both arms, $I(x,y)$, was captured with both paths open. Since any phase measurement is strongly dependent on ambient and source fluctuations, a set of 100 images was captured and averaged. The respective background intensities (captured when the PM was not in the path and when both imaging paths are blocked) were subtracted from ${I}_{\text{ref}}(x,y)$, ${I}_{\mathrm{PM}}(x,y)$, and $I(x,y)$. The cosine of the phase between the reference path and the PM path was computed by the conventional formula [21]:

#### C. Experimental Setup and Image Acquisition

The WFE imaging system was implemented on the side imaging path of a commercial microscope using a 4F optical setup, as shown in Fig. 2(b). The back pupil plane of the objective lens is conjugated to the Fourier plane of the 4F configuration where the fabricated PM is placed. As mentioned earlier, this implementation is an improvement over our previous LC-SLM implementation of the PM [14] because the transmissive fabricated PM allows straight propagation throughout the imaging path without the need for polarizers.

To evaluate the engineered imaging system, the PSF of the SQUBIC WFE system [Fig. 2(b)] was determined by imaging a sample with 170 nm in diameter spherical beads (Invitrogen Molecular Probes, Eugene, Oregon, USA) dried on a coverslip and mounted in optical cement ($\mathrm{RI}=1.56$). The intensity of the measured 3D PSF was corrected by first subtracting a background image (acquired without the fluorescent bead) and then by filtering high-frequency noise using a filter with a cutoff frequency greater than the bandwidth of the optical transfer function.

To assess the advantage of the WFE system [Fig. 2(b)] over the conventional one, images from a test sample were acquired using a $63\times /1.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{NA}$ oil immersion objective lens at 515 nm emission wavelength. The sample was prepared by mounting 6 μm in diameter spherical shells with shell thickness equal to 1 μm (Invitrogen, Molecular Probes, FocalCheck microspheres, 6 μm fluorescent green ring stain/blue throughout) in ProLong Diamond ($\mathrm{RI}=1.47$). The sample was cured at room temperature for 24 h in order to fix the spherical shells inside the mounting medium at different random locations. Images were captured from four different depths below the coverslip at 3, 15, 30, and 65 μm where the centers of spherical shells were determined to be located. To determine the location of the spherical shells, 170 nm in diameter spheres dried on the coverslip were used as a reference of the coverslip-mounting medium interface, which we set as the 0 μm depth. Images from the same sample were also captured using the conventional imaging system, and the SQUBIC ($A=20$) and SQUBIC ($A=50$) WFE imaging systems implemented using a LC-SLM [14] for a comparative study.

All images were captured using a $63\times /1.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{NA}$ oil immersion objective lens for which the design parameter of the fabricated SQUBIC PM ($A=50$) was selected. The CCD camera used in this experiment (Zeiss Axiocam MRm) has a $6.45\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}\times 6.45\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$ pixel, which produces a $0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}\times 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$ pixel in the image space. For capturing a 3D volume, the object was axially scanned at 0.1 μm intervals, thereby providing a cubic voxel dimension ($0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}\times 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}\times 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$) in the acquired image.

#### D. PSF and Forward Image Computation

In this study a $63\times /1.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{NA}$ oil immersion objective lens (similar to the experiment) was simulated. PSFs for this lens were computed using the Gibson and Lanni optical path difference model [4] and the vectorial field approximation theory [22] implemented in MATLAB (Mathworks, Inc., Massachusetts, USA). The PSFs were computed on a $512\times 512\times 1000$ grid, and then cropped to a $400\times 400\times 800$ grid in order to avoid edge wrap around in the fast Fourier transform (FFT) operation required for the WFE PSF computation [Eq. (2)] to use in restoration. Both the axial and lateral spacing of the PSFs were set equal to 0.1 μm, and the emission wavelength used was 515 nm to match the experimental conditions. For the quantitative comparison of experimental and simulated PSFs, PSF values were first scaled to sum up to the same integrated intensity and then normalized so that the maximum value is equal to 1.

Images of 6 μm in diameter spherical shells (with shell thickness equal to 1 μm) located at various depths from 0 μm (just below the coverslip) to 100 μm, at an interval of 10 μm, were simulated to mimic the images acquired from the experimental test sample. WFE system images were simulated by computing Eq. (3) over a $400\times 400\times 1000$ grid using PSFs at respective depths. The amount of SA aberration induced at these depths due to the RI mismatch is within the range of 0 to $-3.8$ (which corresponds to $-1.96\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}/\lambda $) computed using the ${W}_{40}$ term in Eq. (4) of [10].

Noisy images were generated using Eq. (4), where $b$ was determined empirically by averaging the peak to background intensity ratio computed in five measured images, and found to be equal to 1.1% of the peak image intensity. In order to match the noise level in the experimental image, an appropriate scale factor was introduced in the simulated images to apply Poisson noise using MATLAB’s *imnoise* built-in function. The signal-to-noise ratio (SNR) of the noisy simulated images was found to be equal to 25.5 dB, calculated by comparing the energy of the noiseless image to the energy of the noise as follows:

*psnr*.

#### E. Simulated and Experimental Image Restoration

SI image restoration of noisy-simulated and experimental images was used to evaluate performance stability of the WFE-COSM imaging system in the presence of different amounts of SA. Stability was quantified using the correlation coefficient, which is a measure of structure similarity between two images [23]. Iterative and non-iterative regularized algorithms implemented in the *CosmEstimation* module of our open-source software package COSMOS [24] were used in this study.

The non-iterative fast, regularized linear least-square (RLLS) algorithm was used in a study that investigates the performance of the WFE system in the presence of PSF mismatch (due to the imaging depth and the PSF depth used for the SI restoration). Eleven noisy images of spherical shells simulated over a depth range of 0–100 μm at an interval of 10 μm were reconstructed using a single PSF computed at a 50 μm depth below the coverslip. An appropriate value of the regularization parameter ($R$) used in the RLLS algorithm was obtained empirically by computing the similarity (quantified by the correlation coefficient) between the true object and the restored image for different values of $R$ [8]. For this study $R$ was found to be equal to $8\times {10}^{-6}$ for images from the conventional system; $1\times {10}^{-5}$ for images from the SQUBIC ($A=20$) WFE system; and $3\times {10}^{-6}$ for images from the SQUBIC ($A=50$ and 80) WFE system. Because the SQUBIC PM reduces the effective cutoff frequency of the optical transfer function for the WFE system, the impact of regularization becomes more significant especially as the design parameter $A$ increases. Thus, a different value of $R$ was determined for each WFE system.

Experimental images were processed using the iterative SI expectation maximization (SIEM) algorithm instead of the RLLS algorithm, because it has the capability of restoring more details of the unknown object by frequency band extrapolation [25]. To regularize the SIEM results, the Good’s roughness penalty [26] was used with a regularization parameter equal to ${10}^{-5}$, found using the methodology described above. Results after 1000 iterations are displayed and analyzed.

## 4. RESULTS

#### A. Evaluation of Fabricated Mask

The capability of the fabricated PM in producing the designed SQUBIC ($A=50$) PM was evaluated by interferometric measurement, and compared with the LC-SLM-based PM implementation previously reported {Fig. 3(c) of [14]}. The experimental setup and the method of computing the phase are described in Section 3.A. Figure 3 summarizes a comparison of the theoretical PM [Fig. 3(a)] with the wrapped phase computed from the measured experimental data from the fabricated PM [Fig. 3(b)] and the SLM-based PM [Fig. 3(c)]. The correlation coefficient computed between the phase profiles [Fig. 3(c)] of the theoretical and the fabricated PM is 0.75, while the one computed between the SLM-implemented PM and the theoretical phase is 0.12. This demonstrates the improved fidelity achieved by the fabrication process over the SLM implementation of the PM.

Phase value profiles along the horizontal axis through the center of the PMs are also compared [Fig. 3(d)] to quantify the differences between the three PMs. For better visualization of these profile differences, two different regions are magnified and shown in Figs. 3(e) and 3(f). It is seen in Fig. 3(e) that the phase wraps in the case of the SLM PM implementation occur ahead of those in the theoretical phase, whereas the phase wraps of the fabricated PM are consistent with the ones of the theoretical phase (shown by the arrows). This is due to the evident elliptical shape of the PM [Fig. 3(c)], which is a direct result of the oblique illumination used in the SLM implementation [14]. The ratio between the major axis and the minor axis of the retrieved phase in Fig. 3(c) yields a $\mathrm{cos}(\theta )=0.9397$, which corresponds to the angle $\theta =20\xb0$ between the propagation axis and the surface normal of the SLM used in the SLM implementation setup [14]. On the other hand, the fabricated SQUBIC PM preserves the circular shape of the phase [Fig. 3(b)] as it is a transmissive element and the light propagates through it along the axis of propagation. Note that, to account for the effect of the oblique illumination, the pupil aperture in the WFE system was imaged and was found to have 238 pixels along the horizontal direction and 222 pixels along the vertical direction. The phase mask values were fitted inside the image of the pupil aperture, which resulted in an elliptical PM. Notwithstanding, the obliquity of the SLM-based PM is not the main drawback of its performance. The major limitation of the SLM implementation is that high-frequency PMs cannot be sampled properly due to the finite pixel size of the SLM.

An important feature of the SQUBIC PM design is that it has several phase wraps at its edge. In the case of $A=50$, four phase wraps are evident in the range of [1.4–1.7 mm] of the theoretical and fabricated PM, marked by the arrows in Fig. 3(f). However, in the SLM-based PM only three phase wraps are evident because of the LC-SLM’s limitation to implement high spatial frequencies (see [14] for further details). In summary, these results show that the fabricated PM approximates the desired phase function more accurately than the PM implemented with a LC-SLM.

As mentioned earlier, the fabricated SQUBIC PM used in this study was fabricated for a wavelength of 515 nm. Different emission wavelengths can affect the performance of the phase mask. Because the physical height of the fabricated SQUBIC PM is fixed, variation in the wavelength can be understood as an effective change of the effective $A$-parameter (${A}_{\text{eff}}$) of the PM. This means that for any emission wavelength one can compute ${A}_{\text{eff}}$ so that the phase profile produced by the PM coincides with the phase produced by the equation mentioned in Section 3.A. In particular, for the visible spectrum, 400–700 nm, the $A$-value changes from 65 to 37, respectively, in the case of the designed PM. Then, one can conclude that the longer the emission wavelength, the smaller the design $A$-parameter of the SQUBIC PM and, consequently, the effective strength of the PM is reduced. In one previous study [13], the authors verified that the higher the value of the SQUBIC design parameter, $A$ is, the more invariant the WFE system becomes to spherical aberration, so it seems that the imaging performance will be distorted for a longer wavelength. However, it is important to realize that for a longer wavelength the imaging system presents less depth-induced SA [10]; therefore, the imaging system could produce a similar performance. To verify this statement, we studied how the reconstructed image quality is affected for different wavelengths, in both noisy and noiseless simulations. The results of this investigation study are presented in the next section (Section 4.B).

#### B. Evaluation of WFE Imaging System Depth Invariance in the Presence of SA and Noise

To demonstrate the stability of the proposed WFE-COSM system to depth-induced SA, we show results from a comparative study where images of spherical shells (located at different depths in a test sample) simulated for the WFE system and the conventional system in the presence of SA (${W}_{40}=-3.8$) and noise ($\mathrm{SNR}=25.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$) were restored using SI restoration with a single PSF at a depth not necessarily equal to the shell’s depth. In this study we compare the performance of the WFE-COSM system using SQUBIC PMs with a different design parameter $A=20$, 50, and 80 to show the impact of $A$ on system performance. Noisy images from beads located at different depths were simulated following the methods described in Section 3.C and reconstructed as described in Section 3.D. The restored images are compared using the correlation coefficient, $r$, computed between the XZ meridional sections of the restored image of a spherical shell centered at a depth of 50 μm (which is taken as the reference because restoration was performed using the PSF at a depth of 50 μm) with the corresponding section from images of shells at other depths [Fig. 4(a)]. Because of this choice, $r=1$ at a depth of 50 μm and it monotonically decreases as a function of depth change from this reference point.

It is seen in Fig. 4(a) that the curves are not symmetrical around the reference $r=1$ point; the curves show more variability within the range of 0–50 μm in all system cases with the exception of the WFE system for SQUBIC $A=20$. This is because the theoretical PSF varies more rapidly in shallow (0–30 μm) depths [15]. It is observed that the response of the WFE imaging system is more symmetric than that of the conventional system over the same depth range. Based on the curves of Fig. 4, the conventional imaging system shows a 35% variability in terms of SI image restoration with a depth-mismatched PSF over the 100 μm depth range, whereas the variability in the SQUBIC WFE system with $A=50$ and $A=80$ is 14% and 13%, respectively.

To demonstrate the performance of the SQUBIC PM design at different emission wavelengths, we computed the intermediate forward image of a 6 μm spherical shell located at a depth of 50 μm below the coverslip, and the image was restored using two PSFs, computed at the object depth (at 50 μm) and at a 0 μm depth, respectively, to quantify the maximum depth mismatch effect on the restoration. The correlation coefficient was calculated between those two restored images to quantify the restoration accuracy. Figure 4(b) shows the result of this study. From this figure, one can realize that the image quality remains invariant with the variation of the wavelength for both cases. Also, as it is expected, we can observe that the restoration accuracy is higher for the noiseless case than for the noisy simulation. The standard deviations between the correlation coefficients in the case of noiseless and noisy images are 0.0086 and 0.0124, respectively, which is a variability of 3% in the performance of the imaging system over the entire visible spectrum (from 400 to 700 nm). In addition, for our emission filter, whose bandwidth is from 515 to 565 nm, the noisy simulation shows a variability of less than 1% over its entire bandwidth. These findings suggest that the performance of the fabricated SQUBIC PM does not change significantly with the wavelength. The use of the effective $A$-parameter value in the computation of the SQUBIC WFE PSFs for restoration would prevent artifacts due to model mismatch. In our experiments the test object has a peak emission at 515 nm; hence, for the restoration we used the design $A=50$ to compute the SQUBIC WFE PSFs.

#### C. Evaluation of the Experimental WFE PSF

The SQUBIC WFE imaging system with design parameter $A=50$ (Fig. 5) is evaluated by comparing the experimentally measured PSF with the theoretical PSF. A comparison to the SLM PSF is also shown in order to show the benefit of using the fabricated PM over the SLM-based PM in the implementation of the WFE system. Experimental PSFs were measured by capturing images of 170 nm fluorophores (as described in Section 3.B) from just below the coverslip. The meridional sections (XZ view images) of the theoretical WFE PSF and the experimentally acquired PSFs when the fabricated mask and the SLM are used for the WFE implementation are shown in Figs. 5(a)–5(c), respectively. The normalized intensity distributions of the PSFs are compared by the axial intensity profiles taken through the center of the XZ view images of the PSFs [Fig. 5(d)]. A close examination of the intensity profiles shows that although the first four lobes of the theoretical PSF are captured by both experimental PSFs, the next five lobes of the theoretical PSF are evident only in the WFE PSF from the fabricated PM. This comparison shows that although quantitative differences are evident between theoretical and experimental PSF intensities, the experimental WFE PSF based on the fabricated PM captures more accurately the main trends predicted by theory than the SLM-based PSF; the increase in the correlation between profiles was found to be equal to 4%.

#### D. Evaluation of Experimental WFE Intermediate Images from a Test Sample

The performance of the experimental implementation of the SQUBIC ($A=50$) WFE imaging system is also evaluated in terms of the intermediate images of a 6 μm in diameter spherical shell with its center located at 3 μm below the coverslip (Fig. 6). For this evaluation, XZ sections of the simulated WFE image of this test object [Fig. 6(a)] and experimental images acquired from both the fabricated-PM-based [Fig. 6(b)] and the SLM-based [Fig. 6(c)] implementation of the WFE system are compared. Qualitatively, it is observed that the intensity distribution is similar in both the simulated image [Fig. 6(a)] and the experimental image when the WFE system is implemented using the fabricated PM [Fig. 6(b)], whereas the SLM-implemented WFE system produces artifactual intensities shown by the arrow in Fig. 6(c). Quantitative comparison between the experimental and simulated images is performed using the intensity profiles taken through the center of the XZ section images [Fig. 6(d)]. The profiles show that there is a better match between the simulated image and the experimental image when the fabricated PM is used to implement the WFE system than the LC-SLM implementation. The correlation coefficient values between the intensity distributions in the simulated image and the experimental images are equal to 0.99 and 0.93 in the case of the fabricated PM and the SLM-based WFE implementation, respectively.

The stability to SA of the WFE imaging system implemented using the fabricated mask is evaluated by comparing intermediate images acquired from similar spherical shells centered at depths of 3, 30, and 65 μm. The correlation coefficient value between the image of the shell located at a depth of 3 μm and the images of shells located at depths of 30 and 65 μm is found equal to 0.99 and 0.97, respectively. An intensity profile comparison from these images also demonstrates the high degree of similarity between the WFE images of the shells at different depths [Fig. 6(e)]. Because the XZ view images at depths of 30 and 65 μm are so similar to the image of the test object centered at 3 μm [Fig. 6(b)], they are not displayed here for simplicity. The close similarity between the images, despite the increasing depth of the object, demonstrates the depth invariance property of the SQUBIC ($A=50$) WFE imaging system over a 65 μm depth range below the coverslip.

#### E. Experimental Image Restoration of WFE Images

The performance of the experimental SQUBIC WFE imaging system was also evaluated through restoration. As mentioned in Section 3.E, for the restoration the SIEM algorithm and the ideal unaberrated PSF (computed at a depth of 0 μm) were used (for the case of the conventional system, an aberrated PSF was also used for restoration). Images of spherical shells (located at depths of 3, 15, 30, and 65 μm) captured with four different systems were processed in the same way using the SIEM algorithm. Restoration results from this study are summarized in Fig. 7.

The stability achieved in the restored images of shells located at different depths is quantified by the correlation coefficient ($r$) shown at the bottom of each image. The restored images of the shells located at 3 μm were used as the reference in computing $r$ and thus, $r=1$ in these cases except for the conventional system [Fig. 7(a)]. For the conventional system, the image of the shell located at 3 μm restored using an aberrated PSF [Fig. 7(e)], to compensate the restoration artifact evident in Fig. 7(a), was used as the reference instead. For each system, images were compared to the reference image to quantify how robust the system is with respect to SA.

This study demonstrates that the SQUBIC WFE with the design parameter $A=50$ doubles the range of depth in which the SIEM restoration results are artifact-free compared to those obtained with the $A=20$ WFE system [see Figs. 7(b) and 7(d)]. In fact, in the case of SQUBIC WFE ($A=50$), images can be restored with reasonable accuracy ($r=0.81$) up to a 65 μm depth, whereas for the $A=20$ PM, the WFE system has similar performance ($r=0.89$) up to a 30 μm depth. The restoration accuracy reduces significantly at a depth of 65 μm ($r=0.51$) in the case of $A=20$. It is worth mentioning that the SLM-based implementation of the SQUBIC PM for $A=50$ improves the restoration at a depth of 65 μm ($r=0.67$) as observed in Fig. 7(c) compared to $A=20$ [Fig. 7(b)]; however, the restored image in this case suffers from artifacts due to limitations in the SLM implementation of the designed PM (see Section 4.A) and thus, the SLM-based WFE system does not achieve the same performance as the WFE system based on the fabricated PM for $A=50$.

It is worth mentioning that even though results obtained from the WFE system with $A=50$ show more robustness to object depth than in the case of $A=20$, the restored images achieved with $A=50$ show more distortion than those with $A=20$. The fundamental reason for this lies in the basic principle of the wavefront encoding, which modifies the pupil by inserting the SQUBIC phase mask, which introduces a known aberration. The design parameter $A$ of the SQUBIC phase mask controls the amount of the aberration applied to the WFE system. The higher the value of $A$, the stronger the applied aberration is, which subsequently increases the robustness of the imaging system to depth-induced aberrations. The cost of increasing $A$ is the undesirable effective reduction in the cutoff frequency of the optical transfer function, which results in decreased SNR in the forward imaging process and imposes some distortion in the final reconstructed images. Thus, there is a trade-off between high fidelity restoration and high imaging stability over depth in the SQUBIC WFE system.

Because of a residual aberration at the coverslip in our system [15], in the case of the conventional image restoration, the SIEM with the ideal PSF fails to reconstruct the shape of the object Fig. 7(a). By proper determination of the aberration present in the imaging system, it was possible to obtain a better reconstruction [Figs. 7(e)–7(f)]. However, this process is not practical for an arbitrary biological sample because the aberration needs to be determined and accounted for at every imaging depth via depth-variant image reconstruction.

To show a quantitative comparison of system performance in the presence of SA, normalized axial intensity profiles through the center of restored images of shells located at a 65 μm depth are plotted in Fig. 7(h). The intensity profile from a simulated 6 μm spherical shell [Fig. 7(g)] is also included as ground truth. It is observed that the result obtained with the SQUBIC WFE system using the fabricated ($A=50$) PM is the closest to the ground truth over all the other results.

## 5. SUMMARY AND CONCLUSION

In this study, the implementation of a WFE-COSM system using a fabricated SQUBIC ($A=50$) PM is presented. The fabricated PM allows a compact, cost-effective implementation of the WFE system, free from LC-SLM specific limitations encountered in our previously demonstrated experimental implementation of a WFE-COSM system [15]. Results presented here show that the use of the fabricated PM increases the performance of the WFE-COSM imaging system in terms of achieving artifact-free SI restoration in the presence of SA (up to ${W}_{40}=-3.8$) and noise ($\mathrm{SNR}=25.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$).

The fabricated PM was evaluated through interferometric measurement, and was compared with the LC-SLM implementation of the PM, as well as with the theoretical phase of the SQUBIC PM. Our results (Fig. 3) show 75% fidelity in the implementation of the fabricated PM compared to the theoretical design, which represents a 63% increase in fidelity compared to that achieved with the LC-SLM implementation of the PM.

The utility of the fabricated PM was evaluated by comparing experimentally measured WFE PSFs with the theoretical one (Fig. 5). The high implementation accuracy of the fabricated PM is evident in the agreement of major features of the PSFs. In addition, comparison of intermediate simulated and experimental WFE images [Figs. 6(a)–6(d)] shows a 99% agreement, which further validates the implementation. The insensitivity of the WFE imaging system to SA is demonstrated by the high similarity of experimental images of the same test objects located at depths of 3–65 μm [Fig. 6(e)]; the computed variability between the intensity distributions of these images is only 3%.

Experimental evaluation of the implemented WFE-COSM system through restoration shows several advantages. The performance of the WFE-COSM system based on the SQUBIC ($A=50$) fabricated PM is increased by up to 14% over the SLM-based implementation for the same $A$. The stability of the implemented imaging system for $A=50$ is demonstrated over a 65 μm depth below the coverslip, which is more than double the depth of 30 μm achieved with a SQUBIC ($A=20$) PM in a previous implementation. This promising result encourages the fabrication of a new PM with a higher $A$-value. A study to determine the optimal value of the SQUBIC design parameter $A$ for WFE-COSM will be reported in a future publication.

## Funding

National Science Foundation (NSF) (0844682, 0852847, 1353904); University of Memphis (UoM); Herff Graduate Fellowship, Dept. of EECE, University of Memphis; Ministerio de Economia y Competitividad (DPI2015-66458-C2-1-R); Generalitat Valenciana (PROMETEOII/2014/072).

## Acknowledgment

The authors are thankful to S. V. King for help with the experimental setup and for acquiring the SLM-based WFE PSF; M. Martínez-Corral (Department of Optics, University of Valencia, Spain) for useful technical discussion; and L. H. Schaefer (Advanced Imaging Methodology Consultation, Ontario, Canada) for PSF MATLAB code.

## REFERENCES

**1. **D. A. Agard, “Optical sectioning microscopy: cellular architecture in three dimensions,” Annu. Rev. Biophys. Bioeng. **13**, 191–219 (1984). [CrossRef]

**2. **D. A. Agard, Y. Hiraoka, P. Sha, and J. W. Sedat, “Fluorescence microscopy in three dimensions,” Methods Cell Biol. **30**, 353–377 (1989). [CrossRef]

**3. **J. R. Swedlow, J. W. Sedat, and D. A. Agard, “Deconvolution in optical microscopy,” in *Deconvolution of Images and Spectra* (Academic, 1997), Chap. 9.

**4. **S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A **9**, 154–166 (1992). [CrossRef]

**5. **C. Preza and J.-A. Conchello, “Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A **21**, 1593–1601 (2004). [CrossRef]

**6. **M. Arigovindan, J. Shaevitz, J. McGowan, J. W. Sedat, and D. A. Agard, “A parallel product-convolution approach for representing the depth varying point spread functions in 3D widefield microscopy based on principal component analysis,” Opt. Express **18**, 6461–6476 (2010). [CrossRef]

**7. **S. Ben Hadj, L. Blanc-Féraud, G. Aubert, and G. Engler, “Blind restoration of confocal microscopy images in presence of a depth-variant blur and Poisson noise,” in *Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing* (IEEE, 2013), pp. 915–919.

**8. **N. Patwary and C. Preza, “Image restoration for three-dimensional fluorescence microscopy using an orthonormal basis for efficient representation of depth-variant point-spread functions,” Biomed. Opt. Express **6**, 3826–3841 (2015). [CrossRef]

**9. **S. Ghosh and C. Preza, “Three-dimensional block-based restoration integrated with wide-field fluorescence microscopy for the investigation of thick specimens with spatially variant refractive index,” J. Biomed. Opt. **21**, 046010 (2016). [CrossRef]

**10. **G. Saavedra, I. Escobar, R. Martínez-Cuenca, E. Sánchez-Ortiga, and M. Martínez-Corra, “Reduction of spherical-aberration impact in microscopy by wavefront coding,” Opt. Express **17**, 13810–13818 (2009). [CrossRef]

**11. **S. Yuan and C. Preza, “Point-spread function engineering to reduce the impact of spherical aberration on 3D computational fluorescence microscopy imaging,” Opt. Express **19**, 23298–23314 (2011). [CrossRef]

**12. **E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. **34**, 1859–1866 (1995). [CrossRef]

**13. **A. Doblas, S. V. King, N. Patwary, G. Saavedra, M. Martínez-Corral, and C. Preza, “Investigation of the SQUBIC phase mask design for depth-invariant widefield microscopy point-spread function engineering,” Proc. SPIE **8949**, 894914 (2014). [CrossRef]

**14. **S. V. King, A. Doblas, N. Patwary, G. Saavedra, M. Martínez-Corral, and C. Preza, “Spatial light modulator phase mask implementation of wavefront encoded 3D computational-optical microscopy,” Appl. Opt. **54**, 8587–8595 (2015). [CrossRef]

**15. **N. Patwary, S. V. King, G. Saavedra, and C. Preza, “Reducing effects of aberration in 3D fluorescence imaging using wavefront coding with a radially symmetric phase mask,” Opt. Express **24**, 12905–12921 (2016). [CrossRef]

**16. **N. Patwary, S. V. King, H. Shabani, and C. Preza, “Experimental implementation of wavefront encoding in 3D widefield fluorescence microscopy using a fabricated phase mask designed to reduce system depth variability,” in *Classical Optics 2014*, OSA Technical Digest (online) (2016), paper CW2D-3.

**17. **J.-A. Conchello and J. G. McNally, “Fast regularization technique for expectation maximization algorithm for optical sectioning microscopy,” Proc. SPIE **2655**, 199–208 (1996). [CrossRef]

**18. **M. Bertero, P. Boccacci, G. Desiderà, and G. Vicidomini, “Image deblurring with Poisson data: from cells to galaxies,” Inverse Probl. **25**, 123006 (2009). [CrossRef]

**19. **D. L. Snyder, A. M. Hammoud, and R. L. White, “Image recovery from data acquired with a charge-coupled-device camera,” J. Opt. Soc. Am. A **10**, 1014–1023 (1993). [CrossRef]

**20. **S. Bäumer, *Handbook of Plastic Optics* (Wiley, 2011).

**21. **B. E. A. Saleh and M. C. Teich, *Fundamentals of Photonics* (Wiley, 1991).

**22. **O. Haeberlé, “Focusing of light through a stratified medium: a practical approach for computing microscope point spread functions. Part I: conventional microscopy,” Opt. Commun. **216**, 55–63 (2003). [CrossRef]

**23. **Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. **13**, 600–612 (2004). [CrossRef]

**24. **Computational Imaging Research Laboratory, Computational Optical Sectioning Microscopy Open Source (COSMOS) software package; http://cirl.memphis.edu/COSMOS.

**25. **J.-A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A **15**, 2609–2619 (1998). [CrossRef]

**26. **I. J. Good, “Non-parametric roughness penalty for probability densities,” Nat. Phys. Sci. **229**, 29–30 (1971). [CrossRef]