## Abstract

The aliasing effect in the discrete Fourier transform inherent will impose a serious detrimental effect on conventional phase retrieval measurement accuracy with under-sampled intensity. In this Letter, we describe a modal-based nonlinear optimization phase retrieval approach that is capable of retrieving wavefront measurements using under-sampled intensities. The extended Nijboer–Zernike theory is introduced to establish an analytic solution between wavefront phase and intensity image, and then nonlinear optimization is further utilized to solve wavefront aberration coefficients from under-sampled intensity data. The feasibility and accuracy of the algorithm are verified by simulations and experiments. This is a promising method that is especially suitable for full field phase recovery of optical systems with a relatively high numerical aperture.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

Phase retrieval is a simple and effective method to obtain the missing phase from intensity measurements. The computational procedure in which the desired plane is propagated to the measurement plane based on the discrete Fourier transform (DFT) is always involved in iterative phase retrieval algorithms [1]. The numerical propagation for the phase retrieval algorithm, which is generally computed using a fast Fourier transform (FFT), needs to satisfy the Nyquist criterion [2]. The intensity pattern should be adequately sampled by the detector which is determined using a calculated sampling parameter, $Q$, given by [3]

where $\lambda$ is the wavelength of illumination, $ F $-number is $F/\# = L/D$, $ L $ is the propagation distance from the desired plane to the detector plane, $D$ is the testing aperture, and ${d_u}$ is the detector pixel pitch. $Q \ge 2$ is the measure that the intensity pattern is adequately (Nyquist) sampled. For the general optical testing system, this limits the use of most phase retrieval algorithms to measure the full pupil wavefront with a more moderate $F/\#$ (${5}\sim 10$).To overcome this limitation, the transverse translation diverse phase retrieval (TTDPR) [3–5] algorithm is applied to the large numerical aperture (NA) optical element testing. A mask is placed near the aperture of the desired plane to reduce its NA and form an adequately sampled intensity pattern. TTDPR technique has good measurement accuracy and reproducibility for optical metrology [5]. While this method is capable of testing wavefront over the arbitrarily large aperture, it requires additional translation mask devices and longer data collection time.

For under-sampled broadband images, specialized algorithms have been demonstrated without additional hardware [6,7]. For optical wavefront metrology using a coherent optical system, one approach is that making additional point spread function (PSF) measurements with random pointing fluctuations or intermediate amounts of defocus improves testing results accuracy [8]. Another method changes conjugates used in the test to create a beam, which is adequately sampled by the detector array [9]. In addition, the under-sampled intensity via super-resolution preprocessing can also be applied to the traditional phase diversity algorithm [10]. However, they could induce additional errors that degrade measurement accuracy.

The sampling requirement in the DFT matrix calculation is the main limitation for the traditional phase retrieval algorithm to measure the wavefront of elements with relatively large NA. If the Fourier transform can be represented in the analytic form, the diffraction field can be flexibly sampled according to the pixel size of the detector, then the under-sampled restriction problem would be avoided in the phase retrieval testing system. Here, we describe a diffraction calculation approach instead of the DFT matrix calculation for which the extended Nijboer–Zernike (ENZ) basis function is used for modeling the under-sampled PSFs [11,12]. The ENZ theory derives the analytic expression of the two-dimensional PSF in the focal region through the Debye diffraction integral and the Zernike polynomials expansion of the wavefront function in the desired plane. The phase retrieval solved via a linear system of equations is a general method based on ENZ theory [13,14]. Those methods are proposed for high NA projection lens aberration reconstruction. However, the accuracy algorithm working in a general optical testing system has not been adequately explored in real experiments.

The nonlinear optimization phase retrieval (NOPR) algorithm is superior to the iterative transform algorithm and has high flexibility suitable for different circumstances [15]. In this Letter, the modal-based nonlinear optimization phase retrieval algorithm using the ENZ theory (ENZ-NOPR) is proposed. Compared with the classical NOPR algorithm calculating the gradients with a pair of DFT propagation operations, we derive the gradient calculation suitable for the ENZ-NOPR algorithm, which involves using a nonlinear least-square-error optimization method based on the ENZ basis function among measurement planes.

