OSA's Digital Library

Journal of the Optical Society of America A

Journal of the Optical Society of America A

| OPTICS, IMAGE SCIENCE, AND VISION

  • Editor: Franco Gori
  • Vol. 30, Iss. 10 — Oct. 1, 2013
  • pp: 2002–2011
« Show journal navigation

Iterative linear focal-plane wavefront correction

C. S. Smith, R. Marinică, A. J. den Dekker, M. Verhaegen, V. Korkiakoski, C. U. Keller, and N. Doelman  »View Author Affiliations


JOSA A, Vol. 30, Issue 10, pp. 2002-2011 (2013)
http://dx.doi.org/10.1364/JOSAA.30.002002


View Full Text Article

Acrobat PDF (520 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We propose an efficient approximation to the nonlinear phase diversity (PD) method for wavefront reconstruction and correction from intensity measurements with potential of being used in real-time applications. The new iterative linear phase diversity (ILPD) method assumes that the residual phase aberration is small and makes use of a first-order Taylor expansion of the point spread function (PSF), which allows for arbitrary (large) diversities in order to optimize the phase retrieval. For static disturbances, at each step, the residual phase aberration is estimated based on one defocused image by solving a linear least squares problem, and compensated for with a deformable mirror. Due to the fact that the linear approximation does not have to be updated with each correction step, the computational complexity of the method is reduced to that of a matrix-vector multiplication. The convergence of the ILPD correction steps has been investigated and numerically verified. The comparative study that we make demonstrates the improved performance in computational time with no decrease in accuracy with respect to existing methods that also linearize the PSF.

© 2013 Optical Society of America

1. INTRODUCTION

All optical measurements are subject to optical aberrations either coming from exterior sources or intrinsic to the instrument. If the aberrations can be estimated, they can be compensated for either through adaptive optics during image acquisition or postprocessing. One method that has been used mostly in postprocessing is phase diversity (PD) [1

1. R. G. Paxman, T. J. Schultz, and J. R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” J. Opt. Soc. Am. A 9, 1072–1085 (1992). [CrossRef]

]. PD estimates wavefront aberrations using nonlinear optimization techniques from multiple images of the same unknown scene acquired simultaneously, which contain additional user-introduced aberrations, the latter referred to as diversities. To be able to uniquely estimate wavefront aberrations, more than one in-focus image is needed [2

2. L. M. Mugnier, A. Blanc, and J. Idier, “Phase diversity: a technique for wave-front sensing and for diffraction-limited imaging,” Adv. Imaging Electron Phys. 141, 1–76 (2006). [CrossRef]

], because rotating a wavefront by 180° and flipping its sign produces the same point spread function (PSF) as the original wavefront [3

3. D. R. Gerwe, M. M. Johnson, and B. Calef, “Local minima analysis of phase diverse phase retrieval using maximum likelihood,” The Advanced Maui Optical and Space Surveillance Technical Conference, Kihei, Hawaii (2008).

]. The resulting optimization problem is nonlinear and is known to be computationally complex due to the repetitive evaluations of Fourier transforms. In addition, the method is also prone to converge to local optima [3

3. D. R. Gerwe, M. M. Johnson, and B. Calef, “Local minima analysis of phase diverse phase retrieval using maximum likelihood,” The Advanced Maui Optical and Space Surveillance Technical Conference, Kihei, Hawaii (2008).

]. As a consequence, nonlinear PD has a limited usage in real-time correction algorithms [2

2. L. M. Mugnier, A. Blanc, and J. Idier, “Phase diversity: a technique for wave-front sensing and for diffraction-limited imaging,” Adv. Imaging Electron Phys. 141, 1–76 (2006). [CrossRef]

], and different ideas have been presented to decrease the complexity of the calculations. These ideas can be split up into Fourier domain [4

4. R. W. Gerchberg and W. O. Saxton, “Phase retrieval by iterated projections,” Optik 35, 237–246 (1972).

7

7. F. Martinache, “Kernel phase in Fizeau interferometry,” Astrophys. J. 724, 464–469 (2010). [CrossRef]

] and spatial domain [8

8. S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]

,9

9. C. van der Avoort, J. J. M. Braat, P. Dirksen, and A. J. E. M. Janssen, “Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer–Zernike approach,” J. Mod. Opt. 52, 1695–1728 (2005). [CrossRef]

] techniques. The Gerchberg–Saxton [4

4. R. W. Gerchberg and W. O. Saxton, “Phase retrieval by iterated projections,” Optik 35, 237–246 (1972).

] algorithm is one of the oldest and best known Fourier domain techniques, which is an iterative algorithm for retrieving the phase from intensity measurements. Spatial domain techniques make use of a local model for the PSF, but do not use the Fourier transform. The common idea in decreasing the computational complexity is the approximation of the PSF based on the assumption that the total aberration is small [5

5. R. A. Gonsalves, “Small-phase solution to the phase-retrieval problem,” Opt. Lett. 26, 684–685 (2001). [CrossRef]

,6

6. C. U. Keller, V. Korkiakoski, N. Doelman, R. Fraanje, R. Andrei, and M. Verhaegen, “Extremely fast focal-plane wavefront sensing for extreme adaptive optics,” Proc. SPIE 8447, 844721 (2012). [CrossRef]

,10

10. I. Mocœur, L. M. Mugnier, and F. Cassaing, “Analytical solution to the phase-diversity problem for real-time wavefront sensing,” Opt. Lett. 34, 3487–3489 (2009). [CrossRef]

]. This small-phase assumption is associated in the literature with the Born approximation [5

5. R. A. Gonsalves, “Small-phase solution to the phase-retrieval problem,” Opt. Lett. 26, 684–685 (2001). [CrossRef]

,11

11. W. J. Wild, “Linear phase retrieval for wave-front sensing,” Opt. Lett. 23, 573–575 (1998). [CrossRef]

,12

12. J. Antonello, M. Verhaegen, R. Fraanje, T. van Werkhoven, H. C. Gerritsen, and C. U. Keller, “Semidefinite programming for model-based sensorless adaptive optics,” J. Opt. Soc. Am. A 29, 2428–2438 (2012). [CrossRef]

], which implicitly assumes that the diversity used is small.

