Deep Learning Estimation of Modified Zernike Coefficients and Recovery of Point Spread Functions in Turbulence

Recovering the turbulence-degraded point spread function from a single intensity image is important for a variety of imaging applications. Here, a deep learning model based on a convolutional neural network is applied to intensity images to predict a modified set of Zernike polynomial coefficients corresponding to wavefront aberrations in the pupil due to turbulence. The modified set assigns an absolute value to coefficients of even radial orders due to a sign ambiguity associated with this problem and is shown to be sufficient for specifying the intensity point spread function. Simulated image data of a point object and simple extended objects over a range of turbulence and detection noise levels are created for the learning model. The MSE results for the learning model show that the best prediction is found when observing a point object, but it is possible to recover a useful set of modified Zernike coefficients from an extended object image that is subject to detection noise and turbulence.


Introduction
The point spread function (PSF) refers to the impulse response of an imaging system [15].When the PSF is known, it can be used for the correction of blur and other artifacts in images that are due to the system's response.For example, in a space-invariant imaging situation, a PSF correction might be applied in a deconvolution step.Additionally, the propagation of light through a medium such as the atmosphere introduces wavefront aberrations at the aperture plane that further degrade the images.Therefore, the estimation of the combined system and medium PSF can potentially quantify the aberrations so that image correction can be performed.

