OSA's Digital Library

Optics Express

Optics Express

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

Maximum-likelihood estimation of parameterized wavefronts from multifocal data

Julia A. Sakamoto and Harrison H. Barrett  »View Author Affiliations


Optics Express, Vol. 20, Issue 14, pp. 15928-15944 (2012)
http://dx.doi.org/10.1364/OE.20.015928


View Full Text Article

Acrobat PDF (1697 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A method for determining the pupil phase distribution of an optical system is demonstrated. Coefficients in a wavefront expansion were estimated using likelihood methods, where the data consisted of multiple irradiance patterns near focus. Proof-of-principle results were obtained in both simulation and experiment. Large-aberration wavefronts were handled in the numerical study. Experimentally, we discuss the handling of nuisance parameters. Fisher information matrices, Cramér-Rao bounds, and likelihood surfaces are examined. ML estimates were obtained by simulated annealing to deal with numerous local extrema in the likelihood function. Rapid processing techniques were employed to reduce the computational time.

© 2012 OSA

1. Introduction

Phase retrieval (PR) is a useful method for recovering the phase distribution in the pupil of an optical system from the irradiance distribution in the focal plane. However, the usual PR problem is ill posed, since both distributions are unrestricted 2D functions [1

1. I. S. Stefanescu, “On the phase retrieval problem in two dimensions,” J. Math. Phys. 26(9), 2141–2160 (1985). [CrossRef]

] and a single irradiance measurement does not ensure that the recovered phase is unique [2

2. J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7(3), 412–427 (1990). [CrossRef]

,3

3. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

]. A straightforward solution to avoid the ambiguities is to make multiple irradiance measurements near the focal plane, where the minimum number of planes is two.

Another approach to avoiding the phase ambiguities is to estimate the parameters describing the phase function, instead of the function itself [4

4. G. R. Brady and J. R. Fienup, “Improved optical metrology using phase retrieval,” 2004 Optical Fabrication & Testing Topical Meeting, OSA, Rochester, NY, paper OTuB3 (2004).

6

6. G. R. Brady, M. Guizar-Sicairos, and J. R. Fienup, “Optical wavefront measurement using phase retrieval with transverse translation diversity,” Opt. Express 17(2), 624–639 (2009). [CrossRef] [PubMed]

]. Parameterization of the phase is achieved with a set of expansion functions. Obvious choices for circular pupils are Zernike polynomials, although there are many other options.

In our method, the data consist of irradiance measurements in multiple planes near the focus of an aberrated optical element (Fig. 1
Fig. 1 Data-acquisition system for collecting multiple irradiance patterns near the focus of an optical element.
). From the multifocal data, we estimate the coefficients of phase polynomials in a wavefront expansion, as proposed by Brady and Fienup [4

4. G. R. Brady and J. R. Fienup, “Improved optical metrology using phase retrieval,” 2004 Optical Fabrication & Testing Topical Meeting, OSA, Rochester, NY, paper OTuB3 (2004).

6

6. G. R. Brady, M. Guizar-Sicairos, and J. R. Fienup, “Optical wavefront measurement using phase retrieval with transverse translation diversity,” Opt. Express 17(2), 624–639 (2009). [CrossRef] [PubMed]

], but with various extensions. We use maximum-likelihood estimation (MLE) for the estimation procedure, which enables us to optimize the data-acquisition system by analyzing the Fisher information matrix and associated Cramér-Rao lower bound, as well as likelihood surfaces. Additionally, we have developed our method to handle very large wavefront aberrations. To deal with the high computational demand of the estimation step, we employ rapid-processing techniques using dedicated computer hardware.

2. Propagation algorithm

The form of the diffraction integral, without the Fresnel approximation, used in this paper is
uz(r)=1iλzd2r0u0(r0)exp(ik|rr0|2+z2),
(1)
where r0 and r are two-dimensional vectors in the pupil and image planes, respectively [7

7. H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, 2004).

,8

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

]. It is a version of the Huygens wavelet formulation when the obliquity factor in the integral is neglected and distance is replaced by z in the denominator. Thus, the field at each observation point in the image plane is the sum of an infinite number of secondary waves emanating from the aperture. Since the integrand of the double integral in Eq. (1) can vary rapidly over the integration domain, especially when the aberrations are large, a significant amount of sampling in the aperture is required to accurately predict the irradiance distribution in the image plane. Consequently, numerical evaluations of the Huygens diffraction integral are prohibitively time-consuming and an approximation to this integral is necessary to reduce the computing problem.

