OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 14 — Jul. 2, 2012
  • pp: 15392–15405
« Show journal navigation

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

Xinyue Liu, Liang Wang, Jianli Wang, and Haoran Meng  »View Author Affiliations


Optics Express, Vol. 20, Issue 14, pp. 15392-15405 (2012)
http://dx.doi.org/10.1364/OE.20.015392


View Full Text Article

Acrobat PDF (1395 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

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

1. Introduction

In most cases of focus-diverse phase retrieval and deconvolution, such as optical metrology, adaptive optics, and deconvolution microscopy [1

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

6

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]

], 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

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]

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

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]

10

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.

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

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

11. A. A. Ramos and A. L. Ariste, “Image reconstruction with analytical point spread functions,” Astron. Astrophys . 518, paper A6 (2010).

].

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

12. C. W. McCutchen, “Generalized Aperture and the Three-Dimensional Diffraction Image,” J. Opt. Soc. Am. 54(2), 240–242 (1964). [CrossRef]

], 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

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]

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

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]

]. 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
Fig. 1 Geometry for computing the optical field distribution in the focal region.
, 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

15. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press, 1999).

]
U(P)=iAλexp(ikf)fWU(Q)exp(ikr)rcosθdS,
(1)
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

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.

,12

12. C. W. McCutchen, “Generalized Aperture and the Three-Dimensional Diffraction Image,” J. Opt. Soc. Am. 54(2), 240–242 (1964). [CrossRef]

], the focal field distribution can be expressed in the Cartesian coordinates as
UF(x,y,z)=iAλdf3VUW(u,v,w)A(u,v,w)exp(ikxu+yv+zwf)dudvdw,
(2)
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)=δ(u2+v2+w2f),
(3)
where δ(●) is the Kronecker delta function. Like the two-dimensional case, we can define the generalized aperture function as

P(u,v,w)=UW(u,v,w)A(u,v,w).
(4)

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

16. S. van Haver, “The Extended Nijboer-Zernike diffraction theory and its applications,” PhD Dissertation, Delft University of Technology (2010).

,17

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.

] 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
z¯'=z¯/(1+z¯λfπ(NA)2),r¯'=r¯/(1+z¯λfπ(NA)2),
(5)
where z¯ and r¯ are the cylindrical coordinates of focal field normalized with the axial unit λ/π(NA)2 and lateral unit λ/2π(NA) respectively as

z¯=zπ(NA)2λ,r¯=x2+y22π(NA)λ.
(6)

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

x'=ff+zx,y'=ff+zy,z'=ff+zz.
(7)

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

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]

], 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

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]

], which also derives a scaled three-dimensional Fourier transform expression of focal field distribution by using the correction of Eq. (7) as

UF(x',y',z')=iAλfz'f3exp(ikx'2+y'2+2z'22(fz'))×VUW(u,v,w)A(u,v,w)exp(ikx'u+y'v+z'wf)dudvdw.
(8)

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)F3[P(u,v,w)],
(9)
where F3[●] denotes the three-dimensional Fourier transform. Following the treatment of the ENZ approach [10

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.

], we can expand the generalized aperture function P(u,v,w) with the Zernike polynomials
P(u,v,w)=nmβnmZ'nm(u,v,w),
(10)
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'nm(u,v,w)=Znm(u,v)A(u,v,w),
(11)
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)=nmβnmV'nm(x,y,z),
(12)
where V'(x,y,z) are the three-dimensional basis functions for the expansion of focal field distribution, which are defined by

V'nm(x,y,z)=g(x,y,z)F3[Z'nm(u,v,w)].
(13)

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)nm2imβnmVnm(r,z)exp(imφ),
(14)
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
Vnm(r,z)=01exp(izρ2)Rn|m|(ρ)Jm(2πrρ)ρdρ,
(15)
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'nm(r,φ,z)=2imh(r,φ,z)Vnm(r,z)exp(imφ).
(16)

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)=|nmβnmV'nm(x,y,z)|2.
(17)

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

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]

]. 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)S(x,y,z)+N(x,y,z),
(18)
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

20. C. R. Vogel, T. Chan, and R. Plemmons, “Fast algorithms for phase diversity-based blind deconvolution,” Proc. SPIE 3353, 994–1005 (1998). [CrossRef]

] 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

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]

,6

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]

].

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
E3D=12uvw|I˜(u,v,w)O˜(u,v,w)S˜(u,v,w)|2+γ2uvw|O˜(u,v,w)|2+μuvw[|P(u,v,w)|2|P0(u,v,w)|2]2,
(19)
where ~denotes the Fourier spectrum, |P0(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

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]

]. 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
E2D=12uvz|Iz˜(u,v)O˜(u,v)Sz˜(u,v)|2+γ2uv|O˜(u,v)|2+μuvw[|P(u,v,w)|2|P0(u,v,w)|2]2.
(20)
where Sz 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
E'3D=uvw|I˜*(u,v,w)S˜(u,v,w)|2γ+|S˜(u,v,w)|2+uvw|I˜(u,v,w)|2+μuvw[|P(u,v,w)|2|P0(u,v,w)|2]2,
(21)
E'2D=uvz|Iz˜*(u,v)Sz˜(u,v)|2γ+z|Sz˜(u,v)|2+uvz|Iz˜(u,v)|2+μuvw[|P(u,v,w)|2|P0(u,v,w)|2]2,
(22)
and the minimizers are expressed by