Recently, in [6

6. C. U. Keller, V. Korkiakoski, N. Doelman, R. Fraanje, R. Andrei, and M. Verhaegen, “Extremely fast focal-plane wavefront sensing for extreme adaptive optics,” Proc. SPIE 8447, 844721 (2012). [CrossRef]

], it was shown that using a second-order expansion of the generalized pupil function (GPF), wavefront retrieval algorithms give more accurate results than using the Born approximation, which results from a linear expansion of the GPF. The key assumption of these methods is that the sum of the diversity and the aberration is small. However, as has been shown in [13

13. D. J. Lee, M. C. Roggemann, and B. M. Welsh, “Cramer-rao analysis of phase-diverse wave-front sensing,” J. Opt. Soc. Am. A 16, 1005–1015 (1999). [CrossRef]

], the optimal diversity depends on the present aberration and can generally not be considered small. In the present paper, we overcome this shortcoming by the use of an alternative approximation of the PSF. The linearization of the PSF is done around zero aberration and a (possibly large) diversity and it is suited for small values (0.5 radians (rad) root mean square (rms) [11

11. W. J. Wild, “Linear phase retrieval for wave-front sensing,” Opt. Lett. 23, 573–575 (1998). [CrossRef]

,12

12. J. Antonello, M. Verhaegen, R. Fraanje, T. van Werkhoven, H. C. Gerritsen, and C. U. Keller, “Semidefinite programming for model-based sensorless adaptive optics,” J. Opt. Soc. Am. A 29, 2428–2438 (2012). [CrossRef]

]) of the phase aberration. The iterative manner in which the method is applied compensates for this small-phase assumption. In this context, the use of one image is enough for the uniqueness of the phase estimate [14

14. A. Tokovinin and S. Heathcote, “Donut: measuring optical aberrations from a single extrafocal image,” Publ. Astron. Soc. Pac. 118, 1165–1175 (2006). [CrossRef]

].

The paper is organized as follows. In Section 2, we present the general problem and introduce the PSF of the optical system and the noise model. In Section 3, we review four linear and quadratic PSF approximations and show the advantages and disadvantages of each of them, which we prove in Appendix A. In Section 4, we use the previously mentioned approximations and present the ILPD solution. In Section 5, we discuss results of numerical simulations and compare them to the ones in [8

8. S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]

]. We end with conclusions in Section 6.

Some mathematical notations used are standard: T and * denote transposition and transpose conjugation, respectively, denotes the convolution operator, denotes the vector 2-norm, O() describes the complexity of a function when the argument tends toward a particular value, usually in terms of simpler functions, O(a) is the a-th order Lagrange residue, R and C are the sets of real and complex numbers, respectively, Rm×n and Cm×n are the sets of m×n matrices with elements in the set of real or complex numbers, respectively.

2. OPTICAL SYSTEM

In this section, a model is presented for the image formation of a point source in the presence of phase aberrations ϕRm2×1, approximated using a normalized Zernike basis [16

16. ANSI, “Methods for reporting optical aberrations of eyes,” American National Standard for Ophtalmics ANSI-Z80.28-2004 (2004).

]
ϕ(uj,vj)=Z(uj,vj)Tα,
(1)
where αRn×1 contains the Zernike coefficients corresponding to the unknown aberration and ZRn×m2 is a matrix containing the n Zernike polynomials evaluated in the pupil plane coordinates (uj,vj). Besides the “in-focus” image, PD uses additional images with known diversities. The phase aberration in the i-th diversity image is
ϕi(uj,vj)=Z(uj,vj)T(α+βi),
(2)
where βiRn×1 is a known diversity. These phase aberrations nonlinearly influence the PSF. The incoherent image formation of a point source is given by [17

17. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005).

]
yi,j=μih(sj,tj,α,βi)+ni(sj,tj),
(3)
where yi,j denotes the j-th pixel of the i-th diversity image, μi is the number of photons (the expected arrival rate multiplied with the integration time of the camera), h denotes the spatially invariant PSF expressed in the spatial coordinates (sj,tj) with aberration α and user-introduced diversity βi, and ni(sj,tj) is Gaussian white noise with standard deviation σi,j, which we assume to be equal for all pixels by dropping the index j. If only a defocus aberration is present, the schematic representation of adding a defocus diversity is given in Fig. 1.

Fig. 1. Optical system: focal plane (black), defocus plane (gray), unknown small aberration Zernike coefficients (α), known arbitrary diversity Zernike coefficients (βi).

In Subsection 2.A the aberrated PSF presented in Eq. (3) is derived. Subsequently, in Subsection 2.B the measurement noise is presented.

A. Image Formation