Note that both the Huygens wavelet formation and the Fresnel approximation assume scalar diffraction theory, which we assume to be adequate for our method of wavefront measurement with multifocal data, even when large aberrations are present.

Suppose a converging wave has wavefront error W(r0) in the exit pupil of an optical system. Then the field in the exit pupil is given by
u0(r0)=exp[ikW(r0)]exp[ikRf(r0)]Rf(r0),
(2)
where Rf(r0)=|r0|2+f2 and f is defined as the radius of curvature of the unaberrated wave in the exit pupil (i.e., the distance between the exit pupil and paraxial focus). Inserting Eq. (2) into Eq. (1) leads to
uz(r)=1iλzxpd2r0exp[ikW(r0)]exp[ikRf(r0)]Rf(r0)exp(ik|rr0|2+z2),
(3)
or equivalently,
uz(r)=1iλzxpd2r0exp[ikW(r0)]Rf(r0)exp(ikrr0z)×exp[(ik|rr0|2+z2|r0|2+f2+rr0z)],
(4)
where the last exponential term is part of the integrand. A binomial expansion for z>|rr0|and f>r0 leads to
|rr0|2+z2|r0|2+f2+rr0z=zf+r22z+r022(1z1f)+HOT,
(5)
where HOT denotes the higher-order terms in r0:
HOT=|rr0|48z3+|r0|48f3+|rr0|616z5|r0|616f5+....
(6)
Combining Eq. (4) and Eq. (5) results in
uz(r)=A(r)xpd2r0exp[ikW(r0)]Rf(r0)exp(ikrr0z)exp{ik[r022(1z1f)+HOT]},
(7)
where the function A(r) is given by

A(r)=1iλzexp[ik(zf+r22z)].
(8)

If the higher-order terms in Eq. (7) can be ignored, that equation can be represented by a 2D Fourier transform, known as the Fresnel approximation:
uz(r)A(r)F2{1Rf(r0)exp[ikW(r0)]exp[ikr022(1z1f)]}ρ=r/λz,
(9)
where the spatial frequency is ρ=r/λz. The corresponding irradiance under this approximation is given by

I(r)=|uz(r)|21λ2z2|F2{1Rf(r0)exp[ikW(r0)]exp[ikr022(1z1f)]}ρ=r/λz|2.
(10)

The approximation in Eq. (9) breaks down as the numerical aperture of the optical system increases. Under circumstances when the higher-order terms are not negligible, we might consider including, say, the fourth-order terms in Eq. (6). A limitation is that r and r0 are inseparable in these terms, so including them in the integral means that we cannot take advantage of the FFT. Since a brute-force computation of the diffraction integral is too computationally expensive and impractical, the best we can do is minimize the terms in Eq. (6) by considering planes sufficiently close to nominal focus, so that zf and r is small.

3. Parameterized wavefront description

The unknown function of interest in these equations is the wavefront error W(r0), which we represent in a parameterized form in the exit pupil of the optical system. Our approach is based on the fundamental assumption that the continuous wavefront can be approximated to sufficient accuracy by a finite set of expansion functions. We choose to represent this function by expanding it in some number of Zernike polynomials,
W(r0)n=1NαnZn(r0),
(11)
where Zn(r0) is the nth polynomial with coefficient αn.

The parameters to be estimated are the Zernike coefficients {αn, n = 1,…, N), but an important determination is the number of coefficients necessary for an accurate representation of the wavefront. Even if a small number of terms is used in the expansion, this does not imply that the wavefront aberration is small; representing large wavefront errors simply requires large coefficients. In our approach, we assume that the wavefront is smoothly varying, so that sufficient accuracy can be achieved with a relatively small value of N.