O˜3D(u,v,w)=I˜(u,v,w)S˜*(u,v,w)γ+|S˜(u,v,w)|2,
(23)
O˜2D(u,v)=zIz˜(u,v)Sz˜*(u,v)γ+z|Sz˜(u,v)|2.
(24)

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=klβkβl*V'kV'l*=klβkβl*Wkl,Wkl=V'kV'l*,
(25)
where * denotes the complex conjugate and Wkl 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

S˜=klβkβl*Wkl˜,Sz˜=klβkβl*Wkl(z)˜.
(26)

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
Sz˜βkRe=lβl*Wkl(z)˜+lβlWlk(z)˜,Sz˜βkIm=ilβl*Wkl(z)˜ilβlWlk(z)˜,
(27)
2Sz˜βkReβlRe=2Sz˜βkImβlIm=Wlk(z)˜+Wkl(z)˜,2Sz˜βkReβlIm=2Sz˜βkImβlRe=iWlk(z)˜iWkl(z)˜.
(28)
where Re and Im denote the real part and imaginary part, respectively. Using Eq. (10), |P|2 can be represented as
|P|2=klβkβl*Z'kZ'l=klβkβl*Xkl,Xkl=Z'kZ'l,
(29)
where Xkl are the basis functions for the expansion of |P|2. The gradient and Hessian of |P|2 with respect to the real and imaginary parts can be obtained accordingly

|P|2βkRe=2lβlReXkl,|P|2βkIm=2lβlImXkl,
(30)
2|P|2βkReβlRe=2Sz˜βkImβlIm=2Xkl,2Sz˜βkReβlIm=2|P|2βkImβlRe=0.
(31)

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

E'2D=uv|T|2D+uvz|Iz˜|2+μuvwQ2,T=Iz˜*Sz˜,D=γ+z|Sz˜|2,Q=|P|2|P0|2.
(32)

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

|T|2βkRe,Im=2Re[T*zIz˜*Sz˜βkRe,Im],2|T|2βkRe,ImβlRe,Im=2Re[T*zIz˜*2Sz˜βkRe,ImβlRe,Im+zIz˜*Sz˜βkRe,ImzIz˜Sz˜*βlRe,Im],
(33)
DβkRe,Im=2Re[zSz˜*Sz˜βkRe,Im],2DβkRe,ImβlRe,Im=2Re[z(Sz˜*2Sz˜βkRe,ImβlRe,Im+Sz˜βkRe,ImSz˜*βlRe,Im)],
(34)
Q2βkRe,Im=2[Q|P|2βkRe,Im],2Q2βkRe,ImβlRe,Im=2[Q2|P|2βkRe,ImβlRe,Im+|P|2βkRe,Im|P|2βlRe,Im].
(35)

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

E'2DβkRe,Im=uv[|T|2D2DβkRe,Im1D|T|2βkRe,Im+μQ2βkRe,Im],
(36)
2E'2DβkRe,ImβlRe,Im=uv[2|T|2D3DβkRe,ImDβlRe,Im+|T|2D22DβkRe,ImβlRe,Im1D2|T|2βkRe,ImβlRe,Im+1D2DβkRe,Im|T|2βlRe,Im+1D2|T|2βkRe,ImDβlRe,Im+μ2Q2βkRe,ImβlRe,Im].
(37)

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 [21

21. J. Nocedal and S. J. Wright, Numerical Optimization 2nd ed. (Springer, 2006).

] and the Goldstein’s branch cut method for phase unwrapping [22

22. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley-Interscience, 1998).

], where the chirp z-transform is used to perform the three-dimensional Fourier transform [23

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]

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

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]

] 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
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λ.
, 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 α [11

11. A. A. Ramos and A. L. Ariste, “Image reconstruction with analytical point spread functions,” Astron. Astrophys . 518, paper A6 (2010).

]. 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
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.
, 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
Fig. 5 Layout of the experimental setup.
, 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

Table 1. Experimental results for small to large phase aberrations

table-icon
View This Table
and shown in Fig. 6
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.
, 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
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.
-8
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.
.

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]

OCIS Codes
(010.7350) Atmospheric and oceanic optics : Wave-front sensing
(100.1830) Image processing : Deconvolution
(100.5070) Image processing : Phase retrieval
(100.6890) Image processing : Three-dimensional image processing

ToC Category:
Image Processing

History
Original Manuscript: May 9, 2012
Revised Manuscript: June 16, 2012
Manuscript Accepted: June 18, 2012
Published: June 25, 2012

Citation
Xinyue Liu, Liang Wang, Jianli Wang, and Haoran Meng, "A three-dimensional point spread function for phase retrieval and deconvolution," Opt. Express 20, 15392-15405 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-14-15392


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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. A9(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. A24(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. A19(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. A19(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. A20(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. SPIE3353, 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. Express14(23), 11277–11291 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-23-11277 . [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited