Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Non-propagation fast phase diverse phase retrieval for wavefront measurement with generalized FFT-based basis function

Open Access Open Access

Abstract

Phase retrieval is an attractive optical testing method with a simple experimental arrangement. The sampling grids wave propagation computation based on the FFT operations is usually involved in each iterative process for the classical phase retrieval model. In this paper, a novel non-propagation optimization phase retrieval technique with the FFT-based basis function is proposed to accelerate wavefront measurement. The sampling grids wave diffraction propagation computation is converted to matrix-vector products that have small dimensions to reduce the computational burden. The diffraction basis function based on generalized numerical orthogonal polynomial and two-step Fresnel propagation is deduced, which is suitable for the generally shaped pupil. This paper provides a universal non-propagation framework to accelerate phase retrieval which is applicable to the arbitrarily shaped wavefront measurement.

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

1. Introduction

Phase retrieval (PR), which retrieves the unknown phase of the desired pupil from measured intensities, is an image-based wavefront-sensing technique. Phase retrieval is an attractive method suitable for in situ wavefront measurement due to its relatively simple and flexible experimental arrangement [1]. For the most mainstream Gerchberg-Saxton-Fienup algorithms [2,3] applied in wavefront measurement, a pair of propagation operations typically involving fast Fourier Transforms (FFTs) is always involved. Phase diverse phase retrieval (PDPR) [46], which has an iterative calculation process among different measured planes based on FFT operations, has higher measurement robustness and accuracy. In addition, the nonlinear optimization theory combined with the iterative phase retrieval method to estimate the desired complex field offers more flexibility. The nonlinear optimization phase retrieval method has several advantages [7]: firstly, the phase retrieval process can be implemented based on a more modern nonconvex optimization algorithm, which is beneficial for eliminating the local minimum stagnation; secondly, the convergence speed and robustness of the PR algorithm would be improved by using adaptive step-size [8]; thirdly, the weight function, which is applied to discard the poor signal-to-noise ratio (SNR) pixels for the measured intensity, is added into the error metric to improve the robustness of the PR algorithm.

Generally, the nonlinear optimization phase retrieval algorithm [7] includes the point-by-point optimization algorithm and the basis-function optimization algorithm. The basis-function optimization algorithm is also called the parametric PDPR algorithm [9]. For classical parametric PDPR algorithm, the phase of desired plane is modeled by using the Zernike polynomials, respectively, then some dozens of Zernike coefficients are optimized. Compared with point-by-point optimization, where there may be more than 105 or 106 parameters over which we are optimizing, the parametric PDPR algorithm has less computationally demanding. The parametric PDPR algorithm has excellent performance for measuring large-scale aberration [9]. The previous nonlinear parametric optimization phase retrieval algorithms have higher accuracy for wavefront measurement, which is suitable for the circular pupil with Zernike polynomials [1013].

However, while the PDPR measurement method has a simple experimental arrangement, it would demand many computational resources [14]. The essential FFT operations, which need adequate intensity sampling and match testing system through zero-padding, would cost much time and memory. The large computational burden would seriously slow down wavefront measurement efficiency [15]. The parametric PDPR algorithm, which still needs to calculate polynomial basis functions, would further improve the requirements of computer performance.

Besides, in practical engineering applications, many optical elements with irregular apertures are widely applied in various large optical systems. For example, the large aperture telescope such as the James Webb space telescope (JWST) [16] is made up of several hexagonal mirrors; high-energy laser systems such as inertial confinement fusion (ICF) [17] contain many square optical elements. The conventional parametric PDPR cannot be applied to measure irregularly shaped wavefront. To propose a new parametric PDPR method, which has the characteristics of high speed and low computational complexity, it is essential to achieve wavefront measurement for the arbitrarily shaped pupil.

In this paper, a new non-propagation nonlinear optimization phase diverse phase retrieval method suitable for the generally shaped pupil (GSP) wavefront measurement is proposed to reduce the computational complexity and costing time. Instead of using sampling grid, the complex field of the GSP is modeled via numerical orthogonal polynomials [18] and the diffraction integral is also modeled using the FFT-based basis function. Because the desired field and diffractive field are both modeled via the basis function and the complex coefficients of the basis functions are estimated, the proposed parametric PDPR algorithm is called modal-based PDPR algorithm (mPDPR). The proposed method, which is suitable for the arbitrarily shaped pupil’s reconstruction, is far faster than existing methods. To reduce the computational complexity and costing time, these steps for accelerating the mPDPR method are implemented: firstly, for each diffraction basis matrix, only a small dimensions grid that contains valuable information is extracted to estimate the gradient matrix, which is under-sampled for the classical iterative phase retrieval method; secondly, the gradient calculation kernel matrices are constant without repeating calculation for all iterations; in addition, the analytic gradient in the measured plane is derived by only matrix operation.

This paper is organized as follows. Section 2 introduces the FFT-based diffraction basis function calculation model based on numerical orthonormal polynomials. Section 3 describes the proposed non-propagation phase retrieval algorithm. Section 4 verifies the accuracy and stability of the method through numerical simulations. Section 5 presents some verified experiments. Section 6 concludes the paper.

2. FFT-based diffraction basis function calculation