4. Maximum-likelihood estimation

Immediate advantages of ML estimation are that it is efficient (i.e., unbiased and yields the best possible variance), if an efficient estimator exists, and that it is asymptotically efficient as more or better data are acquired. Furthermore, its performance can be analyzed with the Fisher information matrix (FIM), which can be used to design and optimize the system that acquires the data to be used as input to the estimation routine. One major challenge in ML estimation is that an accurate probability model must be used that includes all sources of randomness (e.g., photon noise and electronic noise).

A significant limitation to be overcome is that the ML estimation step is very time-consuming, so making it practical requires rapid processing techniques and dedicated computer hardware. Improvement of computational time can be achieved by parallelizing the propagation algorithm describing the forward model of the system. Parallel algorithms can be implemented on a variety of hardware platforms, including the graphics processing unit (GPU) in video cards, capable of massively parallel high-performance computing in scientific and engineering fields. Modern GPU programming is accomplished using an extension to the C programming language, CUDA (“Compute Unified Device Architecture”), which provides software development tools and allows functions in C to be implemented on a GPU’s multiple stream processors.

Since the logarithm increases monotonically with its argument, Eq. (12) can be rewritten as
θ^ML=argmaxθln pr(g|θ).
(13)
For a purely Gaussian noise model with independent and identically distributed (i.i.d.) detector elements, it can be shown that ML estimation according to Eq. (13) reduces to nonlinear least-squares fitting between the data and the output of the optical design program:
θ^ML=argminθm=1M[gmg¯m(θ)]2.
(14)
Since the data-acquisition system described in this work was not photon-starved, we assumed that it was dominated by electronic noise, which is adequately described by Gaussian statistics.

5. Fisher information and Cramér-Rao bounds

For Gaussian i.i.d. data, the FIM components reduce to
Fjk=1σ2m=1Mg¯m(θ)θjg¯m(θ)θk,
(16)
which depend on the average data vector evaluated at the set of parameters and are inversely proportional to the variance in the detector elements [7

7. H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, 2004).

,10

10. H. H. Barrett, C. Dainty, and D. Lara, “Maximum-likelihood methods in wavefront sensing: stochastic models and likelihood functions,” J. Opt. Soc. Am. A 24(2), 391–414 (2007). [CrossRef] [PubMed]

].

Inversion of the FIM provides the Cramér-Rao lower bound (CRB) on the variance (i.e. the theoretical minimum possible variance) of the parameter estimates. It is well-documented in the literature that the variance of any unbiased estimate obeys the Cramér-Rao inequality [11

11. C. R. Rao, “Information and accuracy attainable in the estimation of statistical parameters,” Bull. Calcutta Math. Soc. 37, 81–91 (1945).

,12

12. H. Cramér, Mathematical Methods of Statistics (Princeton University Press, 1946).

], written as
[Kθ^]pp=Var{θ^p}[F1]pp,
(17)
where θp is the pth parameter and Kθ^ is the covariance matrix of the estimates. For an efficient estimator, the inverse of the FIM is the covariance matrix of the parameter estimates; thus, its off-diagonal elements represent coupling between these estimates.

In general, the FIM can be computed for any system configuration and therefore be used to design and optimize the system that acquires the data to be used as input to the estimation procedure. Additionally, data from multiple image planes may serve to reduce parametric coupling, which could improve parameter estimability.

6. Numerical results

6.1 Test lens description

For the numerical proof-of-principle system, we chose a large nine-surface unit lens for use in a larger assembly. We operated the rotationally-symmetric lens at finite conjugates by placing an on-axis point source at 113.0 mm from the entrance pupil. According to ZEMAX, the exit pupil diameter and working f-number were Dxp = 88.41 mm and f/#w = 1.796, respectively, at the wavelength λ = 0.6328 μm. The peak-to-valley wavefront error was 149.1λ.