The spatially invariant PSF of the i-th optical path in Eq. (3) is given by [17

17. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005).

]
h(sj,tj,α,βi)=F(Π(u,v)exp(iϕi(u,v)))(sj,tj)×F(Π(u,v)exp(iϕi(u,v)))*(sj,tj)=F((exp(iϕi(u,v))Π(u,v)exp(iϕi(u,v))Π(u,v))(u,v))(sj,tj),
(4)
where F(·) is the Fourier transform, (uj,vj)=(2πsj/(fλ),2πtj/(fλ)), f is the focal length, λ is the wavelength, ϕi is the phase, and Π is the pupil function. Next, we define the GPF as
p(uj,vj,α,βi)=Π(uj,vj)exp(iϕi(uj,vj)).
(5)
Using Eq. (5), the optical transfer function (OTF) is given by
W(uj,vj,α,βi)=(p(u,v,α,βi)p(u,v,α,βi)*)(uj,vj).
(6)

Next, we introduce the short-hand notations
pj(α,βi)p(uj,vj,α+βi),pj(α,βi)p(uj,vj,α+βi),
(7)
and
Wj(α,βi)W(uj,vj,α,βi),hj(α,βi)h(sj,tj,α,βi).
(8)
Using Eqs. (7), (8), and (4), the PSF becomes
hj(α,βi)=F(pj(α,βi)pj*(α,βi))(sj,tj)=F(Wj(α,βi))(sj,tj),
(9)
where, for simplicity, we drop the convolution brackets and its new coordinates and only use the star operator.

B. Measurement Noise

We consider here two main noise processes that are dependent on the exposure time and luminosity of the object, namely the Gaussian read-out noise and the photon counting noise. We approximate the photon noise by an additive zero-mean Gaussian noise, with a variance equal to the flux. The read-out noise is the same for each pixel and follows a Gaussian distribution, which has the property that all the pixels are mutually uncorrelated. The SNR level is given by
SNR=1m2j=1m2μihj(α,βi)1m2j=1m2σi2(=μim2σi),
(10)
where μi and hj(α,βi) follow from Eq. (3) and σi is the standard deviation of the read-out noise in the i-th diversity image. The total noise has zero mean and its variance is the sum of the two variances. The approximation of the photon noise is made only in the theoretical part of the paper. For the numerical simulations, photon noise is modeled using a Poisson distribution.

3. APPROXIMATIONS OF THE PSF

A. First-Order Approximations

The approximations presented here are all based on a linear Taylor expansion of the GPF or of the PSF, respectively. In Subsection 3.A.1, the assumption is that both the wavefront phase and the diversity used are small. We approximate the GPF with a linear expression and compute the coefficients of the resulting quadratic PSF. The approximation given in Subsection 3.A.2 can be used for small wavefronts and a general diversity. This is simply the Taylor expansion of the PSF. It approximates the PSF around the diversity with a linear expression.

1. Small Total Phase Approximation

The Born approximation assumes a small phase, ϕi=ZT(α+βi), such that the GPF can be approximated using only a first-order Taylor expansion around α+βi=0. The consequence is that the GPF can be written as
pj(α,βi)=pj(α,βi)|α+βi=0+pj(α,βi)(α+βi)|α+βi=0(α+βi)+O(α+βi2),=Π(uj,vj)(1+iZT(uj,vj)(α+βi))+O(α+βi2).
(11)

Substituting this first-order (Born) approximation of the GPF in Eq. (4) yields a quadratic PSF
hj(α,βi)=A0,j+A1,j(α+βi)+(α+βi)TA2,j(α+βi)+O(α+βi2),
(12)
where
A0,jF(pj(α,βi)pj*(α,βi))(sj,tj)|α+βi=0=hj(α,βi)|α+βi=0,A1,jF(pj(α,βi)(α+βi)pj*(α,βi)+pj(α,βi)pj*(α,βi)(α+βi)T)(sj,tj)|α+βi=0=hj(α,βi)(α+βi)|α+βi=0,A2,jF(pj(α,βi)(α+βi)pj*(α,βi)(α+βi)T+pj(α,βi)(α+βi)pj*(α,βi)(α+βi)T)(sj,tj)|α+βi=0.
(13)

Property 3.1. The linear term of the approximated PSF in Eq. (12) is invariant in the even aberrations.

2. Small Aberration Approximation

Another first-order model is obtained by directly approximating the PSF in Eq. (4) for small aberrations and nonzero diversities. The first-order Taylor approximation of the PSF in α=0 is given by
hj(α,βi)=B0,j(βi)+B1,j(βi)α+O(α2),
(14)
where
B0,j(βi)hj(α,βi)|α=0,B1,j(βi)hj(α,βi)α|α=0.
(15)

Property 3.2. The linear terms of the first-order approximation of the PSF and of the PSF resulting from the first-order Taylor approximation of the GPF are equal.

Property 3.3. The approximation in Eq. (14) has the property that for ϕ0 the even modes do not cancel in the linear term.

Remark 3.1. Note that this approximation is valid for any diversity. Due to the fact that the linear term is not invariant in the even modes, we can estimate the even and odd modes with just a linear equation, as will be shown in a later section.

B. Second-Order Approximations

1. Small Total Phase Approximation