For the classical phase retrieval algorithm, the two-dimensional Fourier transform, which is implemented by using the FFT algorithm, is very commonly used to model the numerical diffractive propagations. For the proposed mPDPR method, the FFT-based basis function, which does not repeat the calculation for every iteration, instead of the sampling grids is applied to represent the wave propagation computation. Here we introduce two-step diffraction theory combining with polynomials to construct the FFT-based basis function. This model has more flexibility in selecting the grids both in the pupil and the intensity plane.

2.1 Diffraction calculation model

The FFT operator imposes a fixed relationship between the sample spacing in its input and output domains based on the integer transform length. Generally, the array length and sampling relationship can be controlled coarsely through zero-padding. For the optical propagation modeled by using the FFT algorithm, there is a fixed sampled relationship, which is determined using a calculated sampling parameter, $\; Q$, given by

$$Q = \frac{{\lambda {F / \# }}}{{{d_u}}},$$
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 sampled. In our previous work [12], the under-sampled problem is discussed.

However, for the optical testing system with large Q, this would suffer from several disadvantages: firstly, the zero-padding for FFT would waste a lot of computing resources; secondly, the aliasing of discrete circular apertures would undermine the accuracy for calculating diffraction from geometrical shapes [19]; thirdly, for optical propagation represented by FFT operator, Q is not a continuous quantity which could undermine the matching accuracy between the numerical model and the real testing system [20].

In our model, the two-step Fresnel operator [21], which achieves the precise matching between the numerical propagation model and the real measurement system by establishing an intermediate plane, is applied to diffraction calculation. The addition of defocus [22] can improve performance of the phase retrieval method, so the intensity is collected at a properly defocus position. For the wavefront measurement system, the desired field $P({x,y} )$ is numerically propagated to the ${k^{\textrm{th}}}$ measured plane ${U_k}({u,v} )$

$$\begin{aligned} {U_k}(u,v) &= \textrm{exp} \left[ { - i\frac{{2\pi }}{\lambda }\frac{{1 - m}}{{m{z_k}}}({{u^2} + {v^2}} )} \right] \times \\ &\textrm{ }{{\cal F}^{ - 1}}\left\{ \begin{array}{l} \textrm{exp} \left[ {i\frac{{{z_k}\pi \lambda }}{{m{{({D{d_1}} )}^2}}}({{x^2} + {y^2}} )} \right]\\ \times {\cal F}\left\{ {P({x,y} )\textrm{exp} \left( { - i\frac{{2\pi }}{\lambda }\frac{{{x^2} + {y^2}}}{{2f}}} \right)\textrm{exp} \left[ {i\frac{{2\pi }}{\lambda }\frac{{1 - m}}{{2{z_k}}}({{x^2} + {y^2}} )} \right]} \right\} \end{array} \right\}, \end{aligned}$$
where ${i^2} ={-} 1,$ ${\cal F}$ and ${{\cal F}^{ - 1}}$ represent FFT and inverse FFT, $\; f$ is focal length, ${d_1}$ and ${d_2}$ are the pixel pitch of desired pupil and detector, respectively, $m = {d_2}/{d_1}$. For convenience, Eq. (2) is abbreviated as
$${U_k}({u,v} )= 2{{\cal F}}\{{P({x,y} ),{z_k}} \}.$$

2.2 FFT-based diffraction basis function calculation for phase retrieval

Here, the diffraction basis function computation method to replace the sampling grid is proposed, which is not only suitable for various shaped pupils but also can be applicable to representing the complex field (including phase and amplitude) based on numerical orthonormal polynomials and the two-step Fresnel propagation model. The orthogonal property of polynomial implies that no redundant information overlapped between the coefficients of the polynomial, which would reduce the amount of the polynomials for retrieving the wavefront [23]. Therefore, in this paper, the numerical orthogonal polynomials are selected to represent the complex field of the desired plane. The under-tested wavefront $P({x,y} )$ can be modeled via polynomials $\{{{Z^j}({x,y} )} \}$

$$\begin{aligned} P({x,y} )&= {P_R}({x,y} )+ i{P_I}({x,y} )\\ &\textrm{ } = \sum\limits_{j = 1}^{j = J} {({\alpha_R^j + i\alpha_I^j} ){Z^j}({x,y} )} \\ &\textrm{ } = \sum\limits_{j = 1}^{j = J} {{\beta ^j}{Z^j}({x,y} )} , \end{aligned}$$
where $\; {P_R}({x,y} )$ and ${P_I}({x,y} )$ represent the real part and imaginary part of the under-tested wavefront, respectively. $\alpha _R^j$ and $\alpha _I^j$ represent the real part and imaginary part of fitted coefficients, respectively. Finally, ${\beta ^j} = \alpha _R^j$ + $\alpha _I^j$ is a complex value.

In the numerical experiments, when the complex field is represented by polynomials, the fitting accuracy of amplitude is far lower than that of phase, as shown in Fig. 1. The representation error of amplitude directly deduces the inaccurate amplitude reconstruction for phase retrieval. The reason for the representation error is analyzed at first. For the desired wavefront, the amplitude $|{P({x,y} )} |$ and phase $\theta ({x,y} )$ are calculated with

$$\left\{ {\begin{array}{{c}} {|{P({x,y} )} |= \sqrt {{{({{P_R}} )}^2} + {{({{P_I}} )}^2}} ,}\\ {\theta ({x,y} )= \arctan ({{{{P_I}} / {{P_R}}}} ).} \end{array}} \right.$$

From Eq. (5), we can know that the amplitude is the root of the sum of squares of real and imaginary parts for the complex field, while the phase is the $arctan$ value of the ratio of the real and imaginary parts for the complex field. For a complex value grid, the polynomials have equal representation ability for the real part ${P_R}({x,y} )$ and the imagery part ${P_I}({x,y} )$ of the complex wavefront. When the imaginary part and the real part are not fitted properly, it has the greatest influence on the amplitude, but almost no influence on the phase. For example, if the weight factor which represents the fitting accuracy is 0.99, phase is calculated with $\theta = \textrm{arctan}({0.99{P_I}/0.99{P_R}} )= \textrm{arctan}({{P_I}/{P_R}} )$, while the amplitude is $|P |= \sqrt {{{({0.99{P_R}} )}^2} + {{({0.99{P_I}} )}^2}} = 0.99\sqrt {{{({{P_R}} )}^2} + {{({{P_I}} )}^2}} $. It is equivalent to multiplying by a weight factor on the imaginary part and the real part, which is canceled when the phase is calculated. Therefore, the representation accuracy of the phase would hardly be affected. The numerical simulation is shown in Fig. 1. Fortunately, the surface of the optical element is only related to the phase but not to the amplitude, so the fitting error would hardly affect the surface measurement accuracy.

 figure: Fig. 1.

Fig. 1. The fitting results for complex value via polynomials. (a1) and (a2) are fitted phase and amplitude, respectively; (b1) and (b2) are truth phase and amplitude, respectively. The residual RMSE is 0.0081 $\mathrm{\lambda }$ for phase, while that is 0.0817 for amplitude.

Download Full Size | PDF

Here, Zernike polynomials $\{{{Z^j}} \}$ [24] are utilized to determine orthonormal polynomials $\{{{{\hat{Z}}^j}} \}$ suitable for the generally shaped pupil. It is worth noting that other polynomials such as Q-polynomials [25], Legendre polynomials and Chebyshev polynomials [26] can also be used to calculate the orthonormal polynomials. The wavefront $P({x,y} )$ can be modeled via numerical orthonormal polynomials as

$$P({x,y} )= \sum\limits_{j = 1}^J {{\beta ^j}{{\widehat Z}^j}} ({x,y} ),$$
where J is the number of polynomials, $\{{{\beta^j}} \}$ are the complex value coefficients. According to Fourier optics theory [27], the complex amplitudes can be linearly superimposed in coherent optical systems. Therefore, the diffraction integral $\{{{u^j}} \}$ of $\{{{{\hat{Z}}^j}} \}$ superposed with weight $\{{{\beta^j}} \}$ is identical to the diffraction integral ${U_k}({u,v} )$ of $P({x,y} )$. The ${U_k}({u,v} )$ can be expressed as
$$\begin{aligned} {\boldsymbol{U}_k} &= \left\{ {\begin{array}{{c}} \vdots \\ {{{[{P({2{\cal F}({\widehat {{Z_k}}({x,y} ),{z_k}} )} )} ]}^H}}\\ \vdots \end{array}} \right\}\boldsymbol{\beta } = \left\{ {\begin{array}{{c}} \vdots \\ {{{[{P({{u^j}} )} ]}^H}}\\ \vdots \end{array}} \right\}\boldsymbol{\beta },\\ &\textrm{ } = {\boldsymbol{C}_k}\boldsymbol{\beta }, \end{aligned}$$
where ${{\boldsymbol U}_k}$ is the one-dimensional vector of ${U_k}({u,v} ),$ ${\boldsymbol \beta }$ is the set of $\{{{\beta^j}} \}$, $P({} )$ denotes extracted valuable basis function grid, ${[{} ]^H}$ denotes transformation to a one-dimensional vector and ${{\boldsymbol C}_k}$ represents the set of diffraction basis functions.

 figure: Fig. 2.

Fig. 2. Experimental configuration for wavefront measurement. The collimated beam emitted from a ZYGO interferometer illuminates the plate (sample). After focusing the beam, the diffraction intensities are collected at different defocus positions. The intensity image in the blue box denotes the extracted valuable grid. Every diffraction basis function is also extracted valuable grid, correspondingly.

Download Full Size | PDF

3. Non-propagation phase diverse phase retrieval algorithm

After calculating the FFT-based basis function, the mPDPR, which is the nonlinear optimization phase diverse phase retrieval without wave propagation iteration, is performed. The optical testing system is shown in Fig. 2. The nonlinear optimization process can be expressed for optimizing error metric ${E_k}$ by estimating ${\boldsymbol \beta }$

$$\textrm { find } \boldsymbol{\beta} \in \min \sum_{k}\left(E_{k}:\left(\sum \boldsymbol{S}_{k} \left| \widehat{\boldsymbol{I}_{k}}-|\widehat{\boldsymbol{C}_{k} \boldsymbol{\beta}}| \right|^{2}\right)^{1 / 2}\right),$$
where ${{\boldsymbol S}_k}$ is the weighting function used to discard the effect of bad or saturated detector pixels and the pixels with the poor signal-to-noise ratio (SNR) for the ${k^{\textrm{th}}}$ measured intensity.

For classical gradient-search model-based optimization algorithm, a pair of propagation operators are used to calculate gradient [3,7]. The basic gradient optimization method is used to test the proposed algorithm, the performance of the algorithm would be further improved via the advanced gradient optimization method [28]. The forward propagation can be represented via Eq. (7), while the inverse forward propagation is not expressed. The coefficient gradient is supposed as $\partial E_k/\partial \boldsymbol{\beta }$, according to Eq. (3) and Eq. (31) in Ref. 7, $\partial {E_k}/\partial \boldsymbol{\beta }\ast {{\boldsymbol C}_k} = \left( {|{{{\boldsymbol C}_k}{\boldsymbol \beta }} |- \sqrt {\widehat {{{\boldsymbol I}_k}}} } \right)\textrm{arg}({{{\boldsymbol C}_k}{\boldsymbol \beta }} )$ (“ $\ast$ ” denotes forking product), the gradient is calculated as

$$\frac{{\partial {E_k}}}{{\partial \boldsymbol{\beta} }} = \boldsymbol{T}_k^{}\left( {\widehat {|{{\boldsymbol{C}_k}\boldsymbol{\beta} } |} - \sqrt {\widehat {{\boldsymbol{I}_k}}} } \right)\arg ({{\boldsymbol{C}_k}\boldsymbol{\beta} } ).$$

Considering that ${{\boldsymbol C}_k}$ may not be a square matrix, here, the pseudo-inverse matrix ${{\boldsymbol T}_k}$ of ${{\boldsymbol C}_k}$ is calculated with

$${\boldsymbol{T}_k} = {({\boldsymbol{C}_k^\textrm{T}{\boldsymbol{C}_k}} )^{ - 1}}\boldsymbol{C}_k^\textrm{T}$$

As shown in Eq. (9), for every gradient solution, ${{\boldsymbol T}_k}$ and ${{\boldsymbol C}_k}$ are constant at the distance ${z_k}$ without repeating calculation. Therefore, there is only a matrix formulation without a pair of propagation operations. The computational complexity is sharply decreased.

The proposed mPDPR algorithm is summarized in Fig. 3. Before the mPDPR is executed to measure wavefront error, the preparatory work needs to be done:

  • (1) According to the shape of the pupil, numerical orthonormal polynomials $\{{{{\hat{Z}}^j}} \}$ are calculated;
  • (2) Diffraction integral ${\{{{u^j}} \}_k}$ is estimated for each polynomial $\{{{{\hat{Z}}^j}} \}$. The valuable matrix $256 \times 256$ is extracted from the diffraction integral with $512 \times 512$ dimensions.
  • (3) After transforming $u_k^j$ into a one-dimensional vector, the basis function matrix ${{\boldsymbol C}_k}$ and gradient kernel matrix ${{\boldsymbol T}_k}$ are calculated via the diffraction integral ${\{{{u^j}} \}_k}$.

 figure: Fig. 3.

Fig. 3. Flow chart of the reconstruction process. The numerical orthogonal polynomials suitable for the irregular pupil is calculated at first. Then the FFT-based basis function is calculated via two-step Fresnel diffraction. The valuable dimensions diffraction basis function grid is extracted, and the ${{\boldsymbol C}_k}$ and ${{\boldsymbol T}_k}$ is calculated. Using Eq. (9), we can calculate the gradient of polynomial coefficient. For the mPDPR method, the iterative sequence is unordered (such as plane: $1 \to 3 \to 2 \to 1 \to 2 \to 3 \cdots $) to improve the speed of convergence. Finally, the wavefront is fitted by using estimated coefficients $\; {\boldsymbol \beta }$ with polynomials $\{{{{\hat{Z}}^j}} \}.$

Download Full Size | PDF

The box on the right of Fig. 3 lists the specific calculation process for mPDPR, and the nth iteration can be carried out as follows:

  • (1) Calculate corresponding gradient function ${{\boldsymbol D}_n}$ for ${k^{th}}$ error metric ${E_k}$ using Eq. (9).
  • (2) Update estimated coefficients with
    $$\boldsymbol{\beta}_{n}=\boldsymbol{\beta}_{n-1}+h_{n} \boldsymbol{D}_{n},$$
    where ${h_n}$ is the step size.
  • (3) The new estimated wavefront in the ${k^{th}}$ plane is specified by
    $${\boldsymbol{U}_{n + 1,k}} = {\boldsymbol{C}_k}{\boldsymbol{\beta} _n}.$$

During the whole iterative process, the applied sequence of the intensity plane is unordered to provide better nonredundant information for intensity distributions. The convergence speed of the algorithm is improved via the ordered iterative phase retrieval methods. The procedure is iteratively repeated until the “stopping criterion” is satisfied.

4. Simulations

4.1 Applicability of mPDPR method for various shaped pupils

First, the numerical simulations for various shaped wavefront reconstructions are done to prove the applicability of the mPDPR method. The first 36 terms Q-polynomials with weighted coefficients are composed to construct wavefront with differently shaped mask. The first 400 numerical orthogonal polynomials for mPDPR are used to retrieve the complex wavefront. The relationship between the retrieved residual RMSE and the number of polynomials is analyzed. The retrieved results are shown in Fig. 4.

 figure: Fig. 4.

Fig. 4. The retrieved wavefront results (including phase and amplitude) with differently shaped pupils. (a1) ∼ (d1) are the wavefront with hexagon mask; (a2) ∼ (d2) are the wavefront with square mask; (a3) ∼ (d3) are the wavefront with square aperture mask. The number in the lower right corner of every subfigure is the residual RMSE value.

Download Full Size | PDF

As shown in Fig. 4, the proposed mPDPR algorithm is suitable for various shaped wavefront reconstructions. Besides, the phase, which is directly related to the optical surface error, is in good agreement with the ground truth, while the accuracy of reconstructed amplitude is lower than that of the phase. It is also proved that the inaccuracy of amplitude fitting would hardly affect the accuracy of phase (surface error) measurement. Next, the number of orthogonal polynomials that is actually necessary for the wavefront reconstruction is analyzed by establishing the relationship between the residual RMSE and the number of polynomials. The SNRs for all diffraction intensities are 40dB. The RMSE of the reconstructed wavefront is not equal to the difference between it and the real value due to the influence of noise. As shown in Fig. 5, the RMSE curve has fluctuated when the number of used polynomials is fewer. The RMSE is tended to stabilize when the number of orthogonal polynomials is about 400.

 figure: Fig. 5.

Fig. 5. The composite plot of RMSE for different numbers of polynomials. (a) is the RMSE curve for hexagon wavefront reconstruction with different numbers of polynomials; (b) is the RMSE curve for square wavefront reconstruction with different numbers of polynomials; (c) is the RMSE curve for square aperture wavefront reconstruction with different numbers of polynomials.

Download Full Size | PDF

4.2 Comparisons of mPDPR and conventional algorithms

To demonstrate its convergence speed, the proposed mPDPR method is numerically tested for the hexagon wavefront. The aberrations represented by the first 36 terms Q-type polynomials of optical elements is the ground truth. All experimental data and code are implemented on MATLAB 2017b using a desktop PC with an NVIDIA GeForce GTX 1060. The first 160 Zernike polynomials which are orthogonalized are applied to calculate the diffraction basis functions. The SNRs for all diffraction intensities are 30dB. The retrieved surface (relating to the phase) results via mPDPR and modified multiple-image Gerchberg-Saxton (Multi-GS) algorithms [29] are shown in Fig. 6. The mPDPR algorithm can retrieve the phase precisely, which well agrees with the truth values. Considering that the surface of the under-tested object is directly related to the phase of the wavefront rather than the amplitude, there is only a comparison between the retrieved values and truth values as shown in Fig. 7.

 figure: Fig. 6.

Fig. 6. Numerical demonstration wavefront measurement using mPDPR and Multi-GS algorithm.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. The composite plot of phase RMSE for different algorithms. We can see that the mPDPR has a far faster convergence speed and the accuracy of mPDPR for 300 iterations (RMSE 0.0098 $\mathrm{\lambda }$) is higher than Multi-GS for 1000 iterations (0.0123 $\mathrm{\lambda }$).

Download Full Size | PDF

As shown in Fig. 7, we can know that the far faster convergence speed for the wavefront reconstruction is realized by using the mPDPR method. The costing time of 500 times iteration for the mPDPR is only 1.78s, while Multi-GS needs 25.37s. Besides, when the unordered iteration is introduced, the reconstruction speed is far faster than the traditional Multi-GS algorithm. The RMSE for mPDPR is 0.0098 $\mathrm{\lambda }$ when the iterative number is 300 times, while the Multi-GS is 0.0123 $\mathrm{\lambda }$ when the iterative number is 1000 times. Monte Carlo simulations are further performed to validate the effectiveness of the proposed method. Two different noise levels (25dB and 35dB) are considered, respectively. The residual RMSE and costing time for 500 iterations are compared with Multi-GS by using box-plot in Fig. 8. The compared results show that the proposed algorithm has comparable measurement accuracy with Multi-GS while the costing time is sharply decreased.

 figure: Fig. 8.

Fig. 8. Summary of 50 numerical experiments for wavefront reconstruction via different algorithms using box-plot. (a) shows RMSEs versus different algorithms. Suffix ‘−25’ represents that the SNR of diffraction intensity is 25 dB; suffix ‘−35’ represents that the SNR of diffraction intensity is 35 dB. (b) shows the costing times for different algorithms.

Download Full Size | PDF

5. Experiments

We further perform two real experiments to demonstrate the effectiveness of the proposed method. In these experiments, the wavelength $\mathrm{\lambda } = 632.8nm$, pixel size $4.4 \times 4.4\mu m$, and the focal length of convergence lens $1079.41mm$. The first sample is testing wavefront with aperture diameter $\; D = 22.9mm$ (the sampling parameter, $\; Q$, is 6.78). For intensity sampling $512 \times 512$, the effective pupil sampling is $75 \times 75$ for FFT-based propagation. In other words, if the effective pupil sampling is $512 \times 512$, the dimensions of the grid which is applied to numerical propagation are $3471 \times 3471$, which leads to a significant computational burden for the classical Multi-GS algorithm. The first 300 terms orthogonal Zernike polynomials are applied to calculate the diffraction basis functions, correspondingly. The Multi-GS and model-based phase retrieval with Zernike polynomials (Zernike-PR) [10] algorithms are introduced based on the two-step Fresnel propagation. The retrieved results are fitted with modeled using a superposition of 300 standard Zernike polynomials. The desired plane sampling is $512 \times 512$ for mPDPR and Zernike-PR algorithms. Figure 9 shows the reconstruction results with different methods.

 figure: Fig. 9.

Fig. 9. Experimental results for different phase retrieval algorithms. (a1) ∼ (c1) are retrieved phases via Multi-GS, Zernike-PR and mPDPR algorithms, respectively; (a2) ∼ (c2) are residual errors for (a1) ∼ (c1), correspondingly; (d) is the ZYGO testing result.

Download Full Size | PDF

The summary of RMSE and costing time for different algorithms are shown in Table 1. The accuracy of the proposed mPDPR algorithm is close to that of other algorithms. For costing time, the mPDPR, which only needs 2.71s, has an enormous advantage compared with other algorithms. If the Multi-GS with FFT operator rather than two-step Fresnel operator, its costing time is 18.37s. Then the number of orthogonal polynomials that is actually necessary for the circle wavefront reconstruction is analyzed by establishing the relationship between the residual RMSE and the number of polynomials. The residual RMSE for different numbers of polynomials is plotted in Fig. 10. When the number of polynomials is more than 280, the residual RMSE is hardly reduced. The convergence tendency of the RMSE is similar to the numerical experiment shown in Fig. 5.

 figure: Fig. 10.

Fig. 10. The curve of RMSE vs different the number of polynomials for circle wavefront. When the number of polynomials is more than 280, it is sufficient for wavefront reconstruction.

Download Full Size | PDF

Tables Icon

Table 1. Summary of RMSE and costing time for different algorithms

For the second experiment, the wavefront of the light emitted from the ZYGO interferometer is measured. A square mask with aperture diameter $\; D = 70mm$ is placed near the interferometer. The retrieved wavefronts via Multi-GS and mPDPR algorithms are shown in Figs. 11(a) and 11(b). The first 200 terms of Zernike polynomials are orthogonalized to calculate diffraction basis functions. The first 36 estimated coefficients of numerical orthogonalized polynomials are compared in Fig. 11(c). This experiment shows the effectiveness of the mPDPR algorithm for irregular pupil wavefront reconstruction. The costing time for 500 times iteration for Multi-GS is 24.21s, while the mPDPR only needs 1.84s. Then the number of orthogonal polynomials that is actually necessary for the square wavefront reconstruction is analyzed by establishing the relationship between the residual RMSE and the number of polynomials. The residual RMSE for different numbers of polynomials is plotted in Fig. 12. When the number of polynomials is about 200, the residual RMSE is tended to stabilize.

 figure: Fig. 11.

Fig. 11. Experimental results for different phase retrieval algorithms. (a) is retrieved phase via Multi-GS algorithm. (b) is retrieved phase via mPDPR. (c) The composite plot of the coefficient value of wavefront for different algorithms.

Download Full Size | PDF

 figure: Fig. 12.

Fig. 12. The curve of RMS errors vs different the number of polynomials for square wavefront. When the number of polynomials is about 200, it is sufficient for wavefront reconstruction.

Download Full Size | PDF

6. Discussion and conclusion

In summary, this paper establishes a framework to greatly improve the efficiency of phase retrieval wavefront measurement by converting FFT-based diffraction propagation to matrix operation. The main novelty of the proposed method is summarized as follows: firstly, the proposed method can be applied to reconstruct the arbitrarily shaped pupils; secondly, the proposed method is far faster than existing methods. We realize the speed-up of the proposed algorithm through two aspects: firstly, with constant kernel matrix without repeating calculation for each iteration, using matrix operation instead of repeating a pair of wave propagations for every iteration greatly improve the measurement speed; secondly, small dimensions grid for gradient calculation matrix and diffraction intensity is achieved via solving the coefficients of polynomial instead of point-by-point recovery. Two-step Fresnel propagation offers a more flexible sampling for the desired pupil. The measured accuracy of the proposed method is comparable to that of the traditional method, while the time costing and memory costing have been greatly reduced. It is worth mentioning that the proposed method is completely suitable for wavefront detection. To improve the robustness and accuracy of the algorithm, three intensity patterns are used to realize wavefront measurement. This method is adaptive to arbitrary pupil wavefront retrieval based on numerical orthogonal polynomials and especially suitable for irregular shape optical element wavefront and freeform surface measurement.

Besides, we proposed a new PDPR method to measure wavefront, there are still some interesting issues to think further. In this paper, the fixed polynomial terms are chosen to represent the wavefront errors for aberrations measurement. We discuss the influence of the number of polynomials on the measurement accuracy for all experimental data in this paper. For the measurement of low-frequency aberrations, the proposed method with those number of diffraction basis functions is sufficient. When the proposed algorithm is further applied to the mid- and high- spatial frequency wavefront error measurement, the relationship between the number of diffraction basis functions and the spatial frequency band of wavefront needs to be further explored. In the future, the compressed sensing [30] approaches would be explored to improve the flexibility of the algorithm for choosing the number of polynomial basis elements.

Funding

Science Challenge Project (TZ2016006-0502-02); National Natural Science Foundation of China (52075507, 61905241); Laboratory of Precision Manufacturing Technology of CAEP (ZD18005).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. F. Jiang, G. Ju, X. Qi, and S. Xu, “Cross-iteration deconvolution strategy for differential optical transfer function (dOTF) wavefront sensing,” Opt. Lett. 44(17), 4283–4286 (2019). [CrossRef]  

2. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237 (1972).

3. J. R. Fienup, “Phase Retrieval Algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef]  

4. G. R. Brady and J. R. Fienup, “Nonlinear optimization algorithm for retrieving the full complex pupil function,” Opt. Express 14(2), 474–486 (2006). [CrossRef]  

5. G. Ju, C. Yan, D. Yue, and Z. Gu, “Field diversity phase retrieval method for wavefront sensing in monolithic mirror space telescopes,” Appl. Opt. 56(15), 4224–4237 (2017). [CrossRef]  

6. J. Fan, Q. Wu, Z. Lu, X. Li, and B. Chen, “Cramér–Rao lower bound analysis of phase diversity for sparse aperture optical systems,” Appl. Opt. 56(9), 2563–2567 (2017). [CrossRef]  

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

8. C. Zuo, J. Sun, and Q. Chen, “Adaptive step-size strategy for noise-robust Fourier ptychographic microscopy,” Opt. Express 24(18), 20724 (2016). [CrossRef]  

9. J. E. Krist and C. J. Burrows, “Phase-retrieval analysis of pre-and post-repair Hubble Space Telescope images,” Appl. Opt. 34(22), 4951–4964 (1995). [CrossRef]  

10. L. Zhao, H. Yan, J. Bai, J. Hou, Y. He, X. Zhou, and K. Wang, “Simultaneous reconstruction of phase and amplitude for wavefront measurements based on nonlinear optimization algorithms,” Opt. Express 28(13), 19726 (2020). [CrossRef]  

11. G. R. Brady, “Application of phase retrieval to the measurement of optical surfaces and wavefronts,” Ph.D. thesis (University of Rochester, 2008).

12. L. Zhao, J. Bai, Y. Hao, H. Jing, C. Wang, B. Lu, Y. Liang, and K. Wang, “Modal-based nonlinear optimization algorithm for wavefront measurement with under-sampled data,” Opt. Lett. 45(19), 5456–5459 (2020). [CrossRef]  

13. J. Antonello and M. Verhaegen, “Modal-based phase retrieval for adaptive optics,” J. Opt. Soc. Am. A 32(6), 1160–1170 (2015). [CrossRef]  

14. Q. Xin, G. Ju, C. Zhang, and S. Xu, “Object-independent image-based wavefront sensing approach using phase diversity images and deep learning,” Opt. Express 27(18), 26102–26119 (2019). [CrossRef]  

15. A. S. Jurling, M. D. Bergkoetter, and J. R. Fienup, “Techniques for arbitrary sampling in two-dimensional Fourier transforms,” J. Opt. Soc. Am. A 35(11), 1784–1796 (2018). [CrossRef]  

16. P. A. Sabelhaus and J. E. Decker, “An overview of the James Webb Space Telescope (JWST) project,” Proc. SPIE 5487, 550–563 (2004). [CrossRef]  

17. D. J. Trummer, R. J. Foley, and G. S. Shaw, “Stability of optical elements in the NIF target area building,” in Third International Conference on Solid State Lasers for Application to Inertial Confinement Fusion (International Society for Optics and Photonics, 1999), 363–371.

18. J. Ye, W. Wang, Z. Gao, Z. Liu, S. Wang, P. Benítez, J. C. Miñano, and Q. Yuan, “Modal wavefront estimation from its slopes by numerical orthogonal transformation method over general shaped aperture,” Opt. Express 23(20), 26208–26220 (2015). [CrossRef]  

19. S. D. Will and J. R. Fienup, “Algorithm for exact area-weighted antialiasing of discrete circular apertures,” J. Opt. Soc. Am. A 37(4), 688–696 (2020). [CrossRef]  

20. A. S. Jurling and J. R. Fienup, “Phase retrieval with unknown sampling factors via the two-dimensional chirp-Z transform,” J. Opt. Soc. Am. A 31(9), 1904–1911 (2014). [CrossRef]  

21. C. Rydberg and J. Bengtsson, “Efficient numerical representation of the optical field for the propagation of partially coherent radiation with a specified spatial and temporal coherence function,” J. Opt. Soc. Am. A 23(7), 1616 (2006). [CrossRef]  

22. Y. Zhang, G. Pedrini, W. Osten, and H. Tiziani, “Whole optical wave field reconstruction from double or multi in-line holograms by phase retrieval algorithm,” Opt. Express 11(24), 3234–3241 (2003). [CrossRef]  

23. C. Wei, Y. Li, W. Chau, and C. Li, “Trademark image retrieval using synthetic features for describing global shape and interior structure,” Pattern Recogn. 42(3), 386–394 (2009). [CrossRef]  

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

25. G. W. Forbes, “Fitting freeform shapes with orthogonal bases,” Opt. Express 21(16), 19061–19081 (2013). [CrossRef]  

26. W. Sun, L. Chen, W. Tuya, Y. He, and R. Zhu, “Analysis of the impacts of horizontal translation and scaling on wavefront approximation coefficients with rectangular pupils for Chebyshev and Legendre polynomials,” J. Opt. Soc. Am. A 30(12), 2539–2546 (2013). [CrossRef]  

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

28. J. Zhong, L. Tian, P. Varma, and L. Waller, “Nonlinear optimization algorithm for partially coherent phase retrieval and source recovery,” IEEE Trans. Comput. Imaging 2(3), 310–322 (2016). [CrossRef]  

29. J. A. Rodrigo, H. Duadi, T. Alieva, and Z. Zalevsky, “Multi-stage phase retrieval algorithm based upon the gyrator transform,” Opt. Express 18(2), 1510–1520 (2010). [CrossRef]  

30. M. Andersen, J. Dahl, and L. Vandenberghe, “CVXOPT: a Python package for convex optimization, version 1.1.7,” 2014, http://cvxopt.org.

Supplementary Material (1)

NameDescription
Supplement 1       supplemental document for time cost comparisons

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1.
Fig. 1. The fitting results for complex value via polynomials. (a1) and (a2) are fitted phase and amplitude, respectively; (b1) and (b2) are truth phase and amplitude, respectively. The residual RMSE is 0.0081 $\mathrm{\lambda }$ for phase, while that is 0.0817 for amplitude.
Fig. 2.
Fig. 2. Experimental configuration for wavefront measurement. The collimated beam emitted from a ZYGO interferometer illuminates the plate (sample). After focusing the beam, the diffraction intensities are collected at different defocus positions. The intensity image in the blue box denotes the extracted valuable grid. Every diffraction basis function is also extracted valuable grid, correspondingly.
Fig. 3.
Fig. 3. Flow chart of the reconstruction process. The numerical orthogonal polynomials suitable for the irregular pupil is calculated at first. Then the FFT-based basis function is calculated via two-step Fresnel diffraction. The valuable dimensions diffraction basis function grid is extracted, and the ${{\boldsymbol C}_k}$ and ${{\boldsymbol T}_k}$ is calculated. Using Eq. (9), we can calculate the gradient of polynomial coefficient. For the mPDPR method, the iterative sequence is unordered (such as plane: $1 \to 3 \to 2 \to 1 \to 2 \to 3 \cdots $ ) to improve the speed of convergence. Finally, the wavefront is fitted by using estimated coefficients $\; {\boldsymbol \beta }$ with polynomials $\{{{{\hat{Z}}^j}} \}.$
Fig. 4.
Fig. 4. The retrieved wavefront results (including phase and amplitude) with differently shaped pupils. (a1) ∼ (d1) are the wavefront with hexagon mask; (a2) ∼ (d2) are the wavefront with square mask; (a3) ∼ (d3) are the wavefront with square aperture mask. The number in the lower right corner of every subfigure is the residual RMSE value.
Fig. 5.
Fig. 5. The composite plot of RMSE for different numbers of polynomials. (a) is the RMSE curve for hexagon wavefront reconstruction with different numbers of polynomials; (b) is the RMSE curve for square wavefront reconstruction with different numbers of polynomials; (c) is the RMSE curve for square aperture wavefront reconstruction with different numbers of polynomials.
Fig. 6.
Fig. 6. Numerical demonstration wavefront measurement using mPDPR and Multi-GS algorithm.
Fig. 7.
Fig. 7. The composite plot of phase RMSE for different algorithms. We can see that the mPDPR has a far faster convergence speed and the accuracy of mPDPR for 300 iterations (RMSE 0.0098 $\mathrm{\lambda }$ ) is higher than Multi-GS for 1000 iterations (0.0123 $\mathrm{\lambda }$ ).
Fig. 8.
Fig. 8. Summary of 50 numerical experiments for wavefront reconstruction via different algorithms using box-plot. (a) shows RMSEs versus different algorithms. Suffix ‘−25’ represents that the SNR of diffraction intensity is 25 dB; suffix ‘−35’ represents that the SNR of diffraction intensity is 35 dB. (b) shows the costing times for different algorithms.
Fig. 9.
Fig. 9. Experimental results for different phase retrieval algorithms. (a1) ∼ (c1) are retrieved phases via Multi-GS, Zernike-PR and mPDPR algorithms, respectively; (a2) ∼ (c2) are residual errors for (a1) ∼ (c1), correspondingly; (d) is the ZYGO testing result.
Fig. 10.
Fig. 10. The curve of RMSE vs different the number of polynomials for circle wavefront. When the number of polynomials is more than 280, it is sufficient for wavefront reconstruction.
Fig. 11.
Fig. 11. Experimental results for different phase retrieval algorithms. (a) is retrieved phase via Multi-GS algorithm. (b) is retrieved phase via mPDPR. (c) The composite plot of the coefficient value of wavefront for different algorithms.
Fig. 12.
Fig. 12. The curve of RMS errors vs different the number of polynomials for square wavefront. When the number of polynomials is about 200, it is sufficient for wavefront reconstruction.

Tables (1)

Tables Icon

Table 1. Summary of RMSE and costing time for different algorithms

Equations (12)

Equations on this page are rendered with MathJax. Learn more.

Q = λ F / # d u ,
U k ( u , v ) = exp [ i 2 π λ 1 m m z k ( u 2 + v 2 ) ] ×   F 1 { exp [ i z k π λ m ( D d 1 ) 2 ( x 2 + y 2 ) ] × F { P ( x , y ) exp ( i 2 π λ x 2 + y 2 2 f ) exp [ i 2 π λ 1 m 2 z k ( x 2 + y 2 ) ] } } ,
U k ( u , v ) = 2 F { P ( x , y ) , z k } .
P ( x , y ) = P R ( x , y ) + i P I ( x , y )   = j = 1 j = J ( α R j + i α I j ) Z j ( x , y )   = j = 1 j = J β j Z j ( x , y ) ,
{ | P ( x , y ) | = ( P R ) 2 + ( P I ) 2 , θ ( x , y ) = arctan ( P I / P R ) .
P ( x , y ) = j = 1 J β j Z ^ j ( x , y ) ,
U k = { [ P ( 2 F ( Z k ^ ( x , y ) , z k ) ) ] H } β = { [ P ( u j ) ] H } β ,   = C k β ,
 find  β min k ( E k : ( S k | I k ^ | C k β ^ | | 2 ) 1 / 2 ) ,
E k β = T k ( | C k β | ^ I k ^ ) arg ( C k β ) .
T k = ( C k T C k ) 1 C k T
β n = β n 1 + h n D n ,
U n + 1 , k = C k β n .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.