We chose to expand the wavefront in terms of the Fringe Zernike polynomials, which are identical to the standard polynomials, except for the indexing format and the order in which they are listed. The Fringe Zernike coefficients describing W(r0) in the exit pupil were calculated in ZEMAX for N = 37, the maximum number of terms (Table 1

Table 1. Fringe Zernike coefficients {αn, n = 1,…, 37}, peak-to-valley, RMS, and variance, provided by ZEMAX for the highly aberrated test lens. Unlisted coefficients are zero.

table-icon
View This Table
| View All Tables
). Since the system has rotational symmetry, the coefficients for all non-rotationally symmetric Zernike terms are zero; thus non-zero coefficients correspond to piston, defocus, and various orders of spherical aberration. Table 2

Table 2. Fringe Zernike Polynomials {Zn, n = 1,…, 9, 16, 25, 36, 37}.

table-icon
View This Table
| View All Tables
provides the polynomial expressions for {Zn, n = 1,…,9, 16, 25, 36, 37}, corresponding to the rotationally symmetric terms and several low-order terms.

6.2 Irradiance data

A critical design option is the number and location of planes at which to measure the irradiance. It is well known that two-plane measurements are sufficient to determine the pupil phase, so to minimize computation time, we selected two output planes just before paraxial focus, z = z1 = zf – 0.25 mm and z = z2 = zf – 0.43 mm, where the paraxial focal plane is located at z = zf = 157.8 mm and the exit pupil lies in the z = 0 plane. So that each irradiance distribution provided unique information, special care was taken to select planes sufficiently far from each other, yet close enough to the nominal focus for the approximation in Eq. (9) to hold.

For both planes, we determined that the minimum level of pupil sampling to accurately represent the irradiance data was 1024 × 1024. Essentially, output planes farther from focus and larger aberrations in the test lens require greater sampling. We added zero-padding for a final FFT grid size of 2048 × 2048, which was used to compute the irradiance data in Fig. 2
Fig. 2 Detector data for the highly aberrated test lens using a pupil sampling of P = 1024 at image plane: (a) z = z1 and (b) z = z2.
. The detector element spacing was roughly 0.56 μm based on these grid sizes.

6.3 Fisher information and Cramér-Rao bounds

We computed the FIM using Eq. (16) for the parameters we wanted to estimate, the Fringe Zernike coefficients {αn, n = 2,…, 9, 16}, shown in Fig. 3(a)
Fig. 3 (a) FIM and (b) its inverse, for Fringe Zernike coefficients {αn, n = 2,…, 9, 16} in the exit pupil of the highly aberrated test lens. Shown on a logarithmic scale.
on a logarithmic scale. Although the off-axis coefficients included here are zero in this idealized case, in a real system, they may become non-zero if there are any system misalignments. Since α1 corresponds to piston and does not influence the irradiance data, we had no interest in estimating it and disregarded it in the FIM. The variance σ2 in the detector elements was based on a peak SNR of 104. We evaluated the FIM components at θ, based on the values in Table 1. Note that the FIM is a 9 × 9 symmetric matrix and is order-specific, with the parameter indices on the x and y axes. The jkth entry has units of

unitsofFjk=1(unitsofθj)(unitsofθk).
(18)

High diagonal values in the FIM indicated that there was substantial information in the data for each parameter. However, the pronounced off-diagonal structure revealed strong parametric coupling between the pairs (2,7), (3,8), (4,9), (4,16), and (9,16), which would only confound the estimation problem. Note that the degree of coupling between two parameters is proportional to the magnitude of the respective FIM component.

We computed the inverse of the FIM, shown in Fig. 3(b), and read off its diagonal components to determine the CRB for the parameters, whose square-root is provided in Table 3

Table 3. Square-root of the CRB for Fringe Zernike coefficients {αn, n = 2,…, 9, 16} in the exit pupil of the highly aberrated test lens at λ = 0.6328 μm. Units are in waves λ.

