A three-dimensional point spread function for phase retrieval and deconvolution

Xinyue Liu; Liang Wang; Jianli Wang; Haoran Meng

doi:10.1364/OE.20.015392

1. Introduction

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 [7]. 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 [12], 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 [13]. 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 [14]. 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 [15]

U (P) = - \frac{i A}{λ} \frac{\exp (- i k f)}{f} \iint_{W} U (Q) \frac{\exp (i k r)}{r} \cos θ d S,

where λ is the wavelength of the monochromatic light and k is the wavenumber. By using a paraxial approximation with the cosθ = d/f, which is equivalent to the Debye approximation [10,12], the focal field distribution can be expressed in the Cartesian coordinates as

U_{F} (x, y, z) = - \frac{i A}{λ} \frac{d}{f^{3}} ∭_{V} U_{W} (u, v, w) A (u, v, w) \exp (- i k \frac{x u + y v + z w}{f}) d u d v d w,

where A(u,v,w) is the three-dimensional aperture function and is nonzero only on the spherical cap W, which can be defined by

A (u, v, w) = δ (\sqrt{u^{2} + v^{2} + w^{2}} - f),

where δ(●) is the Kronecker delta function. Like the two-dimensional case, we can define the generalized aperture function as

Fig. 1 Geometry for computing the optical field distribution in the focal region.

Download Full Size | PDF

P (u, v, w) = U_{W} (u, v, w) A (u, v, w) .

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

\bar{z}' = \bar{z} / (1 + \frac{\bar{z} λ}{f π {(NA)}^{2}}), \bar{r}' = \bar{r} / (1 + \frac{\bar{z} λ}{f π {(NA)}^{2}}),

where

\bar{z}

and

\bar{r}

are the cylindrical coordinates of focal field normalized with the axial unit λ/π(NA)² and lateral unit λ/2π(NA) respectively as

\bar{z} = z \frac{π {(NA)}^{2}}{λ}, \bar{r} = \sqrt{x^{2} + y^{2}} \frac{2 π (NA)}{λ} .

To be more explicitly, the coordinate transformation in Eq. (5) can be expressed in the Cartesian coordinates as

x' = \frac{f}{f + z} x, y' = \frac{f}{f + z} y, z' = \frac{f}{f + z} z .

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 [18], 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 [14], which also derives a scaled three-dimensional Fourier transform expression of focal field distribution by using the correction of Eq. (7) as