First, we briefly introduce our optical testing model based on the ENZ theory. The complex field $P({\rho ,\theta})$ (both amplitude and phase) in the desired plane can be modally parameterized by the coefficients of the Zernike polynomials:

We can now formulate the phase retrieval wavefront measurement problem based on the ENZ theory. For convenience, we introduce some notation to describe our measurement model with matrices. Column vectors denoted in bold by ${\hat {\boldsymbol I}_k}$, ${{\boldsymbol U}_k}$, and ${\boldsymbol C}_k^{{nm}}$ are raster-scanned from their two-dimensional counterparts, $\hat I_k^{{nm}}({{r_k},{f_k}})$, ${{\boldsymbol U}_k}({{r_k},{\varphi _k},{f_k}})$, and $C_k^{{nm}}({{r_k},{\varphi _k},{f_k}})$. ${\hat I_k}$ is the normalized measured intensity ${I_k}$. $\boldsymbol \beta$ and ${{\boldsymbol C}_k}$ are the set of ${\boldsymbol \beta}_n^m$ and ${\boldsymbol C}_k^{{nm}}$ for the optical testing system. At the $j$th iteration process, the error metric for this wavefront testing system is defined as

The main workflow of the ENZ-NOPR algorithm is illustrated in Fig. 1. The reconstruction calculation is started with a weighted ideal lens ($\beta _0^0 = 0.1$). The algorithm framework has the following procedures:

- (1) Set the minimum threshold of computed error metric $\varepsilon$, the number of patterns $N$, the initial step length ${h_0}$, $j = 1$, and $k = 1$. And the reconstruction is started with ${{\boldsymbol \beta}_0} = {\{{0.1,0,0, \ldots ,0} \}^{T}}$;
- (2) Calculate corresponding gradient function$$\frac{{\partial {{\boldsymbol E}_{{j,k}}}}}{{\partial {{\boldsymbol \beta}_{j}}}} = {\boldsymbol C}_{k}^\dagger {{\boldsymbol S}_{k}}\left({\hat {\boldsymbol U}_{{j},{k}}^{c} - {\boldsymbol U}_{{j},{k}}^{c}} \right),$$with$$\hat {\boldsymbol U}_{{j},{k}}^{c} = \sqrt {\hat {{{\boldsymbol I}_{k}}}} \frac{{{\boldsymbol U}_{{j},{k}}^{c}}}{{\left| {{\boldsymbol U}_{{j},{k}}^{c}} \right|}},$$where ${\boldsymbol C}_k^\dagger = {({{\boldsymbol C}_k^{T}{{\boldsymbol C}_k}})^{- 1}}{\boldsymbol C}_k^{T}$. The search direction ${{\boldsymbol D}_j}$ we used is given via the conjugate gradient method [19] in the first 30 iterations and then via the steepest-descent method [20];
- (5) $k = {\rm rem}({j,N}) + 1,j = j + 1$, repeat steps (2)–(4) until satisfying ${{\boldsymbol E}_{j,k}} \lt \varepsilon ,({k = 1,2, \cdots ,N})$.
- (6) The testing wavefront $W({\rho ,\theta})$ is calculated with Eq. (2).

In our numerical and real experiments, because the number of pixels of used diffraction patterns is far more than the number of terms of basis function, ${\boldsymbol C}_k^\dagger$ exists [21]. Therefore, the derivatives of basis function coefficients can be effectively calculated, rather than encountering the aliasing problem in the DFT inherent. Compared with the classical phase retrieval algorithm, the computational complexity of ENZ-NOPR is degraded because there is no backward propagation to the desired plane, and the dimensions of the intensity grid are smaller for every optimization estimation iteration. Using all the collected intensities in sequence helps to eliminate the stagnation problem and improve the robustness.