It has been shown in [6

6. C. U. Keller, V. Korkiakoski, N. Doelman, R. Fraanje, R. Andrei, and M. Verhaegen, “Extremely fast focal-plane wavefront sensing for extreme adaptive optics,” Proc. SPIE 8447, 844721 (2012). [CrossRef]

] that an additional quadratic term leads to a more accurate PSF approximation than using the Born approximation. This term is obtained using a second-order Taylor expansion of the GPF in ϕ=0 and neglecting the third and the fourth orders of the resulting PSF. The resulting approximation is given by
hj(α,βi)=C0,j+C1,j(α+βi)+(α+βi)TC2,j(α+βi)+O(α+βi3),
(16)
where
C0,jhj(α,βi)|α+βi=0(=A0,j),C1,jhj(α,βi)(α+βi)|α+βi=0(=A1,j),C2,j2hj(α,βi)(α+βi)(α+βi)T|α+βi=0=A2,j+F(2pj(α,βi)(α+βi)(α+βi)Tpj*(α,βi)+pj(α,βi)2pj*(α,βi)(α+βi)(α+βi)T)(sj,tj)|α+βi=0.
(17)

Property 3.4. The expression in Eq. (16) is also obtained when the PSF is approximated using a second-order Taylor expansion around ϕ=0.

2. Small Aberration Approximation

The second-order Taylor approximation of the PSF is given by
hj(α,βi)=D0,j(βi)+D1,j(βi)α+αTD2,j(βi)α+O(α3),
(18)
where
D0,j(βi)hj(α,βi)|α=0(=B0,j(βi)),D1,j(βi)hj(α,βi)α|α=0(=B1,j(βi)),D2,j(βi)2hj(α,βi)ααT|α=0.

Property 3.5. The second-order Taylor approximation of the PSF in ϕ=0 is more accurate than the PSF obtained from the first-order GPF approximation in ϕ=0, while the quadratic form remains.

Property 3.6. The second-order Taylor approximation of the PSF in ϕ0 has the property that the even modes do not cancel in the linear term of the PSF.

4. ITERATIVE LINEAR PHASE DIVERSITY

Apart from Eq. (14), all the other approximations of the PSF derived in Section 3 given by Eqs. (12), (16), and (18) are quadratic in the unknown aberration, as represented by α. Here, we aim to obtain a linear relationship between the measured intensity and the unknown aberration due to the fact that a linear model has low computational complexity and gives rise to fast algorithms. The approximations in Eqs. (12) and (16) are based on the Taylor series expansion of the GPF, first order and second order, respectively, and only the terms of the PSF up to the second order are retained. This means that Eq. (16) is more accurate than Eq. (12), which motivates its preferred use. In order to obtain a linear formulation using the approximation in Eq. (16), we could take the difference of two measurements as done in [12

12. J. Antonello, M. Verhaegen, R. Fraanje, T. van Werkhoven, H. C. Gerritsen, and C. U. Keller, “Semidefinite programming for model-based sensorless adaptive optics,” J. Opt. Soc. Am. A 29, 2428–2438 (2012). [CrossRef]

]. Note that this artifice cannot be performed on Eq. (18), because the coefficients D0,j, D1,j, and D2,j are functions of the diversities βi and they do not cancel when two measurements are subtracted, such that the only solution when using this type of approximation is to retain only the linear term as in Eq. (14). Nevertheless, subtracting two measurements, first of all implies obtaining two measurements, which we want to avoid in this work, and second, the numerical estimation problem is not well conditioned when noise is also subtracted. Also, the SNR decreases when taking differences and the following property states this.

Property 4.1. Taking the difference between two images significantly decreases the SNR.

In what follows, we form a linear system using the approximation in Eq. (14). For this, we use one defocused image. The solution of ILPD using the previously mentioned approximation is compared with the solution given in [8

8. S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]

]. The approximation in Eq. (14) is already linear and we have as follows
Y1=bS+ASα+ΔbS(α)+n1,
(19)
ΔbS(α)O(α2),
(20)
where Yi[yi,1yi,j]T, bS[b˜S,1b˜S,i]T, with b˜S,i[B0,i,1B0,i,j], AS[A˜S,1TA˜S,iT]T, with A˜S,i[Bi,1(βi)TBi,j(βi)T]T, and B0,j(βi), B1,j(βi) are defined in Eq. (15). The index i is kept when defining the quantities above to suggest that this method can easily be generalized to more than one image if the optical system can facilitate this, while still keeping the linear system formulation.

The main advantages of the first-order Taylor approximation of the PSF in Eq. (14) are that it is possible to approximate the PSF at large diversities and that the first-order term is not invariant in the even modes, which makes it possible to estimate them (except for piston). Therefore, we do not have to subtract images, which significantly decreases the SNR.

In ILPD, the residual aberration is repeatedly estimated and compensated for with a DM using the residual aberration estimate. Assuming that the DM can fully compensate for the estimated residual aberration, then, denoting the residual aberration estimate at the k-th correction step by α^k1, and denoting the residual aberration at the k-th correction step by αk1, we obtain
αk=Δαk1=αk1α^k1.
(21)
At the k-th correction step, one image, Y1,k, is recorded with diversity β1 assuming the wavefront aberration does not change. The additional index k of Y denotes the correction step. From the new image, a new estimate of αk is obtained via the solution of a LS problem based on Eq. (19). The algorithm continues until the strength of the aberration decreases under a certain threshold or a finite number of correction steps has been performed. Let Eq. (19) (where the step index k has been added) be rewritten as
b¯S,kΔbS(αk)=ASαk+nk,
(22)
where b¯S,kY1,kbS, nk=n1,k, and
ΔbS(αk)=O(αk2)=CSαk2,
with CS a constant defined by the Lagrange remainder. Then, the LS problem that needs to be solved is
minαkb¯S,kASαk2.
(23)
The solution of Eq. (23) after each correction step k (no correction for at the zeroth correction step) with the DM will be indicated by the ILPD method for joint wavefront estimation and correction.