\begin{array}{l} U_{F} (x', y', z') = - \frac{i A}{λ} \frac{f - z'}{f^{3}} \exp (i k \frac{x'^{2} + y'^{2} + 2 z'^{2}}{2 (f - z')}) \\ \times ∭_{V} U_{W} (u, v, w) A (u, v, w) \exp (- i k \frac{x' u + y' v + z' w}{f}) d u d v d w . \end{array}

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 as

U (x, y, z) = g (x, y, z) F^{3} [P (u, v, w)],

where F³[●] denotes the three-dimensional Fourier transform. Following the treatment of the ENZ approach [10], we can expand the generalized aperture function P(u,v,w) with the Zernike polynomials

P (u, v, w) = \sum_{n m} β_{n}^{m} Z'_{n}^{m} (u, v, w),

where β are the complex Zernike coefficients, and Z'(u,v,w) are the three-dimensional Zernike polynomials on the spherical cap, which can be defined as

Z'_{n}^{m} (u, v, w) = Z_{n}^{m} (u, v) A (u, v, w),

where Z(u,v) are the ordinary two-dimensional Zernike polynomials. With this expansion, we can express the focal field distribution as

U (x, y, z) = \sum_{n m} β_{n}^{m} V'_{n}^{m} (x, y, z),

where V'(x,y,z) are the three-dimensional basis functions for the expansion of focal field distribution, which are defined by

V'_{n}^{m} (x, y, z) = g (x, y, z) F^{3} [Z'_{n}^{m} (u, v, w)] .

Comparing with the analytic expression of normalized focal field distribution by using the ENZ approach in the cylindrical coordinates

U (r, φ, z) = h (r, φ, z) \sum_{n m} 2 i^{m} β_{n}^{m} V_{n}^{m} (r, z) \exp (i m φ),

where h(r,φ,z) is the normalizing factor, and the two-dimensional basis functions V(r,z) can be approximated analytically, which are defined as

V_{n}^{m} (r, z) = \int_{0}^{1} \exp (i z ρ^{2}) R_{n}^{| m |} (ρ) J_{m} (2 π r ρ) ρ d ρ,

where R(●) are the radial parts of the Zernike polynomials and J(●) are the Bessel functions of the first kind. It is obvious that the two representations are equivalent, and the relationship between the two sets of basis functions can be expressed explicitly

V'_{n}^{m} (r, φ, z) = 2 i^{m} h (r, φ, z) V_{n}^{m} (r, z) \exp (i m φ) .

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

S (x, y, z) = {| \sum_{n m} β_{n}^{m} V'_{n}^{m} (x, y, z) |}^{2} .

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 [19]. 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 noise

I (x, y, z) = O (x, y, z) \otimes S (x, y, z) + N (x, y, z),

where I(x,y,z) is the image, N(x,y,z) is the noise, O(x,y,z) is the object, and ⊗ denotes the three-dimensional convolution operator. In this paper, we will use the regularized least-square approach for the optical inverse problem of phase diversity like in [20] 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 as

\begin{array}{l} E_{3 D} = \frac{1}{2} \sum_{u v w} {| \tilde{I} (u, v, w) - \tilde{O} (u, v, w) \tilde{S} (u, v, w) |}^{2} \\ + \frac{γ}{2} \sum_{u v w} {| \tilde{O} (u, v, w) |}^{2} + μ \sum_{u v w} {[{| P (u, v, w) |}^{2} - {| P_{0} (u, v, w) |}^{2}]}^{2}, \end{array}

where ~denotes the Fourier spectrum, |P₀(u,v,w)| is the initial estimation of aperture amplitude according to the prior knowledge, and γ and μ are the object and amplitude regularization parameters, respectively. This is a general model for three-dimensional phase diversity, additional constraints can also be incorporated according to the prior knowledge such as the Z support constraint [6]. 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

\begin{array}{l} E_{2 D} = \frac{1}{2} \sum_{u v} \sum_{z} {| \tilde{I_{z}} (u, v) - \tilde{O} (u, v) \tilde{S_{z}} (u, v) |}^{2} \\ + \frac{γ}{2} \sum_{u v} {| \tilde{O} (u, v) |}^{2} + μ \sum_{u v w} {[{| P (u, v, w) |}^{2} - {| P_{0} (u, v, w) |}^{2}]}^{2} . \end{array}

where S_z is the two-dimensional PSF in a plane of defocused distance z. Since the above error metrics are the quadratic convex functions of the object spectrum, they can be reduced to the object-independent forms by the minimization with respect to the object spectrum as

\begin{array}{l} E'_{3 D} = - \sum_{u v w} \frac{{| {\tilde{I}}^{*} (u, v, w) \tilde{S} (u, v, w) |}^{2}}{γ + {| \tilde{S} (u, v, w) |}^{2}} \\ + \sum_{u v w} {| \tilde{I} (u, v, w) |}^{2} + μ \sum_{u v w} {[{| P (u, v, w) |}^{2} - {| P_{0} (u, v, w) |}^{2}]}^{2}, \end{array}

\begin{array}{l} E'_{2 D} = - \sum_{u v} \frac{\sum_{z} {| {\tilde{I_{z}}}^{*} (u, v) \tilde{S_{z}} (u, v) |}^{2}}{γ + \sum_{z} {| \tilde{S_{z}} (u, v) |}^{2}} \\ + \sum_{u v} \sum_{z} {| \tilde{I_{z}} (u, v) |}^{2} + μ \sum_{u v w} {[{| P (u, v, w) |}^{2} - {| P_{0} (u, v, w) |}^{2}]}^{2}, \end{array}

and the minimizers are expressed by

{\tilde{O}}_{3 D} (u, v, w) = \frac{\tilde{I} (u, v, w) {\tilde{S}}^{*} (u, v, w)}{γ + {| \tilde{S} (u, v, w) |}^{2}},

{\tilde{O}}_{2 D} (u, v) = \frac{\sum_{z} \tilde{I_{z}} (u, v) {\tilde{S_{z}}}^{*} (u, v)}{γ + \sum_{z} {| \tilde{S_{z}} (u, v) |}^{2}} .

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

S = \sum_{k l} β_{k} β_{l}^{*} V'_{k} V'_{l}^{*} = \sum_{k l} β_{k} β_{l}^{*} W_{k l}, W_{k l} = V'_{k} V'_{l}^{*},

where * denotes the complex conjugate and W_kl are the basis functions for the expansion of three-dimensional PSF. With this expression, the Fourier spectrum of both three-dimensional and two-dimensional PSF can be represented as

\tilde{S} = \sum_{k l} β_{k} β_{l}^{*} \tilde{W_{k l}}, \tilde{S_{z}} = \sum_{k l} β_{k} β_{l}^{*} \tilde{W_{k l} (z)} .

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 straightforwardly

\begin{array}{l} \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Re}} = \sum_{l} β_{l}^{*} \tilde{W_{k l} (z)} + \sum_{l} β_{l} \tilde{W_{l k} (z)}, \\ \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Im}} = i \sum_{l} β_{l}^{*} \tilde{W_{k l} (z)} - i \sum_{l} β_{l} \tilde{W_{l k} (z)}, \end{array}