We numerically investigate the performance of this proposed ENZ-NOPR algorithm using under-sampled intensity data at first. Two wavefront cases, consisting of a simple wavefront fitted by the first 36 terms in the Zernike polynomial and a general wavefront fitted by the first 190 terms in the Zernike polynomial, are used as the ground truth. For ${Q} = 0.75$, the accuracy of ENZ-NOPR algorithm is numerically verified. The collected intensity patterns are sampled with $128 \times 128$ pixels, and the tested wavefront is sampled with $256 \times 256$ pixels, respectively. The white Gaussian noise is added and the SNR is 35 dB for every intensity pattern. The intensity value for every normalized measured diffraction pattern below 0.1 is removed using matrix ${{\boldsymbol S}_k}$. When the two wavefront samples are reconstructed using the ENZ-NOPR algorithm, the 231 terms in the ENZ basis function ${\boldsymbol C}_k^{{nm}}$ for a simple wavefront and the 351 terms in the ENZ basis function ${\boldsymbol C}_k^{{nm}}$ for a general wavefront are chosen, respectively. The reconstructed results are shown in Figs. 2(b) and 2(e), respectively. The root-mean-square error (RMSE) is calculated, which quantifies the difference between the computed wavefront and the corresponding ground truth or interferometric data. For the ENZ-NOPR and variable sampling mapping-hybrid diversity algorithm (VSM-HDA) [7], the RMSE for two samples recovery is listed in Table 1. Whether the wavefront measured is a simple or general wavefront, the reconstructed phase matches better than VSM-HDA, indicating that the proposed ENZ-NOPR algorithm is effective in the reconstruction of the wavefront using the under-sampled intensity.

To verify the robustness of ENZ-NOPR, simulations are performed on 60 different wavefront realizations for four different Q values. For all simulations, the SNR is 35 dB for every intensity pattern. The multiple-images amplitude-point and phase-Zernike nonlinear optimization algorithm (multi-NOPZ) [18] is selected as a comparison algorithm, and its retrieved results are shown in the last set of Fig. 3. The rest of the results as shown in Fig. 3, which is retrieved via ENZ-NOPR with different Q values, has good accuracy as well as multi-NOPZ based on adequate intensity sampling. Moreover, the computational time spent by ENZ-NOPR is about a quarter of NOPR based on FFT.

The experiments were carried out to verify the feasibility and accuracy of the ENZ-NOPR method using under-sampled intensities. The layout of the experimental setup is depicted in Fig. 4. The beam with $\lambda = 632.8\;{\rm nm} $ illuminates the plate and a circular aperture stop. The camera, which is a beam profiler (BGP-USB-SP620U) with 12-bit depth and $4.4\; {\unicode{x00B5}{\rm m}} \times 4.4\; {\unicode{x00B5}{\rm m}}$ pixel size, is mounted on a motorized translate stage to allow PSFs to be made in various defocus planes. The condensing lens has a focal length of 151.74 mm and a diameter of 22.9 mm, the effective $F/\#$ of the unobstructed system would be 6.62, yielding ${Q} = 0.953$, according to Eq. (1), which is under-sampled for intensities in this phase retrieval testing system. Three adopted diffraction patterns are measured at defocus distance $\Delta z \in [{- 1, - 2, - 3}]\;{\rm mm}$ along the $z$-axis. And there is a sampling parameter $Q = 0.934$ at defocus distance $\Delta z = - 3\;{\rm mm} $. The used intensity matrix is set as $256 \times 256$ for all measurement calculations.

In optical testing experiments, the first 231 terms of the ENZ function basis ${\boldsymbol C}_k^{{nm}}$ are applied to the ENZ-NOPR algorithm. To provide a reference phase map for the performance estimation, we first used the ZYGO interferometer to measure two samples, as shown in Figs. 5(a1) and 5(b1). The wavefront measured via the ENZ-NOPR method is fitted with 36 terms in the Zernike polynomial for sample 1 and 190 terms in the Zernike polynomial for sample 2. The results are shown in Figs. 5(a2) and 5(b2), respectively. After binning $2 \times 2\;{\rm pixels}$, ${Q} = 0.467$ for testing system is obtained, and the retrieved results are shown in Figs. 5(a3) and 5(b3). After two interferometric data are fitted with equivalent Zernike polynomial terms, the RMSE differences between the various measurements are calculated for different algorithms and are listed in Table 2. The best measured result via the ENZ-NOPR algorithm agrees with interferometric data to $0.0168\lambda$, RMSE.