A. Convergence Analysis

In this section, we study the convergence behavior of the ILPD method in the absence of measurement noise. Using Eq. (21), the relative residue after correction step k is given by
rLSαk+1αk.
(24)
The relative residue has to be smaller than one to ensure that the rms value of the wavefront is reduced. If this is not the case, the rms value increases or remains constant. Therefore, the convergence can be quantified using this quantity. We validate by Monte Carlo simulation that using the PSF approximation proposed in Eq. (14) we converge to an unbiased estimate faster than the method in [8

8. S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]

].

For the linear system in Eq. (22), in the noiseless case, we can compute an approximate upper bound on the relative error in the solution [19

19. G. H. Golub and C. F. Van Loan, Matrix Computations (Johns Hopkins University, 1996).

] as
α^kαkαkΔbS(αk)b¯S,k{2κ(AS)cos(θ)+tan(θ)κ(AS)2},
(25)
where κ(AS) denotes the condition number of AS and θ is the acute angle between the vectors ASα^k and b¯S,k. For a well-conditioned matrix, the bound depends on ΔbS(αk). As αk decreases, ΔbS(αk) decreases and the bound becomes zero in the limit
limαk0CSαk2b¯S,k{2κ(AS)cos(θ)+tan(θ)κ(AS)2}=0.
(26)
Using the approximation in Eq. (14) to formulate our problem, it is clear that the model error in Eq. (23) only depends on the unknown aberration. This would be different when differences of two PSFs modeled by the approximation in Eq. (16) are taken, when the model error also depends on the chosen diversity. Then, a compromise should be made between a small diversity which leads to a small model error and a large diversity which ensures that the difference between two images does not become zero and the information content is lost.

5. SIMULATIONS

In this section, we present numerical simulations for the iterative aberration correction problem using ILPD and LIFT. We first describe the simulation setup. Second, we give one example of ILPD. Next, we analyze the behavior of both methods using a Monte Carlo simulation by varying the noise level and the rms value of the initial aberration. The computer employed for these simulations is a 2.67 GHz quad-core Intel Core (TM)2 Quad CPU Q8400 with 4.0 GB of RAM.