table-icon
View This Table
| View All Tables
. The diminutive values in (CRB)1/2, on the order of 10−10 to 10−9, permit the estimation of the wavefront parameters to very high precision, provided that the forward model is exact.

6.4 Likelihood surfaces

When solving optimization problems, it helps to have a strong sense of the objective function to be optimized. This is particularly useful when fitting nonlinear, multivariate functions that may be plagued with numerous local extrema.

In the context of ML estimation, a plot of pr(g|θ ) versus θ for a particular g is called the likelihood surface for that data vector. Although we are minimizing the sum-of-squares (i.e., negative likelihood), we will refer to this objective function as the likelihood surface.

For a vector set of N parameters, the likelihood surface exists in an N-dimensional hyperspace. However, we are restricted to plotting the likelihood surface while varying up to two parameters at a time. We selected a handful of pairs of parameters for the following plots (Fig. 4
Fig. 4 Likelihood surface plotted along two axes at a time for pairs of parameters: (a) α9 and α16, (b) α4 and α9, (c) α3 and α8, (d) α7 and α9, (e) α7 and α16, and (f) α2 and α3.
) and applied the ranges given in Table 4

Table 4. Range in likelihood surface plots for Fringe Zernike coefficients {αn, n = 2,…, 9, 16} in the exit pupil of the highly aberrated test lens. Units are in waves.λ.

table-icon
View This Table
| View All Tables
. In each plot, the fixed parameters were set to the true values underlying the data.

6.5 Maximum-likelihood estimates

We ran 12 optimization trials, with each trial representing a distinct noise realization for the same wavefront parameters, and implemented the same estimation procedure on each. Figure 5
Fig. 5 12 simulated annealing trials for the estimation of wavefront parameters in the exit pupil of the highly aberrated test lens, plotting the optimal cost function versus iteration number on a log-log scale.
illustrates the optimal cost function versus iteration number for the trials, starting at the end of the first temperature phase. The average final cost and number of temperature phases were 147.7 and 87.3, respectively. This number of phases is equivalent to a final temperature of 116.1.

Final average estimates for Zernike coefficients {αn, n = 2,…, 9, 16} are given in Table 5, along with the respective standard deviation. Both the bias and variance are small for every parameter, while the bias is always within one standard deviation. Based on the small biases, it is not surprising that the true and estimated irradiance patterns are virtually indistinguishable (Fig. 6
Fig. 6 Comparison between the true and estimated irradiance patterns for the highly aberrated test lens.
).

An iteration is defined as one cycle through every parameter, and there were 200 iterations, or 1.8 × 103 forward propagations, per temperature phase. For the specified pupil sampling and FFT grid size, the computation time was 790 ms per forward propagation, which included both output planes. So, the average of 87.3 temperature phases per trial corresponds to 34.5 hours of computation time. We carried out this study using an NVIDIA Tesla C1060 GPU, whose peak double-precision (DP) performance is 78 GFLOPS. Had we used a Tesla C2075 model, offering 515 GFLOPS of DP power, the computation time per trial would have been roughly 5.2 hours. Further improvement can be achieved with a cluster of GPUs. For instance, the VSC455 V8 GPU workstation by Velocity Micro combines 8 Tesla C2075s for over 4 TFLOPS of DP power, which would turn the 5.2 hours into just 40 minutes. Another way to improve computational time involves beamlet-based propagation algorithms, such as CODE V Beam Synthesis Propagation (BSP), which can provide accurate and efficient modeling of the wave propagation.

7. Experimental results

7.1 System configuration