\begin{array}{l} \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Re} \partial β_{l}^{Re}} = \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Im} \partial β_{l}^{Im}} = \tilde{W_{l k} (z)} + \tilde{W_{k l} (z)}, \\ \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Re} \partial β_{l}^{Im}} = - \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Im} \partial β_{l}^{Re}} = i \tilde{W_{l k} (z)} - i \tilde{W_{k l} (z)} . \end{array}

where Re and Im denote the real part and imaginary part, respectively. Using Eq. (10), |P|² can be represented as

{| P |}^{2} = \sum_{k l} β_{k} β_{l}^{*} Z'_{k} Z'_{l} = \sum_{k l} β_{k} β_{l}^{*} X_{k l}, X_{k l} = Z'_{k} Z'_{l},

where X_kl are the basis functions for the expansion of |P|². The gradient and Hessian of |P|² with respect to the real and imaginary parts can be obtained accordingly

\frac{\partial {| P |}^{2}}{\partial β_{k}^{Re}} = 2 \sum_{l} β_{l}^{Re} X_{k l}, \frac{\partial {| P |}^{2}}{\partial β_{k}^{Im}} = 2 \sum_{l} β_{l}^{Im} X_{k l},

\frac{\partial^{2} {| P |}^{2}}{\partial β_{k}^{Re} \partial β_{l}^{Re}} = \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Im} \partial β_{l}^{Im}} = 2 X_{k l}, \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Re} \partial β_{l}^{Im}} = \frac{\partial^{2} {| P |}^{2}}{\partial β_{k}^{Im} \partial β_{l}^{Re}} = 0 .

For further simplifying the representations, the object-independent error metric in Eq. (22) is expressed as

\begin{array}{l} E'_{2 D} = - \sum_{u v} \frac{{| T |}^{2}}{D} + \sum_{u v} \sum_{z} {| \tilde{I_{z}} |}^{2} + μ \sum_{u v w} Q^{2}, \\ T = {\tilde{I_{z}}}^{*} \tilde{S_{z}}, D = γ + \sum_{z} {| \tilde{S_{z}} |}^{2}, Q = {| P |}^{2} - {| P_{0} |}^{2} . \end{array}

The gradient and Hessian of |T|², D and Q² can be expressed using those of S and |P|² as