We compute the corresponding nonlinear PSFs in Eq. (4) and the coefficients of the linear Taylor expansion given in Eq. (15). We consider a pupil of radius r sampled on a 32×32 grid embedded in a 4r×4r image to satisfy the Nyquist sampling criterion. We ensure that the wavefront does not contain jumps larger than π/2, which would be problematic for the sampling process. The SNR corresponding to the read-out noise is calculated over the image (m=32). Also, all treatment is monochromatic. We assume that the DM is able to produce known diversity shapes with an error that is negligible compared to other sources. This assumption motivates our choice to model only the first n=14 modes that can be corrected by the DM. To obtain a PSF of unit surface, the pupil function Π is chosen as
Π(uj,vj)={1Suj2+vj2r20uj2+vj2>r2,
(27)
where S is the physical surface of the pupil.

If the frame rates of the imaging camera and of the DM are sufficiently fast, it is an acceptable approximation that a few sequential wavefronts are assumed to be identical. The static aberration consisting of n modes is generated using the Kolmogorov model [20

20. R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. A 66, 207–211 (1976). [CrossRef]

], which assumes aberrations with zero mean and covariance matrix Cϕ. The parameters used to generate the Kolmogorov model are: diameter D=1[m], outer scale L0=42[m], Fried parameter r0=0.3[m].

First, we give an example of ILPD in Subsection 5.A. Second, in Subsection 5.B, for the same aberration, in the noiseless case, we show the convergence and the corresponding rate of convergence in terms of residual wavefront error and relative residual wavefront error. Subsequently, in Subsection 5.C, we study the convergence properties in terms of the residual error for ILPD and LIFT as a function of increasing read-out noise SNR, photon count, and wavefront rms.

A. One Example of Iterative Phase Diversity

Fig. 2. Convergence in terms of wavefront error: 1 rad rms, no read-out/photon noise, 1000 photons per image. LIFT (top), ILPD (bottom).
Fig. 3. Convergence in terms of wavefront error: 1 rad rms, read-out noise with SNR=3.16, no photon noise, 1000 photons per image. LIFT (top), ILPD (bottom).

Table 1. Rms Values of the Corrected Wavefronts for No Read-Out/Photon Noise, 1000 Photons per Image, and 1 rad Initial rms

table-icon
View This Table
| View All Tables

Table 2. Rms Values of the Corrected Wavefronts for Read-Out Noise SNR=3.16, No Photon Noise, 1000 Photons per Image, and 1 rad Initial rms

table-icon
View This Table
| View All Tables

B. Iterative Linear Phase Diversity Without Noise

In Fig. 4 we plot the residual error in the aberration vector αkα^k at each iteration versus the number of iterations. The residual error decreases with each step for almost all the samples of LIFT and for all the samples considered for ILPD at this particular rms value. In the noiseless case, ILPD converges to a residual error 0. LIFT converges to a small bias different from zero, but not as fast. There are also cases when LIFT diverges.

Fig. 4. Residual error in the aberration vector. On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points not considered outliers, and outliers are plotted individually. The diamond signs represent mean values.

Fig. 5. Relative residue.

The computational time necessary for LIFT to complete five iterations is 10.9978 s on average, while ILPD performs them in 0.0028 s on average. Note that the integration time of the CCD is not included in computing these times. This makes ILPD 3927.8 times faster. When we also count the CCD integration time (of approximately 47.3 ms), ILPD is 40 times faster. For a fair comparison with LIFT, we have used here all the pixels in order to compute the estimate at each step, but the computational time for ILPD further decreases when using just a fraction of the pixels.

C. Error Residue in the Presence of Noise

Fig. 6. Wavefront residual error versus increasing SNR.
Fig. 7. Wavefront residual error versus increasing photon count.
Fig. 8. Wavefront residual error versus increasing rms.

Figure 6 plots the residual error after five correction steps versus increasing read-out noise SNR, considering 1000 photons per image and Poisson photon noise. For each SNR, we repeat the experiment 128 times. The considered initial aberration has 1 rad rms. For ILPD, the residual error decreases with the increase of the SNR, which is what we expected. For LIFT this behavior is not very visible. One reason is the high value of the rms. In our simulations, we have noticed that for smaller rms values, e.g., 0.5 rms, LIFT starts to show this decrease in bias for increasing rms, which shows that LIFT is more appropriate for small rms values. ILPD has a lower error variance and a lower error mean. Clearly, for ILPD, the error variance would converge to zero for an SNR equal to if no Poisson photon noise were considered.

Figure 7 plots the residual error after five correction steps versus increasing photon count. For each photon count, we repeat the experiment 128 times. The considered initial aberration has 1 rad rms. Besides the Poisson photon noise, we also added read-out noise with SNR=3.16. The residual error decreases with the increase of the photon count, which is what we expected. ILPD has a lower error variance and a lower error mean. Furthermore, it is visible in Fig. 7 that at low photon counts LIFT diverges.

The same type of analysis is made in Fig. 8 for increasing rms of the initial wavefront, a constant read-out noise level of 3.16 and Poisson photon noise. Both methods are based on a small-aberration assumption, so the bias of the estimation increases with increasing rms or it takes more iterations to converge. It is visible that LIFT starts to diverge for rms values larger than 0.5 rad, while ILPD still corrects for the aberration. This is due to the fact that with each iteration the aberration becomes smaller and the linear model in Eq. (19) is more and more accurate.

6. CONCLUSIONS

We have presented a novel method for wavefront estimation and correction suitable for several applications in astronomy or microscopy. Under the assumption of small-phase aberrations, which is the typical situation in a control loop, the PSF of an incoherent imaging system has been approximated with a linear model, which can be precomputed if the diversity used is a fixed one. This model allows for arbitrary phase diversities to be introduced in the system. Our iterative approach uses at each step one image of a point-like object, which includes a known phase diversity, and estimates the aberration using a LS approach. In this way we increase the computational speed of phase retrieval methods that linearize the PSF at each iteration around the current estimate of the aberration. Also, as the residual aberration decreases, the precomputed model of the PSF becomes a better fit to the real one. This creates the premises for the method to be applied in a real-time correction system.

APPENDIX A: PROOF OF PROPOSITIONS

In this appendix, the claims made in Sections 3 and 4 are proved. The claims are invariant of the Fourier transform between the PSF and OTF; therefore, to shorten the proofs, the approximations and their properties will be given in terms of the OTF. For clarity we introduce the following short-hand notations
pi,jpj(α,βi),Wi,jWj(α,βi).
(A1)

Proof of Property 3.1. We introduce the short-hand notation γi=α+βi and for the small total phase approximation we have γ0=0. Using Eqs. (12) and (13), the OTF is given by
Wi,j(pi,j+pi,jγiT(γiγ0))(pi,j*+pi,j*γiT(γiγ0))|γi=γ0=pi,jpi,j*|γi=γ0+[pi,jγiTpi,j*+pj(α,βi)pi,j*γiT]|γi=γ0(γiγ0)+(γiγ0)T[pi,jγiTpi,j*γi+pi,jγipi,j*γiT]|γi=γ0(γiγ0).
(A2)
The first-order term is
L[pi,jγiTpi,j*+pj(α,βi)pi,j*γiT]|γi=γ0.
(A3)
To show that Eq. (A3) is invariant in the even modes we reorder γi and ZT by even and odd parts and Eq. (A3) becomes
L=iZjTΠjΠj*ΠjiZjTΠj*=i[Ze,jTZo,jT]ΠjΠj*Πji[Ze,jTZo,jT]Πj*,
(A4)
where the subindexes j and j are short-hand notations for the coordinates (uj,vj) and (uj,vj). Next, because Π is even and real, we have that
L=i[Ze,jTZo,jT]ΠjΠjΠji[Ze,jTZo,jT]Πj=i[02Zo,jT]ΠjΠj.
(A5)

Proof of Property 3.2. The OTF is given by
Wi,jpi,jpi,j*|γi=γ0+[pi,jγiTpi,j*+pj(α,βi)pi,j*γiT]|γi=γ0(γiγ0)
(A6)
and the linear term is equal to Eq. (A3).□

Proof of Property 3.3. We show that Eq. (A3) for γ00 is not invariant for even modes. We again reorder γi and ZT by even modes ZeT and odd modes ZoT to obtain
L=iZjTΠjexp(iZjTγ0)Πj*exp(iZjTγ0)Πjexp(iZjTγ0)iZjTΠjexp(iZjTγ0)=i[Ze,jTZo,jT]Πjexp(iZjTγ0)Πj*exp(iZjTγ0)Πjexp(iZjTγ0)i[Ze,jTZo,jT]Πjexp(iZjTγ0)=i[Ze,jTZo,jT]Π˜jΠ˜j*Π˜ji[Ze,jTZo,jT]Π˜j*,
(A7)
where Π˜jΠjexp(iZjTγ0). Next, because Π˜j is neither even nor real we have that the two terms are different and the even modes do not cancel.□

Proof of Property 3.4. The second-order Taylor approximation of the GPF is
pi,jpi,j|γi=γ0+pi,jγi|γi=γ0(γiγ0)+12(γiγ0)T2pi,jγiγiT|γi=γ0(γiγ0).
(A8)
Dropping terms of order 3 and higher, the resulting approximated OTF reduces to
Wi,jpi,jpi,j*|γi=γ0+[pi,jγiTpi,j*+pj(α,βi)pi,j*γiT]|γi=γ0(γiγ0)+(γiγ0)T[2pi,jγiγiTpi,j*+pj(α,βi)2pi,j*γiγiT+pi,jγiTpi,j*γi+pi,jγipi,j*γiT]|γi=γ0(γiγ0),
(A9)
which is exactly the second-order Taylor approximation.□

Proof of Property 3.5. The difference between the approximated PSF following from a first-order approximation of the GPF in Eq. (A2) and the second-order Taylor approximation in Eq. (A9) is given by
2pi,jγiγiTpi,j*+pj(α,βi)2pi,j*γiγiT.
(A10)
The addition of the missing term from the first-order GPF approximation results in a residue of order O(α3) instead of O(α2). Therefore, the second-order Taylor expansion of the PSF is more accurate than the first-order approximation of the GPF, while the quadratic form remains.□

Proof of Property 3.6. Inspecting Eq. (A10), we observe that the linear term is not affected; therefore, Property 3.1 and Property 3.3 still hold for the linear terms of Eq. (A9).□

Proof of Property 4.1. The intensity of both signals is positive and subtracting two images decreases the mean signal at each pixel. Recall that we assume that all the camera pixels are mutually independent and that the measurement noise is Gaussian distributed
μdiff=E[Δyj]=μ1h1,jμ2h2,j.
(A11)
In Eq. (A11), μdiff is smaller than either μ1h1,j and μ2h2,j. Next, the variance of the signal increases
σdiff2=E[(Δyjμdiff)2]=σ1,j2+σ2,j2.
(A12)
Therefore, the resulting SNR at pixel j is given by
SNRdiff=μ1h1,jμ2h2,jσ1,j2+σ2,j2.
(A13)
If we assume that the noise is the same for each pixel with σi,jσi, the total SNR is equal to
SNRdiff=1m2j=1m2μ1h1,jμ2h2,j1m2j=1m2σ1,j2+1m2jm2σ2,j2=1m21σ12+σ22j=1m2μ1(h1,jμ2h2,j)=1m2μ1μ2σ12+σ22.
(A14)

REFERENCES

1.

R. G. Paxman, T. J. Schultz, and J. R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” J. Opt. Soc. Am. A 9, 1072–1085 (1992). [CrossRef]

2.

L. M. Mugnier, A. Blanc, and J. Idier, “Phase diversity: a technique for wave-front sensing and for diffraction-limited imaging,” Adv. Imaging Electron Phys. 141, 1–76 (2006). [CrossRef]

3.

D. R. Gerwe, M. M. Johnson, and B. Calef, “Local minima analysis of phase diverse phase retrieval using maximum likelihood,” The Advanced Maui Optical and Space Surveillance Technical Conference, Kihei, Hawaii (2008).

4.

R. W. Gerchberg and W. O. Saxton, “Phase retrieval by iterated projections,” Optik 35, 237–246 (1972).

5.

R. A. Gonsalves, “Small-phase solution to the phase-retrieval problem,” Opt. Lett. 26, 684–685 (2001). [CrossRef]

6.

C. U. Keller, V. Korkiakoski, N. Doelman, R. Fraanje, R. Andrei, and M. Verhaegen, “Extremely fast focal-plane wavefront sensing for extreme adaptive optics,” Proc. SPIE 8447, 844721 (2012). [CrossRef]

7.

F. Martinache, “Kernel phase in Fizeau interferometry,” Astrophys. J. 724, 464–469 (2010). [CrossRef]

8.

S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]