Imaging Model
Incoherent imaging of an object can be modeled using a linear space-invariant forward model where I(x, y) is the intensity image of a source object I 0 (x, y), * is the convolution operator, and h(x, y) is the intensity PSF given by h(x, y) ∝ F p(x, y)e jφ(x,y) 2 , ( where F is the Fourier transform, | • | is the modulus operator, p(x, y) is the pupil function, and φ(x, y) is the wavefront phase distortion applied at the pupil plane [17].

Representing Wavefront Distortions
Although the PSF can be parameterized in various ways, to leverage results from disciplines such as adaptive optics, we consider the wavefront distortion φ(x, y) related to the PSF through Eq. ( 2).The phase distortion a nm Z m n (x, y), where Z m n are the Zernike polynomials and a nm are the Zernike coefficients [9,11,10,22].Zernike polynomials form a set of orthogonal functions defined over the unit circle and may be expressed in a double (m, n) indexing scheme [14] as where ρ is the normalized radial distance (0 ≤ ρ ≤ 1), ϕ is the azimuthal angle, n is the radial order, m is the angular frequency with n ≥ m ≥ 0, and R |m| n (ρ) are radial polynomials given by 0, otherwise. (5) For our work, it is more convenient to utilize the single-indexed Zernike polynomials Z q where the index q is related (OSA/ANSI standard) to the double-index by Zernike polynomials can be illustrated in a methodical way using a pyramid structure as given in Fig. 1.In the figure, the indexing starts with q = 0 on the top, then moves down the pyramid while scanning left to right.The rows of the pyramid correspond to the radial order n and the columns correspond to the angular frequency m.Using single-indexed Zernike polynomials, we rewrite Eq. (3) in polar coordinates as Due to the relationships between the intensity image I(x, y) in Eq. ( 1), the PSF h(x, y) in Eq. ( 2), and the Zernike coefficients a q in Eq. ( 7), the Zernike coefficients a q correspond to a particular intensity result and can thus be used to parameterize and model the associated PSF of an aberrated imaging system.
Unfortunately, there is an ambiguity associated with determining the Zernike coefficients from an intensity image alone without knowledge of the field phase.In particular, there are some Zernike polynomials Z q which generate the same PSF intensity image for oppositely-signed Zernike coefficients (details in Appendix).Due to this ambiguity, neural networks and other learning machines which have been extensively used to estimate PSFs for aberrated imaging systems [10, 20,18,28,24] face difficulty in directly predicting Zernike coefficients from intensity images [17,18].In [18], the deep learning model completely fails to predict Zernike coefficients for the in-focus setup from an extended source object due to this intensity image ambiguity (although the ambiguity is not pointed out or identified as the source of failure by the authors).In an attempt to improve performance, the authors in [18] analyze different preconditioners: overexposure, defocus, and scatter, and assess the impact on Zernike coefficient prediction.In [17], the authors recognize the PSF intensity image ambiguity and address it by performing multiple measurements (focused and defocused intensity images) of the source object, passing them through a feature extractor block followed by a neural network (Resnet50).
In our work, we propose a novel approach to address this ambiguity.We begin by identifying the Zernike polynomials which are susceptible to the ambiguity and for these polynomials we only consider the prediction of the absolute value of the associated Zernike coefficients.For the remainder of this work, we will refer to the signless Zernike coefficients for polynomials susceptible to the ambiguity in addition to the signed Zernike coefficients not susceptible to the ambiguity collectively as the modified Zernike coefficients.We use a deep neural network to predict the modified Zernike coefficients directly from intensity images.We show that these predicted coefficients can then be used to model the PSF of the aberrated imaging system.Finally, we consider the performance of our network in predicting the modified Zernike coefficients for both point source and extended source objects using incoherent imaging scenarios with both low and high noise levels (Poisson noise and read-out noise).Note that the proposed method is not a complete wavefront sensing approach, particularly because it does not consider the sign of the Zernike coefficients for polynomials that are susceptible to ambiguity.However, it works unambiguously for recovering the point spread function and provides a new construct and insight for the development of novel wavefront sensing methods [29,5,16,23,6].
The remainder of the work is organized as follows.In Section 2, we describe the methodology, data generation process, and deep learning architecture to predict modified Zernike coefficients.In Section 3, we illustrate our results and discuss the logical interpretation of these results.In Section 4, we conclude our work.

Methodology
As has been pointed out in Section 1.2, there is an ambiguity associated with predicting the Zernike coefficients from intensity images.In particular, we have found via symmetry properties of the Fourier transform [2] (details in Appendix) that angularly even Zernike polynomials which correspond to even values of radial order n and even values of angular frequency m generate the same PSF and thus the same intensity image for oppositely signed Zernike coefficients.In other words, the signs of these coefficients are immaterial for the description of the intensity PSF.Fig. 1 shows the first 28 polynomials where labels of angularly even polynomials are colored red.To address the ambiguity, we predict the signless Zernike coefficients for angularly even polynomials and the signed Zernike coefficients for angularly odd polynomials using a deep neural network.This approach might not be sufficient for applications such as wavefront sensing where the coefficient signs are important for specifying the actual shape of the wavefront, but it is useful for characterizing and correcting intensity images of an aberrated imaging system.

Data Generation
A dataset was generated using the imaging model and Zernike polynomial representation of the phase aberration given in Section 1.More specifically, the wavefront phase distortion due to the aberrated system was introduced in the pupil plane by fitting a Kolmogorov turbulence phase screen that is applied to the pupil.The strength of the phase screen was parameterized by the ratio D/r 0 , where D is the diameter of the pupil and r 0 is the Fried parameter (the coherence diameter of the transverse turbulence field) [25,27].Other simulation values include the wavelength λ = 0.5 µm and pupil diameter D = 0.4 m.A block diagram of our training methodology is given in Fig. 2. The process begins with the generation of a random phase screen φ(x, y) that obeys the Kolmogorov spectrum for a selected D/r 0 value [21].The first 28 Zernike coefficients of the phase screen φ(x, y) are fitted using matrix inversion [26].Other authors have demonstrated the effective estimation of this number of coefficients for turbulence-degraded wavefronts [17].Then, the coefficients are modified by discarding the signs of the coefficients corresponding to angularly even Zernike polynomials resulting in the modified Zernike coefficients {a q }.Additionally, the first three Zernike coefficients are considered as zero because there is no relationship between the Z 0 (piston) polynomial and the PSF intensity image [10], and the Z 1 (tip) and Z 2 (tilt) terms correspond to simple offsets in the image plane that can be easily found using centroiding algorithms or other registration methods [20,4].The modified phase screen φ(x, y) is reconstructed and used along with the pupil function p(x, y) and source object I 0 (x, y) in the imaging model described by Eq. ( 1) to generate the intensity image I(x, y) for neural network input.Finally, the modified Zernike coefficients {a q } are compared to the output of the neural network {â q } to determine the prediction error {Error q } associated with each Zernike coefficient.As an example, Fig. 3 shows a source object I 0 (x, y) and three intensity images I(x, y) corresponding to three different D/r 0 values.It is noted that it is possible to skip the phase screen generation step and use random draws of Zernike coefficient values subject to known variances corresponding to the Kolmogorov spectrum [19] but we find it useful to have a high spatial resolution version of the original phase screen for comparison purposes.
Four different scenarios at nine atmospheric turbulence strengths D/r 0 = 2, 3, . . ., 10 were considered.For each scenario, 54000 intensity images (512 × 512 pixels, 6000 images per D/r 0 value) were generated and then partitioned 45000/4500/4500 for training/validation/testing purposes.The first scenario considers only a point source object with zero noise.The remaining three scenarios consider extended source objects at zero/low/high noise levels, respectively.More specifically, 6000 different extended source objects (5000/500/500 for training/validation/testing purposes) were used for each D/r 0 .The extended source objects (28 × 28 pixels) were taken from the EMNIST database [3] and were zero-padded to 512 × 512 pixels.These objects have relatively simple spatial features, which provide the opportunity to extract the PSF even in the presence of wavefront components due to the object.For the low noise scenario, Poisson noise with peak photon levels of 4000 and read out noise with zero mean, standard deviation 10 was used.For the high noise scenario, Poisson noise with peak photon levels of 15000 and read out noise with zero mean, standard deviation 100 was used.These noise levels are similar to those used in [20].

Model Architecture
We use a deep neural network architecture based on AlexNet [7, 13] to predict the modified Zernike coefficients from PSF intensity images.The model structure is illustrated in Fig. 4. The model is implemented in the Julia language [1] using the Flux package [8].The 512 × 512 pixel PSF intensity images for the data set are centrally cropped down to 100 × 100 pixels and the amplitude is normalized using min-max normalization [7] prior to input to the neural network.The cropping procedure is carried out to exclude the outer zerovalue pixels that do not contain any significant information for further analysis.The model takes cropped (100 × 100 pixels) and amplitude normalized (in range [0, 1]) PSF images as input and outputs predictions for the modified Zernike coefficients that model the PSF of the aberrated imaging system.The network consists of 5 convolution layers, 3 max-pooling layers, and 3 fully-connected layers.Convolution layers 1 and 2 use 32 kernels with kernel size 5 × 5 and convolution layers 3, 4 and 5 use 64 kernels with kernel size 3 × 3.Each convolution layer in the model is followed by a RELU activation.Max pooling layer 1 performs a 5 × 5 kernel operation and max pooling layers 2 and 3 perform 3 × 3 kernel operations.The fully-connected layers 1 and 2 are followed by a RELU activation and have 2000 and 512 nodes, respectively.The final fully-connected layer has 25 nodes followed by a linear activation.We use mean squared error (MSE) as the loss function and the ADAM optimizer [12] with default settings (learning rate η = 0.001 and decay rates B = (0.9, 0.999)) as an optimizer for the model.The neural network is trained for 20 epochs with a batch size of 100 for each scenario described earlier as in Section 2.2.During the training process, the validation MSE typically levels off but does not diverge and the neural network with the lowest MSE for the validation data is saved for each scenario. Full

Results & Discussion
Table 1 shows prediction results for each of the four scenarios described in Section 2.2.Each of these results provides an average MSE for predicting the modified Zernike coefficients from the testing PSF intensity images.Additionally, when we repeated scenario 1, but instead of the modified Zernike coefficients we tried to predict the signed Zernike coefficients for all polynomials, we only achieved an MSE of 0.5673.The results from Table 1 show that scenario 1 with a point source object and no noise produces the lowest average MSE whereas scenario 4 with the extended source objects set and high noise gives the largest MSE.This is expected because of the simplicity of the object and the lack of noise in scenario 1.The introduction of extended source objects without noise in scenario 2 causes an increase in the MSE compared to scenario 1. Scenario 3 uses the same extended source objects as scenario 2 but with low noise.The average MSE increases about 20% compared to scenario 2. The presence of high noise in scenario 4 with the extended

Conclusion
A deep learning model based on a convolutional neural network was trained to predict modified Zernike coefficients in the pupil of an imaging system from a single turbulence-degraded intensity image.The modified Zernike coefficient set differs from a conventional set in that an absolute value is assigned to coefficients of even radial orders due to a sign ambiguity associated with using the intensity image.The modified set was shown to be sufficient for specifying the intensity PSF.Data for the learning model was created with an image simulation of a point object and simple extended objects for a range of turbulence and detection noise levels.The prediction MSE for the learning model shows that it is possible to recover a useful set of modified Zernike coefficients from an extended object intensity image subject to noise and turbulence.
As expected, the results show that the point source object with no noise produces the lowest average MSE whereas the extended source objects with high noise give the largest MSE.In all cases, the MSE increases in a predictable way with turbulence strength (D/r 0 ).Future work could explore the prediction of higher order Zernike terms and the use of more varied source objects.The quality and utility of the PSFs derived from the predicted Zernike coefficients were not investigated in this work and are essential topics for future efforts.
[10] Yuncheng Jin, Yiye   A Ambiguity Associated with Predicting Zernike coefficients from Intensity Images In order to explicitly show the ambiguity associated with predicting Zernike coefficients from intensity images, we utilize the symmetry properties of the Fourier transform [2].We begin by writing an expression for the point spread function h(x, y) in terms of the Fourier transform as h(x, y) ∝ F p(x, y)e jφ(x,y) F * p(x, y)e jφ(x,y) , and expressing the pupil function p(x, y) and wavefront phase distortion φ(x, y) in terms of real/imaginary and even/odd parts as p(x, y) = p re (x, y) + p ro (x, y) + jp ie (x, y) + jp io (x, y) and φ(x, y) = φ re (x, y) + φ ro (x, y) + jφ ie (x, y) + jφ io (x, y) where the subscripts re, ro, ie, and io represent the real-even, real-odd, imaginary-even, and imaginary-odd parts, respectively.Under the assumption that the pupil is a real-even function then only term p re (x, y) in ( 9) is non-zero.Furthermore, assuming the wavefront aberration phase is purely real, then only terms φ re (x, y) and φ ro (x, y) in (10) are non-zero.Next, we consider the two special cases of particular interest.Specifically, the cases when the Zernike polynomial is either angularly even which corresponds to even radial order, or angularly odd which corresponds to odd radial order.

A.1 Angularly Even Zernike Polynomial
An angularly even Zernike polynomial Z ae (x, y) must be a real even symmetric function [14].Thus, for an angularly even Zernike polynomial Z ae (x, y) the wavefront phase distortion φ re (x, y) has real even symmetry.Therefore, for an angularly even polynomial Z ae (x, y), we can write F p(x, y)e jφ(x,y) = F p re (x, y)e jφre(x,y) = F p re (x, y) cos φ re (x, y) + jp re (x, y) sin φ re (x, y) For Zernike coefficient ±a, we have φ re (x, y) = ±aZ ae (x, y) and F p(x, y)e ±jaZae(x,y) = F p re (x, y) cos[±aZ ae (x, y)] + jp re (x, y) sin[±aZ ae (x, y)] = F p re (x, y) cos[aZ ae (x, y)] ± jp re (x, y) sin[aZ ae (x, y)] = A(x, y) ± jB(x, y) where p re (x, y) cos[aZ ae (x, y)] is a real even symmetric function (cosine of a real even symmetric function is also real even symmetric) and jp re (x, y) sin[aZ ae (x, y)] is an imaginary even symmetric function (sine of a real even symmetric function is also real even symmetric).Using symmetry properties of the Fourier transform, A(x, y) = F p re (x, y) cos[aZ ae (x, y)] is real even and jB(x, y) = F jp re (x, y) sin[aZ ae (x, y)] is imaginary even.This leads to h(x, y) ∝ F p re (x, y)e jφre(x,y) 2 = A 2 (x, y) + B 2 (x, y).
Therefore, as can be seen from ( 16), an angularly even polynomial results in same value for h(x, y) irrespective of the sign of the Zernike coefficient a, thus leading to an intensity image ambiguity.

A.2 Angularly Odd Zernike Polynomial
Similarly, an angularly odd Zernike polynomial Z ao (x, y) must be a real odd symmetric function [14].Thus, for an angularly odd Zernike polynomial Z ao (x, y) the wavefront phase distortion φ ro (x, y) has real odd symmetry.Therefore, for an angularly odd polynomial Z ao (x, y), we can write F p(x, y)e jφ(x,y) = F p re (x, y)e jφro(x,y) ( = F p re (x, y) cos φ ro (x, y) + jp re (x, y) sin φ ro (x, y) For Zernike coefficient ±a, we have φ ro (x, y) = ±aZ ao (x, y) and F p(x, y)e ±jaZao(x,y) = F p re (x, y) cos[±aZ ao (x, y)] + jp re (x, y) sin[±aZ ao (x, y)] = F p re (x, y) cos[aZ ao (x, y)] ± jp re (x, y) sin[aZ ao (x, y)] = C(x, y) ± D(x, y) where p re (x, y) cos[aZ ao (x, y)] is a real even symmetric function (cosine of a real odd symmetric function is real even symmetric) and jp re (x, y) sin[aZ ao (x, y)] is an imaginary odd symmetric function (sine of a real odd symmetric function is also real odd symmetric).Using symmetry properties of the Fourier transform, C(x, y) = F p re (x, y) cos[aZ ao (x, y)] is real even and D(x, y) = F jp re (x, y) sin[aZ ao (x, y)] is real odd.This leads to h(x, y) ∝ F p re (x, y)e jφro(x,y) 2 = C(x, y) ± D(x, y) 2 . ( Therefore, as can be seen from ( 22), an angularly odd polynomial results in a different value for h(x, y) respective to the sign of the Zernike coefficient a, thus does not lead to an intensity image ambiguity.

A.3 Example
As an example, consider Figure 8 which shows the PSF intensity images for a point object modeled using Zernike polynomials Z 12 and Z 6 , respectively, with oppositely signed Zernike coefficients.Observe that for the angularly even polynomial Z 12 the intensity image is the same in 8(a) and 8(b) regardless of the sign of the Zernike coefficient.On the contrary, for the angularly odd polynomial Z 6 , the intensity image is different in 8(c) and 8(d) as a result of the sign of the Zernike coefficient and thus does not result in an ambiguity.

Figure 1 :
Figure 1: The first 28 Zernike polynomials arranged in a pyramid structure.The rows follow radial order n and columns follow angular frequency m [14].The angularly even polynomials which are subject to ambiguity in intensity images correspond to even values of n (corresponding labels are colored red).In the colormap, positive values are indicated as , zero values are indicated as , and negative value are indicated as .

Figure 2 :
Figure 2: Block diagram of the proposed training methodology.

Figure 4 :
Figure 4: The neural network architecture used for the estimation consists of several convolutional and pooling layers followed by three fully connected layers.The architecture is based on AlexNet [7, 13].

Fig. 5
shows the average MSE for predicting modified Zernike coefficients from the testing PSF intensity images as a function of D/r 0 ratios for the four scenarios.The MSE grows exponentially as a function of D/r 0 for all scenarios, and the ratios between the MSE values of the different scenarios remain relatively constant.Fig.6illustrates the estimation results for a point source object with zero noise.The diffraction limited point source object is shown in Fig.6(a) and the PSF intensity images I(x, y) for that point source with D/r 0 = 2, 5 are shown in Fig. 6(b) and Fig. 6(c) respectively.Fig. 6(d) and Fig. 6(e) display the actual and predicted Zernike coefficients corresponding to the intensity images in Fig. 6(b) and Fig. 6(c), respectively.As can be seen in Fig. 6(d) and Fig.6(e), the difference between the actual and predicted Zernike coefficients is small for both D/r 0 cases for point source object at zero noise.Similarly, Fig.7illustrates the estimation results for an extended source object with high noise.From Fig.7(d) and Fig.7(e), it may be observed that the deviation between actual and predicted Zernike coefficients increases significantly for both D/r 0 cases for the extended source with high noise, compared to Fig.6(d) and Fig. 6(e).

2 MSEFigure 5 :
Figure 5: The MSE as a function of D/r 0 for ( ) a point source object with no noise, ( ) extended source objects with no noise, ( ) extended source objects with low noise, and ( ) extended source objects with high noise

Figure 6 :
Figure 6: Demonstration of the proposed estimation for a point source object with zero noise.(a) Diffraction limited PSF for the point source, (b) I(x, y) for D/r 0 = 2, (c) I(x, y) for D/r 0 = 5, (d) the actual ( ) and predicted ( ) Zernike coefficients corresponding to the intensity image in (b), and (e) the actual ( ) and predicted ( ) Zernike coefficients corresponding to the intensity image in (c).The contrast was enhanced for display in (a)-(c).10

Figure 7 :
Figure 7: Demonstration of the proposed estimation for an extended source object with high noise.(a) I 0 (x, y) for the the extended source, (b) I(x, y) for D/r 0 = 2, (c) I(x, y) for D/r 0 = 5, (d) the actual ( ) and predicted ( ) Zernike coefficients corresponding to the intensity image in (b), and (e) the actual ( ) and predicted ( ) Zernike coefficients corresponding to the intensity image in (c).

Table 1 :
Average MSE for the four scenarios described in Section 2.2.