\begin{array}{l} \frac{\partial {| T |}^{2}}{\partial β_{k}^{Re, Im}} = 2 Re [T^{*} \sum_{z} {\tilde{I_{z}}}^{*} \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Re, Im}}], \\ \frac{\partial^{2} {| T |}^{2}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} = 2 Re [T^{*} \sum_{z} {\tilde{I_{z}}}^{*} \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} + \sum_{z} {\tilde{I_{z}}}^{*} \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Re, Im}} \sum_{z} \tilde{I_{z}} \frac{\partial {\tilde{S_{z}}}^{*}}{\partial β_{l}^{Re, Im}}], \end{array}

\begin{array}{l} \frac{\partial D}{\partial β_{k}^{Re, Im}} = 2 Re [\sum_{z} {\tilde{S_{z}}}^{*} \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Re, Im}}], \\ \frac{\partial^{2} D}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} = 2 Re [\sum_{z} ({\tilde{S_{z}}}^{*} \frac{\partial^{2} \tilde{S_{z}}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} + \frac{\partial \tilde{S_{z}}}{\partial β_{k}^{Re, Im}} \frac{\partial {\tilde{S_{z}}}^{*}}{\partial β_{l}^{Re, Im}})], \end{array}

\begin{array}{l} \frac{\partial Q^{2}}{\partial β_{k}^{Re, Im}} = 2 [Q \frac{\partial {| P |}^{2}}{\partial β_{k}^{Re, Im}}], \\ \frac{\partial^{2} Q^{2}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} = 2 [Q \frac{\partial^{2} {| P |}^{2}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} + \frac{\partial {| P |}^{2}}{\partial β_{k}^{Re, Im}} \frac{\partial {| P |}^{2}}{\partial β_{l}^{Re, Im}}] . \end{array}

After lengthy but straightforward manipulations, the gradient and Hessian of error metric in Eq. (32) can be obtained by using Eq. (33-35)

\frac{\partial E'_{2 D}}{\partial β_{k}^{Re, Im}} = \sum_{u v} [\frac{{| T |}^{2}}{D^{2}} \frac{\partial D}{\partial β_{k}^{Re, Im}} - \frac{1}{D} \frac{\partial {| T |}^{2}}{\partial β_{k}^{Re, Im}} + μ \frac{\partial Q^{2}}{\partial β_{k}^{Re, Im}}],

\frac{\partial^{2} E'_{2 D}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} = \sum_{u v} [\begin{array}{l} - 2 \frac{{| T |}^{2}}{D^{3}} \frac{\partial D}{\partial β_{k}^{Re, Im}} \frac{\partial D}{\partial β_{l}^{Re, Im}} + \frac{{| T |}^{2}}{D^{2}} \frac{\partial^{2} D}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} - \frac{1}{D} \frac{\partial^{2} {| T |}^{2}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} \\ + \frac{1}{D^{2}} \frac{\partial D}{\partial β_{k}^{Re, Im}} \frac{\partial {| T |}^{2}}{\partial β_{l}^{Re, Im}} + \frac{1}{D^{2}} \frac{\partial {| T |}^{2}}{\partial β_{k}^{Re, Im}} \frac{\partial D}{\partial β_{l}^{Re, Im}} + μ \frac{\partial^{2} Q^{2}}{\partial β_{k}^{Re, Im} \partial β_{l}^{Re, Im}} \end{array}] .

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 W_kl and X_kl and as well as the Hessian of S and |P|² 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 [21] and the Goldstein’s branch cut method for phase unwrapping [22], where the chirp z-transform is used to perform the three-dimensional Fourier transform [23]. For the computational convenience, the axial and lateral coordinates of focal field distribution are normalized with the axial unit 2λ/(NA)² 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 β₀ = 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.

Fig. 2 Simulation results for the two-dimensional phase retrieval (all units are the wavelength), where (a-c) are the cases of circular aperture and (d-f) are the cases of annular aperture. (a) and (d) are the randomly generated phase aberrations with the same RMS = 0.222λ. (b) and (e) are the reconstructed phases with amplitude regularization, where the RMS of residual errors are 0.023λ and 0.032λ, respectively. (c) and (f) are the reconstructed phases without amplitude regularization, where the RMS of residual errors are 0.031λ and 0.068λ, respectively.

Download Full Size | PDF

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 [6] 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.

Fig. 3 Simulation results for the three-dimensional deconvolution, where (a) is the stack of the simulated object, (b) is the stack of the focused image, and (c) is the stack of the reconstructed object. The RMS of randomly generated phase aberrations is 0.182λ, and the RMS of residual error of reconstructed phase is 0.07λ.

Download Full Size | PDF

In the cases of small phase aberrations, there is a coarse relationship Im(β)≈α between the complex Zernike coefficients β scaled by β₀ and the ordinary Zernike coefficients α [11]. 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.

Fig. 4 Zernike coefficients for a generalized aperture function with uniform amplitude and large phase aberration of RMS = 1.398λ, where the charts from top to bottom are the α coefficients and the corresponding imaginary and real parts of β coefficients, respectively.

Download Full Size | PDF

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.

Fig. 5 Layout of the experimental setup.

Download Full Size | PDF

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 .

Table 1. Experimental results for small to large phase aberrations

View Table

Fig. 6 Experimental results for small to large phase aberrations, where the graphs from top to bottom are the RMS of residual errors and the times of optimizing iterations, respectively.

Download Full Size | PDF

Fig. 7 Experimental results of phase retrieval and deconvolution for the phase aberration of RMS = 0.647λ. (a) and (d) are the measured and reconstructed phases respectively, where the RMS of residual error is 0.046λ. (b) and (e) are the measured and reconstructed focal PSF, respectively. (c) and (f) are the simulated focal image and reconstructed object, respectively.

Download Full Size | PDF

Fig. 8 Experimental results of phase retrieval and deconvolution for the phase aberration of RMS = 1.203λ. (a) and (d) are the measured and reconstructed phases respectively, where the RMS of residual error is 0.202λ. (b) and (e) are the measured and reconstructed focal PSF, respectively. (c) and (f) are the simulated focal image and reconstructed object, respectively.

Download Full Size | PDF

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.

5. Conclusion

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

1. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef] [PubMed]

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]