We obtained multifocal irradiance data in the focal region of a relatively benign spherical test lens at finite conjugates and with a HeNe source (λ = 0.6328 μm). The distance from the on-axis point source, created with a 40X microscope objective and 10-μm pinhole, to the test lens was 457 mm. To generate an increase in information by illuminating substantially more detector elements, we magnified the intermediate image with a relay lens placed just before the CCD. This had the added benefit of eliminating CCD saturation without the use of neutral density filters. The imaging lens was a Nikon 40X microscope objective (NA = 0.95), configured at the design conjugates. When used at the proper conjugates, the manufacturer claims the objective lens is corrected for spherical aberration. This particular objective was infinity-corrected with an optical tube length of 160 mm. To maintain these conjugates, we placed the imaging lens and CCD on a translation stage, allowing us to scan through the focal region of the test lens. We used a CCD with a fine pixel size of 4.4 μm for high information yield. In general, doubling the number of pixels also doubles the Fisher information, as indicated in Eq. (16), provided that the detector variance remains the same.

7.2 Test lens description

The test lens for this study was a double-convex spherical lens with a diameter of Dlens = 25 mm and radius of curvature of R1 = -R2 = 76.66 mm. According to ZEMAX, the exit pupil diameter and working f-number was Dxp = 25.39 mm and f/#w = 3.471, respectively. It was imperative that the image-space NA of the test lens (NA = 0.14) was less than that of the imaging lens to avoid information loss. The position of the paraxial focal plane was z = zf = 90.83 mm, where the z = 0 plane contained the exit pupil.

The ZEMAX calculations of the Fringe Zernike coefficients (N = 37) in the exit pupil of the lens are provided in Table 6

Table 6. Fringe Zernike coefficients {αn, n = 1,…, 37}, peak-to-valley, RMS, and variance, provided by ZEMAX for the highly aberrated test lens. Unlisted coefficients are zero.

table-icon
View This Table
| View All Tables
. With a smaller peak-to-valley wavefront error of 30.6λ, we anticipated less stringent requirements on the pupil sampling in our propagation algorithm. The minimum pupil sampling level was 256 × 256 for this lens, which we zero-padded for an FFT grid size of 1024 × 1024. This resulted in oversampling of the irradiance distribution, however, this only increased the computational time and did not affect the outcome of the study.

7.3 Irradiance data

Figure 7
Fig. 7 Experimental data for the spherical test lens for image planes: (a) z = z1 and (b) z = z2. Scale bar corresponds to the intermediate image plane just before the imaging lens.
displays the physical data that we used as input to our optimization algorithm, consisting of two image planes just before paraxial focus, where the scale bar corresponds to the intermediate image plane between the test lens and objective lens. In Section 7.6, we discuss the estimation of nuisance parameters in the system, namely, the exact image plane locations and the true magnification of the objective lens.

7.4 Fisher information and Cramér-Rao bounds

We computed the FIM and its inverse (Fig. 8
Fig. 8 (a) FIM and (b) its inverse, for Fringe Zernike coefficients {αn, n = 2,…, 9, 16} in the exit pupil of the spherical test lens. Shown on a logarithmic scale.
) for the Fringe Zernike coefficients to be estimated, {αn, n = 2,…, 9, 16}, based on the method in Section 6.3. The structure of these matrices is comparable to that for the highly aberrated lens, including strong coupling among the rotationally-symmetric terms, α4, α9, and α16. Since the CRB takes on diminutive values, the variance from detector noise is unlikely to preclude high-precision estimates.

Table 7. Square-root of the CRB for Fringe Zernike coefficients {αn, n = 2,…, 9, 16} in the exit pupil of the highly aberrated test lens at λ = 0.6328 μm. Units are in waves λ.

table-icon
View This Table
| View All Tables

7.5 Nuisance parameters

Nuisance parameters are defined as parameters that affect the data and are therefore fundamental to the probability model, but are not useful to the estimation task. Suppose α denotes the wavefront parameters of interest, β denotes the wavefront parameters of no immediate interest, and χ denotes all other nuisance parameters, so that the entire vector parameter can be written as θ = (α, β, χ) t. Therefore, α is comprised of Zernike coefficients {αn, n = 2,…, 9, 16}, β is comprised of all remaining coefficients, and χ consists of unknown system parameters to be discussed next. We dealt with β by replacing it (up to the ZEMAX limit of N = 37) with the respective design coefficients in Table 6, denoted as β0, then separately estimated χ prior to the primary estimation task. Finally, we set pr(g|α,β,χ)pr(g|α,β0,χ^) and proceeded to estimate α.

