We present a formulation of optical point spread function based on a scaled three-dimensional Fourier transform expression of focal field distribution and the expansion of generalized aperture function. It provides an equivalent but more flexible representation compared with the analytic expression of the extended Nijboer-Zernike approach. A phase diversity algorithm combined with an appropriate regularization strategy is derived and analyzed to demonstrate the effectiveness of the presented formulation for phase retrieval and deconvolution. Experimental results validate the performance of presented algorithm.
© 2012 OSA
In most cases of focus-diverse phase retrieval and deconvolution, such as optical metrology, adaptive optics, and deconvolution microscopy [1–6], the optical point spread function (PSF) is usually modeled by using the Fresnel or Fraunhofer diffraction approximation, which can be expressed as a two-dimensional Fourier transform. It is a convenient representation for the computation providing that the accuracy of approximation is sufficient for the applications. For the cases of high numerical aperture (NA) and large Fresnel number (FN), the Debye diffraction approximation can be adopted in some forms to take into account of vectorial effects . For different defocused planes, the PSF will be computed separately for both two-dimensional and three-dimensional applications.
Based on the Debye diffraction approximation and the Zernike expansion of generalized pupil function, the extended Nijboer-Zernike (ENZ) approach derives the analytic expression of through-focus two-dimensional PSF [8–10]. For the cases of high NA and large FN, the ENZ approach agree well with the more rigorous Rayleigh-Sommerfeld-I (RSI) diffraction integral providing that enough expansion terms are used for approximation. The analytic PSF of the ENZ approach has shown its effectiveness in the applications of both phase retrieval and deconvolution [10,11].
It is also known that the Debye approximation of focal field distribution can be expressed as a three-dimensional Fourier transform of the generalized aperture function , by which the through-focus two-dimensional PSFs can be computed with a single three-dimensional Fourier transform. This is a flexible approach to represent the optical PSF, but it is not valid for the cases of low NA or small FN . Recently, a paraxial approximation of focal field distribution is proposed by using a scaled three-dimensional Fourier transform, which can provide accurate results for the cases of low NA and small FN . We will show that for the computation of PSF this approximation is consistent with the corrected Debye approximation by introducing an axial and lateral coordinate transformation according to the concept of axial optical coordinate. Therefore, the Fourier transform expression of the Debye approximation can be effectively extended to the cases of low NA and small FN.
In this paper, we present a formulation of optical PSF based on the three-dimensional Fourier transform expression of focal field distribution and the expansion of generalized aperture function. This formulation is equivalent to the analytic expression of the ENZ approach, but its representation is more flexible than that of the ENZ approach since different aperture functions and expansion polynomials can be taken into account other than the case of circular aperture and Zernike polynomials. A phase diversity algorithm combined with an appropriate regularization strategy is derived and analyzed to demonstrate the effectiveness of the presented formulation for phase retrieval and deconvolution, which can realize the rather high accuracy and stability even for large aberrations through the experimental validation.
The rest of this paper is organized as follows. Section 2 presents the formulation of three-dimensional PSF, including a brief comparison with the analytic expression of the ENZ approach. Section 3 describes the derivation and analysis of the phase diversity algorithm by using the three-dimensional PSF. Section 4 gives the experimental results for validating the performance of presented algorithm. Section 5 concludes the paper.
2. The three-dimensional PSF
2.1 Three-dimensional focal field distribution
As shown in Fig. 1 , a three-dimensional space where the optical field distribution in the focal region is described in the Cartesian coordinates (x,y,z). A circular aperture W lies in the (u,v) plane that is parallel to the (x,y) plane at a distance d from the focus F which is defined as the origin of the coordinates. When a convergent spherical wave passes through the aperture and propagates to the focal region, the complex amplitude at P is given by the Huygens-Fresnel principle 10,12], the focal field distribution can be expressed in the Cartesian coordinates as
From Eq. (2), it can be seen that the focal field distribution can be expressed as a scaled three-dimensional Fourier transform of the generalized aperture function, but this expression based on the Debye approximation is only valid for the cases of high NA and large FN. We can apply a correcting coordinate transformation [16,17] to Eq. (2) for the cases of low NA or small FN, which is a simple scaling operation to the axial and lateral coordinates of focal field and can be expressed as
To be more explicitly, the coordinate transformation in Eq. (5) can be expressed in the Cartesian coordinates as
By using the above correction to approximate the asymmetric axial field and focal shift which is derived according to the concept of axial optical coordinate , the corrected Debye approximation shows a better agreement with the RSI diffraction integral for the cases of low NA and small FN. This coordinate transformation is the same as the treatment of a paraxial approximation of focal field distribution for the cases of low NA and small FN , which also derives a scaled three-dimensional Fourier transform expression of focal field distribution by using the correction of Eq. (7) as
Comparing the expression of Eq. (8) with that of corrected Eq. (2), we can see that the only difference is the scaling factors of the three-dimensional Fourier transform. When the defocused distance is small, the difference can be neglected for the computation of PSF. We will not differentiate two expressions and denote the scaling factor as g(x,y,z) in the following derivation. With the correcting coordinate transformation, the Debye approximation expressed as a scaled three-dimensional Fourier transform can be effectively extended to the cases of low NA and small FN.
2.2 Formulation of three-dimensional PSF
For the notational convenience, we rewrite the focal field distribution as10], we can expand the generalized aperture function P(u,v,w) with the Zernike polynomials
Comparing with the analytic expression of normalized focal field distribution by using the ENZ approach in the cylindrical coordinates
Although the analytic expression of the ENZ approach will maintain a higher numerical accuracy for the computation, the Fourier transform approach is more flexible to take in account of different aperture functions and adopt other expansion polynomials. For example, for the cases of annular aperture, the Fourier transform approach can approximate the focal field distribution more appropriately with the annular aperture function and the corresponding annular Zernike polynomials. By using the expression of focal field distribution in Eq. (12), the three-dimensional PSF can be represented as
Therefore, the formulation of through-focus PSF is obtained by using the scaled three-dimensional Fourier transform expression of focal field distribution and the Zernike expansion of the generalized aperture function. This formulation is accurate providing that enough expansion terms are used for approximation and is flexible for the extension by adopting other aperture functions and expansion polynomials.
The formulation of three-dimensional PSF provides a convenient representation for the applications of phase retrieval and deconvolution. On one hand, by the direct expansion of generalized aperture function rather than its phase, the amplitude and phase of generalized aperture function can be parameterized and estimated simultaneously. It has been shown that the phase estimation can be improved by the simultaneous amplitude estimation when the aperture amplitude variations are present . This treatment can cope with the aperture amplitude variations due to the scintillation effects, etc. On the other hand, the phase estimation of generalized aperture function can be improved by the phase unwrapping, therefore the estimation can be less susceptible to the wrapping-induced phase ambiguity for large aberrations.
3. The phase diversity algorithm
3.1 Description of the algorithm
For incoherent imaging, the three-dimensional model of image formation can be expressed as the convolution of object and PSF associated with the noise20] to demonstrate the effectiveness of the formulation of three-dimensional PSF, although other approaches such as the maximum likelihood or maximum a posteriori can also be employed [3,6].
Since the amplitude and phase of generalized aperture function can be simultaneously estimated by using the formulation of three-dimensional PSF and the aperture amplitude can be known roughly according to the prior knowledge, bounding the aperture amplitude variations will be an appropriate regularization strategy to stabilize the nonlinear optimization procedure. Based on the imaging model of Eq. (18), we express the error metric for the optical inverse problem of phase diversity in the Fourier domain as6]. For the opaque objects that are usually considered in the two-dimensional cases, we can regard that the object only lies in the plane of z = 0. Therefore, another error metric can be obtained for the two-dimensional opaque objects as
For the notational convenience in the following derivation, the coordinate variables are omitted, and the single index is used for the Zernike polynomials. To adapt the formulation of three-dimensional PSF to the above framework of phase diversity, we reformulate the PSF expression in Eq. (17) as
From now on, we will only give the derivations of two-dimensional case, the expressions of three-dimensional case can be derived accordingly. From the expression of Eq. (26), the gradient and Hessian of PSF spectrum with respect to the real and imaginary parts of β can be obtained straightforwardlyEq. (10), |P|2 can be represented as
For further simplifying the representations, the object-independent error metric in Eq. (22) is expressed as
The gradient and Hessian of |T|2, D and Q2 can be expressed using those of S and |P|2 as
With the gradient and Hessian of error metric, the Newton’s method of locally quadratic convergence can be employed to improve the performance of nonlinear optimization. Using the estimated β coefficients, the generalized aperture function can be obtained. Then the phase is retrieved by unwrapping the aperture phase and the object is reconstructed by using Eq. (24).
By using the formulation of three-dimensional PSF, the above algorithm has certain computational advantages since the basis functions Wkl and Xkl and as well as the Hessian of S and |P|2 can be precomputed. In addition, the regular structural properties of the optimization problem expression can also be exploited to speed the optimization procedure.
3.2 Implementation and analysis
We implement and analyze the above phase diversity algorithm in the MATLAB environment. The algorithm is implemented by using the trust-region Newton-CG method for nonlinear optimization  and the Goldstein’s branch cut method for phase unwrapping , where the chirp z-transform is used to perform the three-dimensional Fourier transform . For the computational convenience, the axial and lateral coordinates of focal field distribution are normalized with the axial unit 2λ/(NA)2 and lateral unit λ/(NA), respectively. Both the circular and annular apertures are used with the initial estimation of aperture amplitude is set to the aperture function. In the nonlinear optimization, the initial β coefficients are set to zero except that β0 = 1, and the regularization parameters are tuned empirically.
Numerical simulations are carried out to analyze the algorithm for both two-dimensional and three-dimensional extended scenes, especially in the cases of low NA and small FN. In the simulations, the RSI diffraction integral is used to propagate the aberrant optical field to comupte the PSFs. Then the object is convolved with the PSFs to form the diversity images. Finally, the object and optical field at the aperture is estimated from the diversity images. The aperture phase aberrations are randomly generated by using the low-order 10 terms of Zernike polynomials except the piston with the uniform distribution. The aperture size is set to 6mm, both the optical fields and the objects are sampled with 64 × 64 pixels, and the diversity images are Nyquist sampled, which corresponds to the size of 16 lateral units.
For the two-dimensional cases, the simulated object is the resolution test chart, and only two diversity images are used including the focused image and the defocused image with the defocused distance of 1 axial units. Typical simulation results with the parameters of FN = 50 and NA = 0.01 are shown in Fig. 2 , where the first 36 terms of both circular and annular Zernike polynomials are adopted for the expansion of generalized aperture function. For comparison, the corresponding results without the aperture amplitude regularization (μ = 0) are also given. In the simulations, we find that the amplitude regularization plays an important role for the accuracy of presented algorithm even without the amplitude variations, which can be seen obviously from Fig. 2. In addition, it also shows that the cases of annular aperture can be accurately approximated with the corresponding annular Zernike polynomials.
For the three-dimensional cases, the simulated object is a stack of five slices cropped from the resolution test chart with the same size of 16 lateral units, which is padded with two empty slices before the first slice and after the last slice, respectively. The distances between the adjacent slices are set to 1 axial units. The object stack is convolved with the simulated three-dimensional PSF to form the image stack with nine slices, including five focused images and four defocus images with the defocused distances of ± 1 and ± 2 axial units, respectively. The Z support constraint according to the prior knowledge  is incorporated to restrict the optimization procedure. Typical simulation results with the parameters of FN = 50 and NA = 0.01 are shown in Fig. 3 , where the aperture is circular and the first 36 terms of Zernike polynomials are adopted for the expansion of generalized aperture function. The effectiveness of presented formulation for this simple three-dimensional problem can be seen from Fig. 3, but there are still noticeable residual ghosts in the background and high frequency artifacts superimposed on the object. To further improve the reconstructed results, more appropriate object models should be considered.
In the cases of small phase aberrations, there is a coarse relationship Im(β)≈α between the complex Zernike coefficients β scaled by β0 and the ordinary Zernike coefficients α . This is also justified in the simulations, which suggests that using certain terms of β comparable to that of α should accurately represent the generalized aperture function. But for large phase aberrations, we find that the above relationship does not hold. As shown in Fig. 4 , by direct decomposition of a generalized aperture function with uniform amplitude and large phase aberrations, whose phase is generated with the random low-order α coefficients, the scaled β coefficients distribute over a rather broad range and the above relationship is not valid anymore. Therefore, enough expansion terms must be used for adequate approximation of large aberrations.
4. Experimental setup and results
In addition to the simulations, we use a simple experimental setup to validate the performance of presented algorithm for two-dimensional extended scenes in the case of low NA and small FN. As shown in Fig. 5 , the collimated beam with λ = 632.8nm emitted from a ZYGO laser interferometer illuminates a reflective liquid crystal spatial light modulator (LC-SLM) from BNS through a beam splitter. The LC-SLM with 256 × 256 pixels and 6.14mm × 6.14mm array size can accurately generate low-order phase aberrations. Large phase aberrations can also be generated by means of the phase wrapping method. For the phase-only modulation, the LC-SLM is adjusted to align the vertical axis to the polarization direction of the incident beam. A circular aperture stop of diameter 6mm is positioned in front of the LC-SLM, and the reflected beam through the aperture stop is sent to the interferometer for phase measurement and to a camera for PSF measurement. The camera is an ANDOR DU-888 EMCCD camera with 14 bit depth and 13μm × 13μm pixel size, which is mounted on a motorized translate stage for the defocus movement.
In the experiments, random phase aberrations are generated with the LC-SLM by using the low-order 10 terms of Zernike polynomials except the piston. The PSF is sampled with 128 × 128 pixels and the size of 16 lateral units (104 × 104 pixels) is used for the computation, where the marginal pixels are the guard band for image registration. The nominal focal length is 400mm, which corresponds to the parameters of NA = 0.0075 and FN = 35.56. Five diversity PSFs including the focused PSF and four defocused PSFs with the defocused distances of ± 1 and ± 2 axial units are measured for each phase aberration. After preprocessing and registration, the measured PSFs are convolved with the simulated object to form the diversity images. Then the diversity images are used for phase retrieval and deconvolution, where the first 136 terms (15 orders) of Zernike polynomials are adopted for expansion. The sizes of both simulated and reconstructed objects are 128 × 128 pixels. The optical field is also sampled with 128 × 128 pixels, which corresponds to 47μm sampling interval for 6mm aperture size. For comparison, the measured phases by the interferometer are resampled and registered to the corresponding reconstructed phases for computing the residual errors. There are six marginal pixels are trimmed for both reconstructed and measured phases to avoid the mismatch at the ragged edges. Some experimental results for small to large phase aberrations are summarized in Table 1 and shown in Fig. 6 , where each data point is averaged over 10 random phase aberrations with roughly the same RMS, and the error bar represents the data range. Typical results of phase retrieval and deconvolution are given in Figs. 7 -8 .
From the experimental results, it can be seen that the rather high accuracy and stability can be realized by the algorithm for the low-order phase aberrations up to 1λ RMS, where the residual phase errors can be kept almost within the diffraction limit. The experimental results can also demonstrate the effectiveness of the formulation of three-dimensional PSF for phase retrieval and deconvolution. More appropriate object and noise models should be considered to improve the algorithmic performance especially under the low signal-noise-ratio conditions, and further experimental validations will be carried out for the practical applications.
In this paper, we present an accurate and flexible formulation of optical PSF based on a scaled three-dimensional Fourier transform and the expansion of generalized aperture function. The formulation is equivalent to the analytic expression of the ENZ approach but is more flexible to take into account of different aperture functions and expansion polynomials. It provides a convenient representation for the applications of phase retrieval and deconvolution. A phase diversity algorithm combined with an appropriate regularization strategy is derived to demonstrate the effectiveness of the formulation of three-dimensional PSF and is analyzed and validated for two-dimensional extended scenes by simulations and experiments, which can realize rather the high accuracy and stability even for large aberrations. The algorithm based on a three-dimensional imaging model is a general approach for both two-dimensional and three-dimensional applications.
References and links
2. R. A. Gonsalves, “Phase retrieval and diversity in adaptive optics,” Opt. Eng. 21, 829–832 (1982).
3. R. G. Paxman, T. J. Schulz, and J. R. Fienup, “Joint estimation of object and aberrations using phase diversity,” J. Opt. Soc. Am. A 9(7), 1072–1085 (1992). [CrossRef]
5. J. B. Sibarita, “Deconvolution Microscopy,” Adv. Biochem. Eng. Biotechnol. 95, 201–243 (2005). [PubMed]
6. G. Chenegros, L. M. Mugnier, F. Lacombe, and M. Glanc, “3D phase diversity: a myopic deconvolution method for short-exposure images: application to retinal imaging,” J. Opt. Soc. Am. A 24(5), 1349–1357 (2007). [CrossRef] [PubMed]
9. J. J. M. Braat, P. Dirksen, and A. J. E. M. Janssen, “Assessment of an extended Nijboer-Zernike approach for the computation of optical point-spread functions,” J. Opt. Soc. Am. A 19(5), 858–870 (2002). [CrossRef] [PubMed]
10. J. J. M. Braat, S. van Haver, A. J. E. M. Janssen, and P. Dirksen, “Assessment of optical systems by means of point-spread functions,” in Progress in Optics, E. Wolf, ed., (Elsevier, 2008), 51, 349–468.
11. A. A. Ramos and A. L. Ariste, “Image reconstruction with analytical point spread functions,” Astron. Astrophys . 518, paper A6 (2010).
12. C. W. McCutchen, “Generalized Aperture and the Three-Dimensional Diffraction Image,” J. Opt. Soc. Am. 54(2), 240–242 (1964). [CrossRef]
13. E. Wolf and Y. Li, “Conditions for the validity of the Debye integral representation of focused fields,” Opt. Commun. 39(4), 205–210 (1981). [CrossRef]
14. J. Lin, X.-C. Yuan, S. S. Kou, C. J. R. Sheppard, O. G. Rodríguez-Herrera, and J. C. Dainty, “Direct calculation of a three-dimensional diffracted field,” Opt. Lett. 36(8), 1341–1343 (2011). [CrossRef] [PubMed]
15. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press, 1999).
16. S. van Haver, “The Extended Nijboer-Zernike diffraction theory and its applications,” PhD Dissertation, Delft University of Technology (2010).
17. J. J. M. Braat, S. van Haver, and S. F. Pereira, “Microlens quality assessment using the Extended Nijboer-Zernike (ENZ) diffraction theory,” presented at EOS Optical Microsystems, Capri, Italy, 27–30 Sept. 2009.
20. C. R. Vogel, T. Chan, and R. Plemmons, “Fast algorithms for phase diversity-based blind deconvolution,” Proc. SPIE 3353, 994–1005 (1998). [CrossRef]
21. J. Nocedal and S. J. Wright, Numerical Optimization 2nd ed. (Springer, 2006).
22. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley-Interscience, 1998).
23. M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus field calculations,” Opt. Express 14(23), 11277–11291 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-23-11277. [CrossRef] [PubMed]