9.

C. van der Avoort, J. J. M. Braat, P. Dirksen, and A. J. E. M. Janssen, “Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer–Zernike approach,” J. Mod. Opt. 52, 1695–1728 (2005). [CrossRef]

10.

I. Mocœur, L. M. Mugnier, and F. Cassaing, “Analytical solution to the phase-diversity problem for real-time wavefront sensing,” Opt. Lett. 34, 3487–3489 (2009). [CrossRef]

11.

W. J. Wild, “Linear phase retrieval for wave-front sensing,” Opt. Lett. 23, 573–575 (1998). [CrossRef]

12.

J. Antonello, M. Verhaegen, R. Fraanje, T. van Werkhoven, H. C. Gerritsen, and C. U. Keller, “Semidefinite programming for model-based sensorless adaptive optics,” J. Opt. Soc. Am. A 29, 2428–2438 (2012). [CrossRef]

13.

D. J. Lee, M. C. Roggemann, and B. M. Welsh, “Cramer-rao analysis of phase-diverse wave-front sensing,” J. Opt. Soc. Am. A 16, 1005–1015 (1999). [CrossRef]

14.

A. Tokovinin and S. Heathcote, “Donut: measuring optical aberrations from a single extrafocal image,” Publ. Astron. Soc. Pac. 118, 1165–1175 (2006). [CrossRef]

