Aberration retrieval by incorporating customized priors for estimating Zernike coefficients

Zernike expansion is an important tool for aberration retrieval in the optical field. The Zernike coefficients in the expansion can be solved in a linear system from those focal region intensity images, which can be modeled by the extended Nijboer–Zernike approach. Here we point out that those coefficients usually follow from different prior distributions, and especially, their variances could be dramatically diverse. To incorporate the prior information, we further introduce customized penalties to those Zernike coefficients and adopt a global adaptive generative adjustment algorithm for estimating coefficients. Based on both simulated and real data, numerical experiments show that our method outperforms other conventional methods, and provides an estimate of Zernike coefficients with a low mean square error.

www.nature.com/scientificreports/ distribution, which is Gaussian with mean zero and variance σ 2 n,m . The prior mentioned in our paper means a knowledge that the prior distributions for various model parameters could have diverse variances. Since we consider the Gaussian prior with zero mean, the prior distribution is determined only by its variance. The "customized priors" mean the customized variances for model parameters.
For incorporating the prior information, we further introduce customized penalties for those Zernike coefficients and adopt a global adaptive generative adjustment algorithm for estimating coefficients. The experiment results on the simulated and real data illustrate that our method, utilizing customized prior variances, provides an estimate with a lower mean squared error compared to other ENZ aberration retrieval methods.

Methods
Extended Nijboer-Zernike diffraction. In this study, we assume that the optical system is monochromatic and its aperture is circular and unobstructed. The generalized pupil function is usually expressed as 14 where A(ρ, θ) and φ(ρ, θ) are the amplitude and phase of the pupil plane, respectively. i = √ −1 and (ρ, θ) are polar coordinates, ρ ∈ [0, 1] , θ ∈ [0, 2π] . The phase can be linearly expressed by Frits Zernike expansion 14 where {Z k (ρ, θ)} K−1 k=0 are Frits Zernike basis functions, and {α} K−1 k=0 are the coefficients of the Zernike basis. From the work 15 , the generalized pupil function can be decomposed by using Zernike coefficients {β m n } as where m and n are integers such that n ≥ 0 and n − |m| is even. And R |m| n (ρ) is the radial polynomial. The first Zernike coefficient β 0 0 is real, and the others {β m n } are complex for m = 0, n = 0 . In the simulation, we generated random phases {α} K−1 k=0 using the method 13 . Thus we also obtained the phase function φ(ρ, θ) by (3). Furthermore,  www.nature.com/scientificreports/ we computed the generalized pupil function by (2) at constant amplitude. Finally, we consider the first 91 terms in the Zernike expansion (4) of the generalized pupil function.We denote all the Zernike coefficients by a coefficient vector β . Thus the pupil function is a linear function on the vector β , which is to be estimated as aberration retrieval. From works 11,15 , the radial polynomials in (4) can be expressed as Hence R |m| n (ρ) can be constructed before aberration retrieval for a specific choice of m, n, and ρ. In the work 8 , the light field on the focal plane is expressed as where the parameter f is the camera distance from the focus plane, and (r, ϕ) are polar coordinates on the image plane. From the work 7 , the Bessel series has the form where J m is a Bessel function of the first kind with order m, and where q = n+|m| 2 , p = n−|m| 2 , and l = 1, 2, . . . ; j = 0, . . . , p.
Aberration retrieval model using the ENZ approach. Using the formula (1), the PSF intensity in the focal region can be expressed as where is the sum of the remaining second order cross terms. In the summand, symbol ′ means the omission of n = 0 terms in the summand, and symbol ′′ means the omission of n 1 = m 1 = 0 or n 2 = m 2 = 0 terms. Re() and Im() denote the real and imaginary parts of a complex number. Symbol * denote the complex conjugate. The PSF intensity with an additive detector readout noise ε forms an aberration retrieval model: where the noise ε(r, ϕ; f ) ∼ N(0, σ 2 ) . The intensity I b (r, ϕ; f ) can be collected to estimate coefficients {β m n } using (8), as shown in Fig. 3.
(2) Assumes that I (k) can be described as linear combinations of the entrance pupil aberrations with coefficients {β m n } . This is equivalent to omitting the cross terms of (8). And then compute {β m n } (k) .
Notice that due to phase wrapping effects occurring in the reconstructed pupil distribution, this retrieval process may fail in case that the aberration magnitude is beyond some large range. So we only considered a small www.nature.com/scientificreports/ wavefront error (pv < 2π) in our experiment. From works 10,12,16,17 , we know that the intensity I (k) in Step (2) can be linearly transformed into a vector L (k) . We further get a classical linear model on the coefficient vector β: where δ is a Gaussian random noise with zero mean. The specific forms of M , β , δ , L (k) can be referred to the work 12 .
Usually, the least square estimate: is used to deal with this retrieval process 10,16,17 . Recently, considering the potential prior information, the work 12 further introduces the penalty mechanism into the retrieval process, and propose a penalized ENZ AR algorithm.
Global adaptive generative adjustment for estimating coefficients. As shown in the "Introduction" section, those Zernike coefficients follow from different prior distributions. Especially, their standard deviation could be dramatically diverse. So we assume that coefficient β j has a prior N(0, τ 2 j ) and the standard deviation τ j could be diverse. From the classical linear model, the observed response y is generated by a linear system Xβ + ǫ , where β = (β 1 , . . . , β p ) T can be viewed as the true signal, the Gaussian noise ǫ ∼ N(0, σ 2 I) and I is an identity matrix.
Considering the posterior distribution of coefficients, we further obtain an objective function with multiple hyperparameters {b j } . These hyperparameters can be viewed as the prior information of model parameters {β j } . In (14), n is the sample size, p is the dimension of the vector β , and the hyperparameter b j , j = 1, . . . , p, provides a customized shrinkage on the coefficient β j . www.nature.com/scientificreports/ We adopt a global adaptive generative adjustment (GAGA) algorithm to recover a true signal β . In Algorithm 1, hyperparameters and the signal are alternatively updated by a data-driven method. The inputs of this algorithm are the response vector y , the design matrix X , the iteration number K. The output of this algorithm is the signal estimate β = GAGA(y, X, K) . The convergency analysis of Algorithm 1 and the large sample properties of its output is discussed in the work 20 .
The following Algorithm 2 combines the ENZ AR process and the GAGA algorithm, which utilizes the customized prior information and updates model parameters and hyperparameters alternatively.

