Reducing effects of aberration in 3 D fluorescence imaging using wavefront coding with a radially symmetric phase mask

In this work, a wavefront encoded (WFE) imaging system built using a squared cubic phase mask, designed to reduce the sensitivity of the imaging system to spherical aberration, is investigated. The proposed system allows the use of a space-invariant image restoration algorithm, which uses a single PSF, to restore intensity distribution in images suffering aberration, such as sample–induced aberration in thick tissue. This provides a computational advantage over depth-variant image restoration algorithms developed previously to address this aberration. Simulated PSFs of the proposed system are shown to change up to 25% compared to the 0 μm depth PSF (quantified by the structural similarity index) over a 100 μm depth range, while the conventional system PSFs change up to 84%. Results from experimental test-sample images show that restoration error is reduced by 29% when the proposed WFE system is used instead of the conventional system over a 30 μm depth range. ©2016 Optical Society of America OCIS codes: (110.0110) Imaging systems; (110.0180) Microscopy; (110.7348) Wavefront encoding; (070.6110) Spatial filtering; (100.3190) Inverse problems; (110.1758) Computational imaging. References and Links 1. J. R. Swedlow, J. W. Sedat, and D. A. Agard, Deconvolution of Images and Spectra (Academic Press, 1997), Chap. 9. 2. C. Preza and J.-A. Conchello, “Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A 21(9), 1593–1601 (2004). 3. 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(10), 3826–3841 (2015). 4. 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(10), 2609–2619 (1998). 5. 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(7), 6461–6476 (2010). 6. A. Wong, X. Y. Wang, and M. Gorbet, “Bayesian-based deconvolution fluorescence microscopy using dynamically updated nonstationary expectation estimates,” Sci. Rep. 5, 10849 (2015). 7. E. A. Mukamel, H. Babcock, and X. Zhuang, “Statistical deconvolution for superresolution fluorescence microscopy,” Biophys. J. 102(10), 2391–2400 (2012). 8. 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(4), 046010 (2016). 9. B. Kim and T. Naemura, “Blind Depth-variant Deconvolution of 3D Data in Wide-field Fluorescence Microscopy,” Sci. Rep. 5, 9894 (2015). 10. E. Maalouf, B. Colicchio, and A. Dieterlen, “Fluorescence microscopy three-dimensional depth variant point spread function interpolation using Zernike moments,” J. Opt. Soc. Am. A 28(9), 1864–1870 (2011). 11. L. Silvestri, L. Sacconi, and F. S. Pavone, “Correcting spherical aberrations in confocal light sheet microscopy: A theoretical study,” Microsc. Res. Tech. 77(7), 483–491 (2014). #263064 Received 13 Apr 2016; revised 23 May 2016; accepted 25 May 2016; published 3 Jun 2016 © 2016 OSA 13 Jun 2016 | Vol. 24, No. 12 | DOI:10.1364/OE.24.012905 | OPTICS EXPRESS 12905 12. A. Masson, P. Escande, C. Frongia, G. Clouvel, B. Ducommun, and C. Lorenzo, “High-resolution in-depth imaging of optically cleared thick samples using an adaptive SPIM,” Sci. Rep. 5, 16898 (2015). 13. 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(23), 23298–23314 (2011). 14. G. Saavedra, I. Escobar, R. Martínez-Cuenca, E. Sánchez-Ortiga, and M. Martínez-Corral, “Reduction of spherical-aberration impact in microscopy by wavefront coding,” Opt. Express 17(16), 13810–13818 (2009). 15. 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(29), 8587–8595 (2015). 16. S. Tucker, W. T. Cathey, and E. Dowski, Jr., “Extended depth of field and aberration control for inexpensive digital microscope systems,” Opt. Express 4(11), 467–474 (1999). 17. 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). 18. 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(1), 154–166 (1992). 19. N. Patwary, S. V. King, and C. Preza, “3D microscope imaging robust to restoration artifacts introduced by optically thick specimens,” Proc. SPIE 9330, 93300O (2015). 20. N. Patwary, A. Doblas, S. V. King, and C. Preza, “Reducing depth induced spherical aberration in 3D widefield fluorescence microscopy by wavefront coding using the SQUBIC phase mask,” Proc. SPIE 8949, 894911 (2014). 21. N. Patwary and C. Preza, “Wavefront encoded computational optical sectioning microscopy reduces depth variability in 3D imaging,” in Imaging and Applied Optics 2015, OSA Technical Digest (online) (Optical Society of America, 2015), paper CM2E.4. 22. M. Persson, D. Engström, and M. Goksör, “Reducing the effect of pixel crosstalk in phase only spatial light modulators,” Opt. Express 20(20), 22334–22343 (2012). 23. J.-A. Conchello and E. W. Hansen, “Enhanced 3-D reconstruction from confocal scanning microscope images. 1: Deterministic and maximum likelihood reconstructions,” Appl. Opt. 29(26), 3795–3804 (1990). 24. J.-A. Conchello and J. G. McNally, “Fast regularization technique for expectation maximization algorithm for optical sectioning microscopy,” Proc. SPIE 2655, 199–208 (1996). 25. I. J. Good, “Non-Parametric Roughness Penalty for Probability Densities,” Nature 229(1), 29–30 (1971). 26. 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(1), 55–63 (2003). 27. Computational Imaging Research Laboratory, “Computational Optical Sectioning Microscopy Open Source (COSMOS) software package,” http://cirl.memphis.edu/COSMOS]. 28. 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(4), 600–612 (2004).