15.

A. Polo, S. F. Pereira, and P. H. Urbach, “Theoretical analysis for best defocus measurement plane for robust phase retrieval,” Opt. Lett. 38, 812–814 (2013). [CrossRef]

16.

ANSI, “Methods for reporting optical aberrations of eyes,” American National Standard for Ophtalmics ANSI-Z80.28-2004 (2004).

17.

J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005).

18.

F. Gustafsson, Statistical Sensor Fusion (Holmbergs i Malmö AB, 2010).

19.

G. H. Golub and C. F. Van Loan, Matrix Computations (Johns Hopkins University, 1996).

20.

R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. A 66, 207–211 (1976). [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(010.7350) Atmospheric and oceanic optics : Wave-front sensing
(100.5070) Image processing : Phase retrieval
(110.1080) Imaging systems : Active or adaptive optics

ToC Category:
Atmospheric and Oceanic Optics

History
Original Manuscript: May 16, 2013
Revised Manuscript: August 19, 2013
Manuscript Accepted: August 21, 2013
Published: September 13, 2013

Citation
C. S. Smith, R. Marinică, A. J. den Dekker, M. Verhaegen, V. Korkiakoski, C. U. Keller, and N. Doelman, "Iterative linear focal-plane wavefront correction," J. Opt. Soc. Am. A 30, 2002-2011 (2013)
http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-30-10-2002


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. R. G. Paxman, T. J. Schultz, and J. R. Fienup, “Joint estimation of object and aberrations by using phase diversity,” J. Opt. Soc. Am. A 9, 1072–1085 (1992). [CrossRef]
  2. L. M. Mugnier, A. Blanc, and J. Idier, “Phase diversity: a technique for wave-front sensing and for diffraction-limited imaging,” Adv. Imaging Electron Phys. 141, 1–76 (2006). [CrossRef]
  3. D. R. Gerwe, M. M. Johnson, and B. Calef, “Local minima analysis of phase diverse phase retrieval using maximum likelihood,” The Advanced Maui Optical and Space Surveillance Technical Conference, Kihei, Hawaii (2008).
  4. R. W. Gerchberg and W. O. Saxton, “Phase retrieval by iterated projections,” Optik 35, 237–246 (1972).
  5. R. A. Gonsalves, “Small-phase solution to the phase-retrieval problem,” Opt. Lett. 26, 684–685 (2001). [CrossRef]
  6. C. U. Keller, V. Korkiakoski, N. Doelman, R. Fraanje, R. Andrei, and M. Verhaegen, “Extremely fast focal-plane wavefront sensing for extreme adaptive optics,” Proc. SPIE 8447, 844721 (2012). [CrossRef]
  7. F. Martinache, “Kernel phase in Fizeau interferometry,” Astrophys. J. 724, 464–469 (2010). [CrossRef]
  8. S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: a focal-plane wavefront sensor for real-time low-order sensing on faint sources,” Opt. Lett. 35, 3036–3038 (2010). [CrossRef]
  9. C. van der Avoort, J. J. M. Braat, P. Dirksen, and A. J. E. M. Janssen, “Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer–Zernike approach,” J. Mod. Opt. 52, 1695–1728 (2005). [CrossRef]
  10. I. Mocœur, L. M. Mugnier, and F. Cassaing, “Analytical solution to the phase-diversity problem for real-time wavefront sensing,” Opt. Lett. 34, 3487–3489 (2009). [CrossRef]
  11. W. J. Wild, “Linear phase retrieval for wave-front sensing,” Opt. Lett. 23, 573–575 (1998). [CrossRef]
  12. J. Antonello, M. Verhaegen, R. Fraanje, T. van Werkhoven, H. C. Gerritsen, and C. U. Keller, “Semidefinite programming for model-based sensorless adaptive optics,” J. Opt. Soc. Am. A 29, 2428–2438 (2012). [CrossRef]
  13. D. J. Lee, M. C. Roggemann, and B. M. Welsh, “Cramer-rao analysis of phase-diverse wave-front sensing,” J. Opt. Soc. Am. A 16, 1005–1015 (1999). [CrossRef]
  14. A. Tokovinin and S. Heathcote, “Donut: measuring optical aberrations from a single extrafocal image,” Publ. Astron. Soc. Pac. 118, 1165–1175 (2006). [CrossRef]
  15. A. Polo, S. F. Pereira, and P. H. Urbach, “Theoretical analysis for best defocus measurement plane for robust phase retrieval,” Opt. Lett. 38, 812–814 (2013). [CrossRef]
  16. ANSI, “Methods for reporting optical aberrations of eyes,” (2004).
  17. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005).
  18. F. Gustafsson, Statistical Sensor Fusion (Holmbergs i Malmö AB, 2010).
  19. G. H. Golub and C. F. Van Loan, Matrix Computations (Johns Hopkins University, 1996).
  20. R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. A 66, 207–211 (1976). [CrossRef]

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