4. J. R. Fienup, “Phase-retrieval algorithms for a complicated optical system,” Appl. Opt. 32(10), 1737–1746 (1993). [CrossRef] [PubMed]

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]

7. B. M. Hanser, M. G. L. Gustafsson, D. A. Agard, and J. W. Sedat, “Phase-retrieved pupil functions in wide-field fluorescence microscopy,” J. Microsc. 216(1), 32–48 (2004). [CrossRef] [PubMed]

8. A. J. E. M. Janssen, “Extended Nijboer-Zernike approach for the computation of optical point-spread functions,” J. Opt. Soc. Am. A 19(5), 849–857 (2002). [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.

18. C. J. R. Sheppard and P. Török, “Focal shift and the axial optical coordinate for high-aperture systems of finite Fresnel number,” J. Opt. Soc. Am. A 20(11), 2156–2162 (2003). [CrossRef] [PubMed]

19. S. M. Jefferies, M. Lloyd-Hart, E. K. Hege, and J. Georges, “Sensing wave-front amplitude and phase with phase diversity,” Appl. Opt. 41(11), 2095–2102 (2002). [CrossRef] [PubMed]

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]

RMS of measured phases (λ)	RMS of residual errors of reconstructed phases (λ)			Times of optimizing iterations
RMS of measured phases (λ)	average	minimum	maximum	average	minimum	maximum
~0.2	0.030	0.027	0.033	19	15	23
~0.4	0.042	0.036	0.047	35	29	41
~0.6	0.051	0.045	0.055	40	27	52
~0.8	0.075	0.068	0.089	41	36	48
~1.0	0.110	0.089	0.122	55	39	58
~1.2	0.214	0.180	0.251	57	48	66

A three-dimensional point spread function for phase retrieval and deconvolution

Abstract

1. Introduction

2. The three-dimensional PSF

2.1 Three-dimensional focal field distribution

2.2 Formulation of three-dimensional PSF

3. The phase diversity algorithm

3.1 Description of the algorithm

3.2 Implementation and analysis

4. Experimental setup and results

5. Conclusion

References and links

Cited By

Figures (8)

Tables (1)

Equations (37)

Optics Express