Introduction
Wide-field fluorescence microscopy is preferable in live-cell imaging because it allows low light for excitation and fast data acquisition thereby protecting the sample from photobleaching.Although three-dimensional (3D) wide field images are corrupted by out-of-focus light and spherical aberration (SA) due to the refraction of light at the interface of inhomogeneous media, these intermediate images can be significantly improved using computational methods [1][2][3].This approach, known as computational optical-sectioning microscopy (COSM), has facilitated 3D high resolution imaging with a wide-field microscope [2,[4][5][6][7].
When imaging a point light source located deep within the sample, the microscope's point-spread function (PSF) changes with imaging depth (due to depth-induced SA) resulting in a depth-variant (DV) imaging process [2].Attempts to restore resolution with spaceinvariant (SI) methods introduce artifacts.To properly account for this variability with computational methods, image reconstruction algorithms for COSM have been developed based on DV imaging models that use information from multiple DV PSFs [2,3,[8][9][10].However, these DV reconstruction algorithms are costly in terms of both the memory requirement and data processing time compared to the initially-developed SI algorithms, which use only a single PSF [2,3].Reducing the imaging system's depth variability can enable SI restoration of images from thick samples with reduced restoration artifacts for quantitative COSM.
A benefit of the approach presented here is that it provides an alternative method for addressing the effect of SA in applications for which the computational requirements of depth-variant algorithms are not practical for biological investigation, such as live cell imaging.Another potential application of the approach presented here for addressing effects of SA, not explicitly modeled in this study, is in light sheet microscopy of cleared tissue samples [11,12].The use of depth-variant algorithms for restoration is also impractical in this application, because the macroscopic size of the samples in such studies results in large data volumes.Furthermore, a known effect of SA is an axial shift that causes a difference between the apparent versus the actual depth of a point source below the coverslip.By reducing the effect of SA this difference is also reduced, thereby increasing the accuracy of quantitative 3D image analysis in studies of material (e.g.cell and tissue) structure.
Wavefront encoding was recently applied to COSM to reduce the effect of depth induced SA [13][14][15], although it was initially developed for extended depth-of-field (EDF) microscopy [16].Simulation results by Yuan and Preza in [13], and Saavedra et al in [14] introduced the idea of wavefront encoded (WFE) COSM to address depth-induced SA.In this study, WFE system performance evaluation through restoration of experimental and simulated wide field microscope images is presented.The WFE system is implemented using a redefined squared cubic (SQUBIC) phase mask (suitable for experimental implementation), which (SQUBIC) was initially developed to reduce the effect of depth induced SA in 3D confocal microscope, and was demonstrated in a proof-of-concept study (without any experimental verification) using simulated confocal microscopy images [14].The design methodology of the SQUBIC phase mask was described by Doblas et al in [17], and its comparison with other phase masks, such as, the cubic phase mask (CPM) [18] and the generalized cubic mask (GCPM) [13] was described by King et al in [15].The study presented here includes experimental results for imaging model validation for SQUBIC-WFE wide-field microscopy.This is an important step for WFE COSM, which relies on algorithms based on an accurate image formation model that accounts for imaging system attributes in experimental setup.
This paper presents results to show the usability and performance of the WFE system using the SQUBIC phase mask experimentally, which we refer to as SQUBIC-WFE system.The aim is to investigate the ability of WFE-COSM based on the SQUBIC phase mask to reduce the effect of SA on image restoration of 3D wide field microscopy images.The design of the SQUBIC phase mask depends on the microscope objective (MO) lens' numerical aperture (NA) and the refractive index (RI) of the immersion medium [14]; therefore, this study investigates system performance with two types of objective lenses: a 20X/0.8NA dry and a 63X/1.4NA oil immersion, and quantifies performance for different imaging conditions and different amounts of SA.Preliminary results from the study were presented in conference publications [19][20][21].Simulated SQUBIC-WFE images of a 3 µm in diameter solid sphere at different depths within the embedding medium are restored using a single PSF computed at depth 30 µm.The effect of using a single PSF to restore images from a thick sample is investigated using a 40 µm thick object using the space-invariant expectation maximization (SIEM) algorithm.The SQUBIC-WFE system is validated by comparing: 1) experimentally acquired images with simulated intermediated images of test objects results; and 2) restorations from these images.We used a liquid crystal spatial light modulator (LC-SLM) to experimentally implement the SQUBIC phase mask within a commercial microscope.The SQUBIC-WFE restored images were also compared to corresponding conventional restored images.
The paper is organized as follows: Section 2 reviews the theory and mathematical background of the SQUBIC-WFE system, Section 3 describes the methods used in the investigation studies, Section 4 describes the results obtain from different investigation studies, and Section 5 summarizes findings.

Background
This section reviews the mathematical description of the SQUBIC phase mask, the WFE PSF, the intermediate image formation model, and the reconstruction algorithm used in this study.A brief overview of limitations to SLM-based implementation of the SQUBIC phase mask is also given.

Theory and implementation of SQUBIC phase mask
The phase function of the SQUBIC phase mask ( , ) x y f f ϕ is defined as the following equation [20]: where, x f  and y f  are normalized pupil coordinates; ( ) where n is the RI of the lens' immersion medium, NA is the numerical aperture of the MO lens; and A is a design parameter that determines the strength of the phase mask.As evident from Eq. ( 1) the SQUBIC phase mask is radially symmetric and it depends on imaging system parameters in order to directly compensate for the presence of SA in the imaging system.The shape of the SQUBIC PSF varies less with depth-induced aberration when the value of A is increased as predicted by analytical results [17].However, limitations related to the characteristics of the liquid crystal spatial light modulator (LC-SLM), used in this study for the implementation of the phase mask, and the configuration of the WFE imaging path (discussed in detail by King et al. in [15], and Persson et al in [22]) degrade the fidelity of the PSF and the WFE-COSM system, and thus restrict effective experimental implementation of the SQUBIC phase mask to a maximum value of A = 20 when a 63X/1.4NA oil lens is used.In this study, SQUBIC phase masks with different design parameter (A) values and imaging lenses (e.g., 20X/0.8NA dry lens and 63X/1.4NA oil lens) were used.The choice of these values in this study is related to theoretical predictions, the pixel size of the LC-SLM and its effective ability to implement the phase mask.The highest possible value of A (A = 88) determined based on the sampling limit of the LC-SLM [15] was used with the low NA lens.For the case of the high NA lens the value of A was set to 20 in order to be able to compare simulated results with experimental results, which were acquired with the phase mask design for A = 20.Figure 1(a) shows the unwrapped SQUBIC phase function displayed along the pupil plane for the 20X/0.8NA dry objective lens with A = 88 and, Fig. 1(d) shows the corresponding wrapped phase pattern.Note that, phase functions need to be wrapped from zero to 2π when applied to the generalized pupil using an LC-SLM.
It is seen from the profiles [Fig.1(c)] that, for A = 88 (NA = 0.8) the gradient of the phase mask is higher than in the case of the A = 20 (NA = 1.4) mask, which is consistent with the higher special frequency content of its corresponding wrapped phase mask [Fig.1(d)] compared to that of the A = 20 mask [Fig.1(e)].Simulation results for values of A that could not be reproduced experimentally with the current SLM implementation, were also investigated in this study because they provide a benchmark for the expected performance of the mask at these design parameters by future experimental implementations based on alternative configurations [for example through a fabricated phase mask].Please note that, as the experimental limitation of this study to implement phase masks with higher design parameter (A) is related to the finite pixel size of the LC-SLM, this will not be the case for a fabricated phase mask.In the case of a static phase mask fabricated of glass or polymer plastic [i.e.cyclo-olefin polymer] the required optical path difference will be directly proportional to the slowly varying unwrapped phase function [Fig.1(a-c)] rather than a modulo 2π phase distribution [Fig.1(d-e)].In contrast, when the LC-SLM is used, a wrapped phase function [i.e.0 to 2π], which has high spatial frequencies [Fig. 1 (d-e)], is projected directly onto the SLM surface.

SQUBIC WFE-PSF
The WFE PSF is defined in terms of the conventional PSF [18] through a modification of the generalized pupil function with the SQUBIC phase mask, ( , ) x y f f ϕ , as in the Fourier optics formulation [13]: , ( , ) ( , ) , where 1  (3)

Image formation and restoration
For an extended 3D object, the intensity distribution of the image is described by the superposition integral [2]: where ( , , ) is a 3D point in the image space; ( , , ) is a 3D point in the object space O; and ( ) f o x is the object intensity.Equation (4) has successfully approximated by the strata-based forward imaging model, which uses only a finite number of PSFs over the entire thickness of the sample [2].
The inverse imaging problem based on Eq. ( 4) is depth-variant and its solution requires multiple PSFs [2].As the proposed SQUBIC-WFE system is designed to have reduced depth sensitivity over the conventional system, a reasonable approximation is to assume that, the PSF does not change significantly over depth.Consequently, this allows the use of a single PSF with SI restoration algorithms for the solution of the inverse problem.In this study, image restoration was accomplished using the SI expectation-maximization (SIEM) algorithm [23].In the SIEM algorithm, the image intensity at the (k + 1) th iteration is estimated by updating the current estimate, ( , , ) k s x y z , by a back projection to the object space of the measured image, ( , , ) d x y z , divided by the model prediction at the k th iteration, ( , , ) k g x y z calculated using Eq. ( 4), as in the following iteration step [23]: where H is the sum of intensity in the 3D PSF, ( , , ; ) o h x y z z .Because the inverse imaging problem is ill-posed, regularization was introduced to stabilize the solution obtained with the SIEM algorithm.This was achieved by maximizing a penalized log-likelihood function instead [24]: which includes the Good's roughness penalty functional, [ ] x [25], and γ is the regularization parameter that determines the contribution of the penalty to the penalized maximum log-likelihood and controls the amount of smoothness applied to the restored image.The value of γ trades off restoration error due to noise present in the experimental images with loss of resolution in the restored image.In practice, this is a user defined tuning parameter, which needs to be experimentally adjusted for different types of image restoration problems.

Methods
This section describes the methods followed for sample preparation, experimental image acquisition, computation of conventional and SQUBIC-WFE PSFs, generation of simulated intermediate images from numerical objects, and image restoration.Metrics used to quantify performance comparisons are also discussed.

Experimental image acquisition and system calibration
For the data acquisition, a modified ZEISS AxioImager.Z1 was used.This commercial upright microscope has two imaging paths.Conventional imaging was performed through the top imaging path, while WFE imaging was implemented through a 4F imaging system at the side imaging path, which was equipped with a reflective LC-SLM, as described in [15].The Fourier plane of the 4F setup was set conjugate to the back focal plane of the microscope objective lens and it was equipped with an LC-SLM, enabling the application of the implemented SQUBIC phase mask.In this study, the SQUBIC-WFE system was implemented experimentally for the design parameter A = 20 due to limitations of the LC-SLM to adequately capture high spatial frequencies present in a phase mask with a larger value of A [15].The test object used was a set of microspheres, each having a 6 µm in diameter spherical shell with shell thickness equal to 1 µm (Invitrogen, Molecular Probes®, FocalCheck microspheres, 6 µm fluorescent green ring stain/blue throughout).These micro-shells have two fluorescent markers: blue that labels them throughout, and green that labels only the outer shell.A solution of the beads was diluted in deionized water and dried both on a slide and on a ZEISS high precision cover slip (RI = 1.522).Smaller beads (nano-spheres), 170 nm in diameter, were also dried on the coverslip to be used as depth reference.The beads were mounted in ProLong® Diamond Antifade Mountant (Invitrogen, Molecular Probes®) with RI = 1.47 when cured.During the mounting process, the beads on the slide were mixed into the uncured ProLong Diamond® material thereby distributing them to different depths.Images were captured using the 63X/1.4NA oil immersion (RI = 1.518) objective lens from depths d 0 , d 0 + 12 µm, and d 0 + 27 µm where d 0 is approximately 3 µm.The depths were measured from the relative distance between the center of the spherical shells and the nano-spheres' layer (170 nm beads).
Ideally, there should not be any SA just below the coverslip (i.e. at 0 µm depth), for an oil immersion objective lens as there is no RI mismatch [18]; however, in practice this may not always be the case.For example, in the commercial microscope used in this study (Zeiss Axio Images.Z1), some aberration was observed at the ideal imaging condition (i.e. at 0 µm depth), which we refer to as residual aberration (RA).Due to this RA, images of objects just below the coverslip show aberration similar to SA, which depends on the lens used in the experiment.To measure this RA for an oil objective lens, we dried 170 nm beads on the coverslip and mounted them in the same immersion oil used for the lens.The RA is lens dependent and requires careful computation for each microscope objective lens.Images of nano-spheres just below the coverslip, which are essentially the PSFs of the system at the coverslip, were captured using the 63X/1.4NA oil immersion objective lens [see Fig. 2].To calculate the RA, the experimental image was compared with simulated PSFs computed at different depths (d) using different amounts of RI mismatch ( RI Δ ).The structural similarity index measure, defined in Section 3.5, was calculated between the simulated and the experimental images to determine the best match between the experimental and simulated PSF.
Figure 2(a) shows the XZ sectional view of the experimental PSF measured at 0 µm depth and without a RI mismatch between the immersion medium of the MO lens and the object's mounting medium.Therefore, there should be no sample-induced SA, and the PSF should be symmetric along the axial direction [as shown in Fig. 2(b)].The RA introduces a phase distortion at the back focal plane of the MO lens which can be modeled as SA.In this case, the phase distortion was found to be equivalent to that produced by SA in the case of 0.068 RI Δ = and depth, d = 8 µm.The appropriate phase distortion was added to the generalized pupil of the PSF using the Gibson and Lanni optical path difference model and the corresponding d and RI Δ [18].The computed PSF using the estimated phase distortion [Fig.2(c)] shows qualitative agreement with the experimental PSF [Fig.2(a)].The intensity profile comparison shown in Fig. 2(d) supports this observation quantitatively.

Computation of conventional and WFE PSFs
To calculate the SQUBIC WFE-PSFs, complex valued 3D conventional PSFs were computed over a clear circular aperture (CCA) (i.e., no wavefront distortion is introduced at the aperture).3D conventional PSFs were computed using the Gibson-Lanni optical path difference (OPD) model for fluorescence microscopy [18] and vectorial field approximation theory [26] implemented in Matlab (MathWorks, Inc.).The field components were computed using the Legendre-Gauss quadrature numerical integration method with n th order Legendre polynomials where n = 1000.In order to avoid aliasing in the WFE-PSFs, conventional PSFs were computed over a large grid using oversampling.The phase mask was then applied to each axial plane of the PSF in accordance with Eq. ( 2) to obtain the corresponding SQUBIC-PSF plane.The discretized phase mask patterns [Fig.1] were resized to match the pupil aperture size by determining the number of pixels along the mask's diameter, equal to 2 / NA λ , given by: 2( ) , where x Δ is the lateral sampling of the PSF and κ is the number of pixels along the lateral direction.Finally, the phase mask patterns were padded with zeros to match their dimension to the grid size used in the computations.Examples of the conventional and the SQUBIC WFE PSFs are shown in Fig. 3 (a-f).

Simulated object and intermediate image
To demonstrate the performance of the engineered PSF with the SQUBIC phase mask in simulation, three different types of numerical objects were simulated.The first object, referred to as object 1 [Fig.4(a)], has one 3-µm in diameter spherical solid bead.The bead's depth location was varied within a 20-40 µm range with an interval of 2 µm and simulated images were computed for each case in order to investigate accuracy in restoration obtained using a PSF at depth 30 µm.Test object 2 [Fig.5(a)] has five 3-µm in diameter spherical solid beads located at different depths within a 10-50 µm depth range with an interval of 10 µm.Objects 1 and 2 were simulated on a 336 × 336 × 600 grid where each voxel size is 0.1 µm × 0.1 µm × 0.1 µm.Because the SQUBIC-PSF extends significantly in the axial direction [Fig.3(c-f)], the objects were zero padded to a final grid size equal to 336 × 336 × 1250 voxels to completely capture the intensity spread in the blurred images.Both objects were smoothed using a Gaussian kernel to reduce sharp edges (which are not present in real objects) before their simulated images were computed.The forward images were generated using the strata-based model [2] with five strata using six DV PSFs computed in the depth range of 0 -50 µm at a 10 µm depth step) using the COSMOS software package [27] which implements the strata-based forward imaging model.The simulated forward images for the conventional and SQUBIC system are shown in Fig. 5(b & c), respectively.In the case of object 1 and object 2, the design parameter A = 88 was used to implement the SQUBIC phase mask [Eq.( 1)] for a 20x, 0.8NA lens.A is chosen according to the theoretical maximum limit of the SLM due to finite pixel size in order to investigate the maximum capability of the phase mask within theoretical design limit of the SLM used for the implementation [15].The third object used in this study approximates the actual test sample constructed using 6-µm in diameter spherical shells with a shell thickness equal to 1 µm (Sec.3.1).The centers of the objects were simulated sequentially at depths 3 µm, 12 µm and 27 µm.The selection of this particular test object and object depths are chosen to match with the experiments (Sec.3.1) so that the simulated results can be validated.As object 3 is bigger compared to the object 1 and object 2, it was simulated on a bigger grid (512 × 512 × 1500).The forward images [Fig.6] were simulated following the same methodologies described above using three strata (i.e.four DV PSFs in the depth range of 0 µm to 30 µm at a 10 µm depth step).As the sample thickness is smaller in this case (30 µm) compared to the other two objects (50 µm), fewer PSFs were used to simulate the forward image keeping the size of each stratum constant (10 µm) [2].WFE image of these objects were simulated with different values of the SQUBIC design parameters (A = 20, 50, and 80); however, only A = 20 was compared with the experiment [due to limitations in the experimental implementation explained in Section 2.1].The different investigation studies are summarized in Table 1.

Image restoration
The images were restored using the SIEM algorithm [4] described in Section 2.3.This algorithm is implemented in the CosmEstimation module of the COSMOS software package [27].The first type of object from different depth locations (objects 1 &2 were restored with the PSF computed at 30 µm depth.Object 3 (6 µm in diameter spherical shell) was restored using PSFs computed at depths 3 µm, 12 µm and 27 µm.The noiseless simulated images were restored without applying any regularization penalty, whereas the experimental restorations were regularized by the roughness penalty with the regularization parameter γ equal to 0.005.We set the number of iteration to 1000 as in this study we found that the shape of the restored images does not change significantly with subsequent iterations.For a 3D image of size 512x512x1500 the SIEM algorithm took 1.44 minutes/iteration on a Quad Opteron 6274 16C, 2.2 GHz processor.Processing time scales linearly with grid size.

Performance evaluation metrics
To investigate the change of PSFs over depth in the case of conventional and SQUBIC imaging systems, the quantitative comparison was performed by means of the structural similarity index measure (SSIM), which compares the luminescence, contrast, and structure between two 3D images [28].The luminescence and contrast are assessed by measure of the mean and variance, respectively, of the collection of image pixel intensities.From these statistical parameters, the SSIM between two 3D volumetric images X and Y is computed as [28]: where , σ σ are the mean and variance of the collection of pixel values in the corresponding images; XY σ is the covariance between the pixel values in the images; and 1 c and 2 c are two constants that restrict the value of SSIM from becoming infinite.The values of SSIM range from −1 to + 1, however, in the case of fluorescence microcopy, positive intensities are captured, and thus SSIM cannot be negative.Therefore, in this study, the minimum and maximum limits of the SSIM are 0 and 1, which refers to a 0% or 100% match between the intensity distributions of two images, respectively.
To investigate the restoration performance of the conventional and the SQUBIC-WFE systems, the correlation coefficient (ρ) was used to quantify the structural variation between images.The test objects used in this study are spheres or spherical shells which have uniform intensity throughout the solid structure; therefore, there is little luminescence variation and an approximately binary contrast in the test object.For this reason, the correlation coefficient was used to compare XZ sections taken from the 3D images.In simulation, direct comparison between the intensity distribution of the test object and the restored image data is possible.However, for the experimental cases, as the true object is not known, the restored image that appears to be closest to a micro ring (judge by visual inspection) was used as the reference image to compute ρ with other restored images.The correlation coefficient between two 2D images X and Y is given by: where 3 c is a constant chosen to avoid instability when X Y σ σ is very close to zero.
Numerical objects were simulated on a large grid in order to reduce the effect of finite support in the computation.Both ( , ) X Y Γ and ( , ) X Y ρ are computed over a region of interest (ROI) where the object intensities were concentrated.Both ( , )  X Y Γ and ( , ) X Y ρ are computed on the region of interest (ROI) of the images to reduce the effect of background on the performance metrics.The ROI of the images to be compared were registered properly in order to avoid unwanted bias in the results.Besides ( , ) X Y Γ and ( , ) X Y ρ , intensity profiles taken through the center of XZ sectional images are also compared for a 1-D spatial comparison of intensity distribution.

Results and discussion
In this Section we present different results obtained from the investigation studies that include: reduction of depth sensitivity of the proposed SQUBIC-WFE imaging system compared to the conventional imaging in terms of the PSFs of the respective cases.Effect of the strength parameter (A) of the SQUBIC phase mask on the depth sensitivity of the proposed system is also discussed.In addition, validation of the proposed method is demonstrated by restoring simulated and experimentally acquired images of the micro spherical shells at different depths.The XZ sectional images of the 3D volumes are displayed and analyzed as other planes do not carry additional information in this study.The images are shown on their own color scale to visualize the fine details; however, quantitative comparisons are based on the true values of the image intensities.

Effect of depth-induced SA on SQUBIC-WFE PSFs
Results that compare the depth sensitivity of the conventional PSFs with respect to the SQUBIC-PSFs, and demonstrate the reduced depth sensitivity of the SQUBIC-PSF, are summarized in Fig. 3.The simulated XZ sectional view images of the PSFs computed at depth 0 µm and 100 µm for a 63X/1.4NA oil immersion objective lens are shown in Fig. 3  (a-f).The conventional PSFs and SQUBIC-PSFs for A = 20 and A = 50 are shown in Fig. 3  (a-b), (c-d), and (d-f), respectively.It is observed that the intensity distributions of the conventional PSFs change more rapidly at different depths compared to the SQUBIC-PSFs [Fig.3(a-f)].The change of PSFs with depths is quantified [Fig.3 (g)] by the SSIM between the PSF computed at 0 µm depth and the PSFs computed at other depths.The change of SSIM is shown with respect to the depth for the conventional and SQUBIC-PSFs for different values of A. It is noticed that, the SSIM reduces to 0.16 in the case of conventional PSF at 100 µm depth, whereas in the case of SQUBIC-PSFs (A = 80) it only drops to 0.75.This observation indicates an 86% change in the PSFs over the range of 100 µm in the conventional case, whereas only 25% change in PSFs occur in the SQUBIC-WFE case; thus the WFE system shows 61% less variability compared to the conventional system.In addition, it is seen that, with the increase of the value of A, the change of the SSIM with respect to depth is reduced.

Performance evaluation of the SQUBIC-WFE through restoration of simulated images
To show performance of the SQUBIC-WFE imaging in terms of its ability to reduce the impact of depth-induced SA in the case of a 20X/0.8NA dry lens, simulated images of 3 µm in diameter spheres located at different depths are restored (object 1 as described in Sec.3.3) using a single PSF.The location of the sphere is varied from depth 20 µm to 40 µm at an interval of 2 µm, and images are restored by a PSF computed at 30 µm.The XZ sectional view of the object is shown in Fig. 4(a), and the correlation coefficients between the object and the restored images are shown in Fig. 4(b).The correlation coefficient as a function of depth quantifies the restoration performance of the conventional and the SQUBIC-WFE systems for these images.As the images are restored using the PSF computed at 30 µm, the object at 30 µm is restored with the maximum accuracy (highest correlation) in the both cases.The advantage of the SQUBIC system over the conventional is that, the restoration accuracy reduces slowly compared to the conventional system.This result characterizes the reduced depth sensitivity of the system (i.e. the correlation coefficient changes 29% in the case of conventional system, and 16% in the case of SQUBIC-WFE system).To compare the capability of recovering the exact locations of the beads along with the restoration accuracy of the 3D intensity distribution in the case of a thick sample, we restored a 40 µm thick simulated image with five 3 µm in diameter spherical beads located at depths between 10 µm to 50 µm at an interval of 10 µm (object 2) as shown in Fig. 5(a).Figure 5(b)  and c show the simulated forward conventional and SQUBIC-WFE images of the object in Fig. 4(a), respectively.From a qualitative comparison between the forward images, it is seen that in the conventional case, the images of the beads vary at different depths; whereas, in the case of SQUBIC-WFE, all the spheres appear to have the same shape.Figure 5 (d) and (e) show the restored conventional and SQUBIC-WFE images where respective PSFs computed at 30 µm depth are used in the SIEM restoration process.To quantify the restored images, an intensity profile is drawn through the center of the beads shown in Fig. 5(f).It is observed from the intensity profile that, in the both cases, the middle beads (beads at 30 µm depth) are restored the best (ρ = 0.97 and 0.96 for the conventional and SQUBIC-WFE, respectively) compared to the beads at other depth locations.Note that a small shift in the axial location of the restored beads is evident in all cases except for the case of the bead located at depth of 30 µm; this is because of the effective model mismatch when a single PSF computed at 30 µm depth is used in the restoration.The correlation coefficients between the true and restored microspheres at 10 µm and 50 µm depths are 0.79, 0.51 and 0.92, 0.83 in the case of conventional and SQUBIC-WFE system, respectively.Therefore, the variability in the correlation coefficients in respective cases are 0.46 and 0.13, which indicates 33% reduction in depth variability of the SQUBIC-WFE system compared to the conventional system.These results support that, in the case of SQUBIC-WFE, beads at 10 µm and 50 µm are also restored with better accuracy compared to that of the conventional cases.The intensity profile also supports the previous observation that, the location and shape of the bead at 30 µm is restored well in the both cases.The difference between the centers of the FWHM locations between the object and restored images are 2.4 µm and 1.5 µm for the conventional and SQUBIC-WFE system, respectively in the case of beads at both 10 µm and 50 µm depths.This observation suggests that object locations are better identified in the case of SQUBIC-WFE imaging system.

Evaluation of the SQUBIC-WFE imaging system forward model
We observe that the intermediate images at different depths vary more in case of the conventional system than that of the SQUBIC-WFE system.The comparison between the intermediate images in the case of conventional and the SQUBIC-WFE (design parameter A = 20) imaging is summarized in Fig. 6, where the XZ sectional views of the experimentally acquired [Fig.6 (a & d)] and simulated images [Fig.6 (b, c & e, f)] of a 6 µm in diameter spherical shell are shown.The center of the micro shell is located at 3 µm depth in the experimental case, and at 3 µm and 27 µm in the simulated images.We observe that, in the

Evaluation of the SQUBIC-WFE system in simulation
Restoration results from simulated images of the 6 µm in diameter spherical shell, in the case of conventional and SQUBIC-WFE imaging are summarized in Fig. 7.This study was performed for different values of the strength parameter A of the SQUBIC phase mask, i.e. namely A = 20, 50, and 80 in order to show the effect of the design parameter (A) value on the image restoration; however, restored images are shown only for the A = 80 case.The XZ sectional view of these images restored from depths 3 µm, 12 µm, and 27 µm using the PSFs computed at depths 3 µm, 12 µm, and 27 µm are shown in Fig. 7(a) and Fig. 7(b) for conventional and SQUBIC-WFE imaging, respectively.The correlation coefficients ( ρ ) values computed between the true object and the restored images are shown with the respective images.We observe that in the case of the conventional system, the images are restored best when a depth-matched PSF is used.However, artifacts appear when the images are restored using other PSFs.On the other hand, in the case of SQUBIC-WFE imaging, the spherical shape of the restored images is adequate (judged by the visual observation and quantified by the values of ρ ) not only when a depth-matched PSF is used, but also when the PSF does not correspond to the actual depth of the object.The overall variation of the value of ρ in the case of SQUBIC-WFE and conventional system are 0.06 and 0.17, respectively which suggests 11% overall improvement in the reduction of depth variability in SQUBIC-WFE system within 30 µm depth.
To quantify the differences between the restoration results from the different cases, the standard deviation of all the correlation coefficient values for all object depth cases, ( ) σ ρ ,was computed [Fig.7(c)].Note that, the smallest standard deviations are obtained in the case where the spherical shell is at depth 12 µm, because the relative depth mismatch between the object location and the depth at which the restoring PSF was computed is minimum at this location compared to the other cases.It is seen that, ( ) σ ρ as a function of object depth is smaller in the case of SQUBIC-WFE imaging compared to conventional imaging, and it decreases as the value of A increases.The mean value of the ( ) σ ρ is 0.1502 in the case of the conventional system and 0.0023 in the case of the SQUBIC-WFE (A = 80) system; this quantifies the amount of improvement in reducing the depth-sensitivity of the proposed imaging system.
To visualize the quantitative improvement in the image restoration, an intensity profile comparison between the conventional and the SQUBIC-WFE imaging is shown in Fig. 7(d).
The intensity profiles are drawn though the center of the XZ sectional images shown in Fig. 7 (a & b) where the correlation coefficients are the minimum in the respective cases (i.e. ).These images are chosen to compare the worst case result for both systems.It is seen from the intensity profiles that, although the true object intensity has only two peaks (between 0 µm to 6 µm in the axial direction), a third peak is observed due to an artifact in the conventional restored image intensity, which is not part of the object.On the other hand, in the case of SQUBIC-WFE restored image, although the location of the second peak does not coincide with the true location in the object, artifacts are minimized despite the use of a depth-mismatched PSF in the restoration.This result demonstrates that the SQUBIC-WFE system reduces the effect of depth aberration on restoration artifacts.

Experimental evaluation of the SQUBIC-WFE system
For this experimental study, an A = 20 SQUBIC phase mask was used due to limitations of the SLM used [15].The experimental restored images in the case of conventional and SQUBIC-WFE systems are shown in Fig. 8(a) and Fig. 8(b), respectively.The center of the beads are at depths d 0 , d 0 + 12 µm, and d 0 + 27 µm, where d 0 ~3 µm.All the images were restored using the unaberrated PSF (computed at 0 µm) and the PSFs at depths 12 µm and 27 µm.To quantify the depth sensitivity of the conventional and SQUBIC-WFE imaging system separately, the following restored images were used as reference to compute the correlation coefficients (ρ): in the conventional case, the image at depth d 0 restored with the PSF computed at 27 µm depth; and in the SQUBIC-WFE case, the image at d 0 + 12 µm restored using the PSF computed at 0 µm depth.It is observed from the values of ρ that the stucture of the object is restored better in all the cases of the SQUBIC-WFE images compared to the conventional images.Comparison of the values of ρ reveals that the change of ρ, ∆ρ is equal to 0.36 in the case of conventional and 0.07 in the case of SQUBIC-WFE imaging.This observation demostrates that, in the depth range of 30 µm, image restoration performed with a depth-mismatched PSF produces an error of 36% conventional imaging, while the error reduces to 7% in the case of SQUBIC-WFE imaging which indicates a 29% inprovemnt in the restoration.The means of the standard deviations beween the correlation coefficients at different depths are 0.0778 and 0.0183 in the case of conventional and SQUBIC-WFE images, respectively, which also reflect the reduced depth sensitivity achieved in the SQUBIC-WFE imaging system.It is worth mentioning that the restorations obtained in the case of our conventional imaging system [Fig.8(a)] are not consistent with the simulation [Fig.7(a)] due to the residual aberration (RA) in the MO lens described in Section 3.1.As evident in Fig. 8(a), the best restoration for depth d 0 , was obtained with the PSF at depth 27 µm because the RA is partially compensated by the aberrant PSF.On the contrary the RA does not affect restoration achieved with the proposed SQUBIC-WFE imaging system [as evident in Fig. 8(b)], which is shown to provide robust restoration over the investigated imaging conditions.

Summary and conclusion
In this study, a wavefront-encoded based imaging system was presented in which the pupil of the imaging system is modified by a SQUBIC phase mask.The WFE system reduces the effect of optical aberration on restoration by reducing the change of the PSF at different depths due to depth-induced SA.Different SQUBIC phase masks were implemented by varying its strength parameter A and their impact on restoration was investigated.
As the SQUBIC phase mask design depends on the imaging system parameters [Eq.( 1)], designs for both the low NA (20X/0.8NA dry) and high NA (63X/1.4NAoil immersion) objective lenses were investigated.Results from a PSF comparison [Fig.3(g)] support the theoretical prediction of reduced depth variability in the SQUBIC-WFE imaging system over the conventional system.The comparison was quantified using the SSIM, which was found to be 61% less for the SQUBIC (A = 80) WFE system compared to the conventional system (using a 63X/1.4NA oil-immersion) within a depth interval of 100 µm [Fig.3(g)].
Simulation results show that in the case of the proposed SQUBIC-WFE system [Fig.4 and Fig. 5], image restoration accuracy is less dependent on the depth of the PSF used for restoration, compared to the conventional system.The variability of the correlation coefficient used to quantify simulation results was found to be 33% less in the case of SQUBIC-WFE system over the conventional system, within a depth range of 20 to 40 µm.
The performance of the SQUBIC-WFE system in the case of the high NA objective lens was evaluated by restoring simulated and experimental images of samples with 6 µm in diameter spherical shells distributed at different depths below the coverslip.Simulation results shown here for phase mask design parameter A = 80 demonstrated the reduced depth sensitivity of the system as predicted by theory [Fig.7].In this simulation study, the variability in the error between the restored and true object intensities quantified by the correlation coefficients was found to be reduced by 11% in the WFE system over the conventional system.
Results from experimental data [Fig.8] were found to be consistent with simulations results.In this case, the variability in the error between the restored and true object intensities quantified by the correlation coefficients was found to be reduced by 29% in the WFE system over the conventional system over a depth range of 30 µm.
In summary, use of the proposed WFE system enables restoration with reduced artifacts over the conventional system of 3D images suffering from aberration, such as sampleinduced aberration in thick tissue, using a space-invariant image restoration algorithm and a single SQUBIC-WFE PSF.The proposed system reduces the necessity to precisely determine the imaging depth for the computation of PSFs used in restoration and it relaxes the requirement for high precision SA correction of microscope objective lenses.Future, studies will report on experimental validation of the proposed system with phase mask design parameters higher than A = 20 using a fabricated phase mask instead of the SLM-based phase mask presented here, which will allow further improvements over larger imaging depths (>30 µm investigated here) as predicted by theory and simulation.

Fig. 1 .
Fig. 1.SQUBIC phase mask.Unwrapped mask (gray level presentation of the phase values onto the pupil) for: (a) NA = 0.8, A = 88; (b) NA = 1.4,A = 20.(c) Phase depth for low NA and high NA MO; (d) and (e) show the wrapped phase masks corresponding to (a) and (b), respectively.

63X/ 1 A
A values (low, medium, and high) are chosen to show the effect of the design parameter (A) on the image restoration.is chosen accordingly in order to show model validation through experimental verification.SQUBIC phase mask with A = 20 implements accurately by the SLM.

Fig. 4 .
Fig. 4. Impact of system depth sensitivity on restoration of simulated images from SQUBIC-WFE and conventional systems.(a) XZ sectional view of a simulated object with a 3 µm in diameter sphere centered at 30 µm; (b) Correlation coefficients between the simulated true object and the corresponding restored images in the case of conventional and SQUBIC system (beads centered at depths, o d d m d = ± Δ where,

Fig. 5 .
Fig. 5. Reduction of image restoration artifacts while imaging deep into the sample using the SQUBIC-WFE imaging system in simulations: XZ cut view images from: (a) an object with five 3 µm in diameter spherical beads in a depth range from 10 µm to 50 µm; (b) Conventional simulated image; (c) SQUBIC-WFE simulated image; Restored images using the PSFs computed at 30 µm depth in the case of (d) conventional imaging, and (e) SQUBIC-WFE imaging.(f) Intensity profiles comparison through the center of the object and the restored images as show in Fig. (a).Lens: 20X/0.8NAair-immersion, specimen embedding medium is water (RI = 1.33), and the emission wavelength is 515 nm.SQUBIC design parameter A = 88.

Fig. 6 .
Fig. 6.Comparison between the conventional and SQUBIC-WFE intermediate images.The XZ sectional view images of a 6 µm in diameter spherical shell having shell thickness equal to 1 µm in the case of: (a)-(c) Conventional system, and (d)-(f) SQUBIC (A = 20) WFE system.(a) and (d) are experimentally acquired images from depth 3 µm; whereas (b-c) and (e-f) are corresponding simulated images at depths 3 µm and 27 µm, respectively.Lens: 63X/1.4NA oil immersion; Specimen embedding medium is prolong diamond (RI = 1.47 when cured); Emission wavelength: 515 nm.conventional case the intensity distributions of the images at depths 3 µm and 27 µm are different in simulation; on the other hand, in the SQUBIC-WFE case the intensity distribution of the simulated images not only match with the experiment, but the intensity distribution is also similar at depths 3 µm and 27 µm.

Fig. 7 .
Fig. 7. Performance comparison between the conventional and the SQUBIC-WFE system in terms of image restoration through simulations.(a) XZ sectional views of the conventional restored images of the spherical shells centered at depths 3 µm, 12 µm and 27 µm an restored using PSFs at 3 µm, 12 µm and 27 µm; (b) Corresponding restored images as (a) in the case of SQUBIC (A = 80) system.The values on the images are the correlation coefficients between the true object and the corresponding restored images.(c) Quantitative comparison in terms of the standard deviation ( ) σ ρ of the correlation values.Mean value of ( ) σ ρ is indicated by color coded (corresponding to each case) horizontal lines on the bar plot.(d) Intensity profile comparison between the conventional and the SQUBIC images at 27 µm, restored with the corresponding PSFs at 3 µm.Lens: 63X/1.4NA oil immersion; RI of the specimen embedding is equal to 1.47; emission wavelength is 515 nm.

Fig. 8 .
Fig. 8. Performance comparison between the conventional and the SQUBIC-WFE system through the experimental image restoration.(a) Restored images of three different objects located at depths d 0 µm, d 0 + 12 µm, and d 0 + 27 µm using the PSFs at 0 µm (unaberrated PSF), 12 µm and 27 µm depths.(b) Corresponding restored images in the case of SQUBIC-WFE system.The numbers on the images are the correlation coefficients.The corresponding best restored images (from a visual inspection) are used as reference images in each case to compute the correlations.Lens: 63X/1.4NA oil immersion; Specimen embedding medium is prolong diamond (RI = 1.47 when cured); Emission wavelength: 515 nm; Design parameter of the SQUBIC phase mask A = 20.