Two of the critical nuisance parameters in the system were the absolute image-plane positions, z1 = zf + Δz1 and z2 = zf + Δz2. During acquisition of the detector data, it was difficult to accurately pinpoint the paraxial focal plane, z = zf, which served as the reference plane in our studies. Another nuisance parameter was the true magnification of the microscope objective just before the detector, denoted as M1 and M2 for the first and second image, respectively. To achieve the design magnification of 40X, the detector must be placed 160 mm from the rear principal plane of the objective lens, however, we had no knowledge of the exact location of this plane.

Since the detector element size was 4.4 μm, the effective pixel size after magnification was Δxd = 0.1101 μm and Δxd = 0.1098 μm. This resulted in sampling the irradiance distribution below the diffraction limit of approximately 0.5 λ/NA = 0.33 µm (i.e., oversampling by a factor of three). However, the excess magnification also allowed us to avoid detector saturation without the use of additional optics, as mentioned in Sec. 7.1.

7.6 Maximum-likelihood estimates

To compare the physical data with the output of the optical design program during optimization, we first interpolated the data, so that its coordinates matched that of the FFT output. Then we normalized the irradiance pattern in both the data and FFT output, which is analogous to normalizing the power of the source beam.

Final average estimates and standard deviations for Zernike coefficients {αn, n = 2,…, 9, 16} are given in Table 8. Without knowledge of the true values or a gold standard, we have no basis for evaluating biases in the estimates, though the estimated values are within one standard deviation from the design values. Once again, large departures from the CRB may have resulted from the numerous local basins in the cost function.

There were 250 iterations or 2.25 × 103 forward propagations in each temperature phase. For the given FFT grid size, the computation time was 308 ms per forward propagation. On average, there were 105 temperature phases per trial, corresponding to 20.2 hours of computation time with a single NVIDIA Tesla C1060. This amounts to 3.1 hours with a Tesla C2075 model, or 23.3 minutes with a VSC455 V8.

8. Conclusions

In both the numerical and experimental studies, the data-acquisition method involved multiple irradiance patterns collected near the focus of an optical system for the purpose of estimating the pupil phase distribution. We considered various approaches for a suitable propagation algorithm to accurately model the wave propagation and developed the mathematical framework for an aberrated wave emerging from the exit pupil, where we parameterized the wavefront using expansion functions, particularly, Fringe Zernike polynomials. To substantially reduce the computation time, we implemented parallel processing with a state-of-the-art GPU.

After discovering numerous local extrema in the likelihood surfaces, we chose simulated annealing for the estimation of selected Zernike coefficients. In the numerical study with the highly aberrated lens, the estimate biases were negligible, probably because of the lack of systematic errors in the estimation procedure. Although the variances were fairly small, they were far from the CRB due to entrapment in local basins in the cost function.

In the experimental study, we used a benign test lens with significantly less aberrations and a larger working f-number of f/#w = 3.47. Since the imaging (i.e., relay) lens was much faster with f/# ≈0.53 (NA = 0.95), it should not have resulted in information loss from the suppression of high spatial frequencies. However, this may become an issue if the f-number of the test lens is low enough. To avoid having to include the imaging lens in the forward model here, one solution is to use the translation stage to place a ground glass in any plane where an irradiance measurement is desired, then relay the incoherent irradiance, instead of the field, to the detector. Note that it is necessary to have a rotating ground glass to decorrelate the resulting speckle noise. An alternative to a rotating diffuser is a liquid-crystal diffuser operated in a dynamic scattering mode.

Acknowledgments

This work was supported by Canon, Inc. and the National Institutes of Health under grant number R37 EB000803. We would also like to thank Jin Park, Luca Caucci, and Akinori Ohkubo for their valuable input.

References and links

1.

I. S. Stefanescu, “On the phase retrieval problem in two dimensions,” J. Math. Phys. 26(9), 2141–2160 (1985). [CrossRef]

2.