Moreover, if enough ENZ basis functions are used in wavefront testing, the ENZ-NOPR algorithm has very high accuracy, theoretically. However, in practice, because of the influence of noise on complex amplitude, as shown in results of two experiments, the actual measurement accuracy is lower than the numerical experiment accuracy. Moreover, with the rapid development of high SNR cameras and lower coherent noise of laser illumination, the higher accuracy measurement for wavefront is possible to achieve using the ENZ-NOPR algorithm.

In summary, we have proposed a modal-based NOPR algorithm based on the ENZ theory. Using a common phase retrieval experimental arrangement, the full complex pupil wavefront over the relatively large NA system is tested with under-sampled intensity data. Compared with the conventional phase retrieval method, our proposed method for determining the phase in the measurement planes and then calculating the complex field in the desired plane has less experimental and computational complexity. Compared with the data measured by the interferometer, the reconstruction accuracy of the simple wavefront is close to $\lambda /60$, RMSE, whereas that of the general wavefront is better than $\lambda /30$, RMSE. Considering that the ENZ-NOPR breaks through the restriction of sampling rate, this proposed method would be a very effective wavefront measurement approach for large NA elements such as the freeform lens.

## Funding

Science Challenge Project (TZ2016006-0502-02); Laboratory of Precision Manufacturing Technology of CAEP (ZD18005).

## Acknowledgment

Author Lei Zhao expresses thanks to Dongming Sun for valuable suggestions.

## Disclosures

The authors declare no conflicts of interest.

## REFERENCES

**1. **J. R. Fienup, Appl. Opt. **21**, 2758 (1982). [CrossRef]

**2. **J. W. Goodman, *Introduction to Fourier Optics*, 4th ed. (WH Freeman, 2017).

**3. **G. R. Brady, M. Guizar-Sicairos, and J. R. Fienup, Opt. Express **17**, 624 (2009). [CrossRef]

**4. **A. M. Michalko and J. R. Fienup, Opt. Lett. **43**, 1331 (2018). [CrossRef]

**5. **A. M. Michalko and J. R. Fienup, Opt. Lett. **43**, 4827 (2018). [CrossRef]

**6. **J. R. Fienup, J. Opt. Soc. Am. A **16**, 1831 (1999). [CrossRef]

**7. **J. S. Smith, D. L. Aronstein, B. H. Dean, and D. S. Acton, Proc. SPIE **7436**, 7436D (2009). [CrossRef]

**8. **S. T. Thurman and J. R. Fienup, J. Opt. Soc. Am. A **26**, 2640 (2009). [CrossRef]

**9. **G. R. Brady and J. R. Fienup, Appl. Opt. **48**, 442 (2009). [CrossRef]

**10. **E. A. Shields, Opt. Lett. **37**, 2463 (2012). [CrossRef]

**11. **A. J. E. M. Janssen, J. Opt. Soc. Am. A **19**, 849 (2002). [CrossRef]

**12. **J. Braat, P. Dirksen, and A. J. E. M. Janssen, J. Opt. Soc. Am. A **19**, 858 (2002). [CrossRef]

**13. **P. Dirksen, J. Braat, A. J. E. M. Janssen, and C. Juffermans,J. Micro/Nanolith. MEMS MOEMS **2**, 61 (2003). [CrossRef]

**14. **C. Van der Avoort, J. J. M. Braat, P. Dirksen, and A. J. E. M. Janssen, J. Mod. Opt. **52**, 1695 (2005). [CrossRef]

**15. **M. Guizar-Sicairos and J. R. Fienup, Opt. Express **16**, 7264 (2008). [CrossRef]

**16. **R. J. Noll, J. Opt. Soc. Am. **66**, 207 (1976). [CrossRef]

**17. **A. J. E. M. Janssen, J. Braat, and P. Dirksen, J. Mod. Opt. **51**, 687 (2004). [CrossRef]

**18. **L. Zhao, H. Yan, J. Bai, J. Hou, Y. He, X. Zhou, and K. Wang, Opt. Express **28**, 19726 (2020). [CrossRef]

**19. **W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *Numerical Recipes: The Art of Scientific Computing* (Cambridge University, 1986), Chap. 10.

**20. **J. R. Fienup, Appl. Opt. **32**, 1737 (1993). [CrossRef]

**21. **J. C. A. Barata and M. S. Hussein, Braz. J. Phys. **42**, 146 (2012). [CrossRef]