Results
We suggest a method adopting a global adaptive generative adjustment (GAGA) algorithm for estimating the ENZ coefficients. We call this method the GAGA ENZ AR. In a previos work 12 , the least absolute shrinkage and selection operator (Lasso) algorithm can also be applied to the aberration retrieval. So we compared the GAGA ENZ to the Lasso ENZ AR and the conventional ENZ AR 17 in the simulation with synthesized data. The characteristics of the optical system are shown in Table 1.
The simulations were implemented in three steps: (1) We simulated three PSFs (images intra, in, and extra focus) from (8) with the first Q leading terms {β m n } , where β 0 0 = 1 . In this paper, we set Q = 91 in the AR process. We further added Gaussian white noise to PSFs and simulated four noise levels (40 dB, 35 dB, 30 dB, 25 dB) measured using Signal-Noise Ratio (SNR).
(2) Using the simulated images with noise, we estimated {β m n } by ENZ AR, Lasso ENZ AR and GAGA ENZ AR separately.
(3) We divided the estimates {β m n } by β 0 0 to ensure that the first component of {β m n } is 1, and then assessed the experimental results on the residual square error �β − β� 2 2 , where β was the true parameter vector, and β was obtained by ENZ AR, Lasso ENZ AR and GAGA ENZ AR separately.
Different noise levels were chosen, and the above procedures were executed hundreds of time for each noise level. The empirical MSE of the true parameter β were calculated by (15): where i means the ith test. N is the total number of test. The standard deviation of MSE can be further computed by:  Fig. 4. These phases at constant amplitude further generated 100 generalized pupils, which can be expressed by (4) for Q = 91 in ENZ. Then we get 100 groups of {β m n } as the true parameters in experiments. We compared the empirical MSE (15) for the output of ENZ AR, Lasso ENZ AR and GAGA ENZ AR. In Fig. 5 we show the error bar of the empirical MSE at each given SNR, where the length of the bar is two times the standard deviation STD (16). Our method GAGA ENZ AR produced superior parameter estimates at each noise level compared to ENZ AR and Lasso ENZ AR. Though our algorithm works well for a small wavefront error (pv < 2π) , it may fail to handle the AR for strong wavefront aberrations in the atmospheric disturbances.   www.nature.com/scientificreports/ Real data example. The phase shown in Fig. 6a was observed from the measurement of non-common path aberrations from a 1.23 m adaptive optics telescope in Changchun China. More detailed descriptions on this optics telescope can be found in these works [21][22][23] . The field with constant amplitude and phase (Fig. 6a, rms = 0.14π ) has NZ coefficients shown in Fig. 6b,c.
We repeated the simulation one hundred times and showed the error bar of the MSE in Fig. 7.
When the noise SNR is 40, the conventional ENZ AR performed marginally than LASSO ENZ AR. However, GAGA ENZ AR outperformed ENZ AR and LASSO ENZ AR at each noise level. Moreover, GAGA ENZ AR also kept the empirical MSE in a low value even at a high noise level.

conclusion
We find that to represent an optical field, Zernike coefficients have their customized priors. For incorporating the customized prior information, we adopt a global adaptive generative adjustment method for estimating coefficients in the ENZ aberration retrieval. In the simulated and real data experiments, our algorithm provides better performance on the MSE compared to ENZ AR and LASSO ENZ AR.