J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7(3), 412–427 (1990). [CrossRef]

3.

M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

4.

G. R. Brady and J. R. Fienup, “Improved optical metrology using phase retrieval,” 2004 Optical Fabrication & Testing Topical Meeting, OSA, Rochester, NY, paper OTuB3 (2004).

5.

G. R. Brady and J. Fienup, “Phase retrieval as an optical metrology tool,” Optifab: Technical Digest, SPIE Technical Digest TD03, 139–141 (2005).

6.

G. R. Brady, M. Guizar-Sicairos, and J. R. Fienup, “Optical wavefront measurement using phase retrieval with transverse translation diversity,” Opt. Express 17(2), 624–639 (2009). [CrossRef] [PubMed]

7.

H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, 2004).

8.

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

9.

A. Corana, M. Marchesi, C. Martini, and S. Ridella, “Minimizing multimodal functions of continuous variable with the ‘simulated annealing’ algorithm,” ACM Trans. Math. Softw. 13(3), 262–280 (1987). [CrossRef]

10.

H. H. Barrett, C. Dainty, and D. Lara, “Maximum-likelihood methods in wavefront sensing: stochastic models and likelihood functions,” J. Opt. Soc. Am. A 24(2), 391–414 (2007). [CrossRef] [PubMed]

11.

C. R. Rao, “Information and accuracy attainable in the estimation of statistical parameters,” Bull. Calcutta Math. Soc. 37, 81–91 (1945).

12.

H. Cramér, Mathematical Methods of Statistics (Princeton University Press, 1946).

OCIS Codes
(100.5070) Image processing : Phase retrieval
(120.3940) Instrumentation, measurement, and metrology : Metrology
(120.5050) Instrumentation, measurement, and metrology : Phase measurement

ToC Category:
Image Processing

History
Original Manuscript: February 22, 2012
Revised Manuscript: May 30, 2012
Manuscript Accepted: June 6, 2012
Published: June 28, 2012

Citation
Julia A. Sakamoto and Harrison H. Barrett, "Maximum-likelihood estimation of parameterized wavefronts from multifocal data," Opt. Express 20, 15928-15944 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-14-15928


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. I. S. Stefanescu, “On the phase retrieval problem in two dimensions,” J. Math. Phys.26(9), 2141–2160 (1985). [CrossRef]
  2. J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A7(3), 412–427 (1990). [CrossRef]
  3. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am.73(11), 1434–1441 (1983). [CrossRef]
  4. G. R. Brady and J. R. Fienup, “Improved optical metrology using phase retrieval,” 2004 Optical Fabrication & Testing Topical Meeting, OSA, Rochester, NY, paper OTuB3 (2004).
  5. G. R. Brady and J. Fienup, “Phase retrieval as an optical metrology tool,” Optifab: Technical Digest, SPIE Technical DigestTD03, 139–141 (2005).
  6. G. R. Brady, M. Guizar-Sicairos, and J. R. Fienup, “Optical wavefront measurement using phase retrieval with transverse translation diversity,” Opt. Express17(2), 624–639 (2009). [CrossRef] [PubMed]
  7. H. H. Barrett and K. J. Myers, Foundations of Image Science (Wiley, 2004).
  8. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company Publishers, 2005).
  9. A. Corana, M. Marchesi, C. Martini, and S. Ridella, “Minimizing multimodal functions of continuous variable with the ‘simulated annealing’ algorithm,” ACM Trans. Math. Softw.13(3), 262–280 (1987). [CrossRef]
  10. H. H. Barrett, C. Dainty, and D. Lara, “Maximum-likelihood methods in wavefront sensing: stochastic models and likelihood functions,” J. Opt. Soc. Am. A24(2), 391–414 (2007). [CrossRef] [PubMed]
  11. C. R. Rao, “Information and accuracy attainable in the estimation of statistical parameters,” Bull. Calcutta Math. Soc.37, 81–91 (1945).
  12. H. Cramér, Mathematical Methods of Statistics (Princeton University Press, 1946).

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