Research on Adaptive Optics Image Restoration Algorithm by Improved Expectation Maximization Method

and Applied Analysis 3 where θ(⋅) is the wavefront phase error caused by the error of focal length or optical deviation. So the mathematical model of PSF can be written as h (y; θ) = 󵄨󵄨󵄨󵄨󵄨󵄨 ∫ p(u)e −j⋅(2π/λl)⋅y⋅u du 󵄨󵄨󵄨󵄨󵄨󵄨 2 = 󵄨󵄨󵄨󵄨󵄨󵄨 ∫ P(u)e jθ(u) e −j⋅(2π/λl)⋅y⋅u du 󵄨󵄨󵄨󵄨󵄨󵄨 2


Introduction
Adaptive optics (abbreviated as AO) is developed in order to overcome the interference of atmospheric turbulence on the telescope imaging, and it can measure and correct the optical wavefront phase of atmospheric turbulence in realtime.The technology matures; its application is extending to civilian fields besides the fields of large telescopes and laser engineering [1].In recent years, with the development of adaptive optics technology, many large ground-based optical telescopes are equipped with adaptive optics system [2] in the world.
When surface-to-air object imaging (e.g., astral target) is observed by ground-based optical telescopes, adaptive optics system can only realize part of the wavefront error correction, and the information of high frequency on target suppression and attenuation seriously [3,4].There are many factors such as atmospheric turbulence, the AO system error and photon noise on imaging path.Therefore, we acquire more clear images carrying on image restoration before adaptive optics correction.At present, there are many image restoration algorithms to overcome the influence of atmospheric turbulence, such as iterative blind deconvolution algorithm (abbreviated as IBD) [4,5], IBD based on Richardson-Lucy iteration (abbreviated as RL-IBD), IBD based on the Wiener algorithm (abbreviated as Wiener-IBD), and the algorithm based on the expectation maximization algorithm (abbreviated as EM) [6,7].
This paper proposes a blind image restoration method by multiframe iteration with improved adaptive optical images based on expectation maximization algorithm.Section 2 introduces the method of adaptive optics image restoration in detail including constructing adaptive optics model of imaging process, deducing model of point spread function on AO image, and researching the denoise processing method of AO image.Finally, it is given a deconvolution method using multiframe AO images by combining AO imaging system parameters with the regularization technique on the basis of expectation maximization algorithm.Section 3 verifies the validity of the algorithm in this paper by using simulate AO images, actual observation binary AO images, and stellar AO images.Section 4 concludes the paper.

Restoration Algorithm of Adaptive Optics Images
Adaptive optics system can measure, control, and correct the error of optical wavefront in real-time; Figure 1 shows the principle diagram of an adaptive optical system for imaging compensation.The system includes wavefront sensor (abbreviated as WFS), the deformable mirror (abbreviated as DM), high speed processor for wavefront, and other parts.The system can measure the distortion of optical wavefront in real-time which is caused by the influence of atmospheric turbulence and convert it to corresponding controlling signal adding to the wavefront correction element, which can compensate distortion of wavefront in real-time, and the optical telescope would be able to get close to target imaging for diffraction.AO images are corrected by AO system, but its compensation or correction is not sufficient; in order to get clearer target images, AO images restoration algorithm is researched.

The Imaging Process of Adaptive
Optics System.Considering the normal optical imaging model, first of all, the proper adaptive optics imaging model is set up, which lays the foundation of restoring high-quality adaptive optical images.The acquisition process of adaptive optics image model can consist of a degradation function and an additive noise [5]; the imaging process of AO in spatial frequency can be expressed as follows: where () is the original image for energy concentration, ℎ() is the point spread function (abbreviated as PSF) or impulse response for image system, () is output observation images for the AO system, () is the additive noise, and Ω is the image area.
In application, we use multiple frames of short exposure images {  }  =1 for one target to restore the target image ;  is the number of frames.The model of multiframe short exposure image with degradation can be written as where   () and ℎ  () represent the th frame of short exposure AO image with degradation for the target and its corresponding PSF, respectively.This paper adopts the selection techniques of frames referred to in [6] as the image preprocessing, which picks up the images with better quality from image sequence of short exposure images, in order to improve the solution of blind convolution and convergence of iteration.

Model of Point Spread Function on AO Images.
In adaptive optics system, the residual error of wavefront correction mainly consists of incomplete correction error by turbulence and noise of closed-loop system.
When AO system works in closed-loop, the transfer error function is the transfer function of wavefront phase perturbation caused by atmospheric turbulence, and the noise of system introduced by wavefront sensor is closed-loop transfer function [4].Theoretically, the mathematical model of the PSF after closed-loop correction of adaptive optics system is shown as follows: In the formula,  () = { 1, within the lens aperture, 0, others is the pupil function of telescope;  is the coordinates for the pupil plane;  is the central wavelength of imaging beam;  is the focal length of the optical system;  is the coordinates of two-dimensional focal plane.Due to the error of focal length and optical deviation in optical system, the pupil function should be rectified as follows: where The wavefront phase error is caused by the focal length error or optical deviation due to the influence of atmospheric turbulence; the PSF model changes with the focal length error or optical deviation over time as follows: In the formula,   (⋅) is the phase error caused by atmospheric turbulence changing over time.
In terms of AO system, the exposure time of observing images is shorter than turbulent fluctuation.If the intensity of the target is kept unchanged during exposure time, the model for series images is In the formula,   is the phase error of turbulence for the th frame short exposure AO image.

The Denoising Processing of AO Images.
In order to reduce the noise of AO images, Chinese researchers put forward the method to suppress noise based on total variation [8,9].This paper proposes a method to suppress noise that is based on the power spectral density (abbreviated as PSD) of image and the supporting domain field of constraints images.Reference [10] analyzed the method to minimize the noise in symmetrical support domain field of AO images.The analysis for this method is based on the variance ratio of real images.In this paper, the assessment about denoise result using noise changed function.The denoising process of AO image is actual a kind of image noise corrected process in Fourier spectral domain.If the AO image has the symmetrical support field, its noise mainly consist by CCD read-in noise and AO noise, and then noise power spectral density can be researched.Therefore, suppress noise for observation images (), and   () is constrained by field of support domain: In the formula, () is weighting function in support field, the image spectrum   () is: where  is the spatial frequency of image, () = () * (), () is the signal frequency spectrum function, () is the noise spectrum function, and bandwidth function () = () * .In the limited support domain field, () is very narrow, which nearly equals the  times of the width of ().
According to the definition of J. M. Conan, the PSD of an image is the Fourier transformation of its covariance, defining the power spectrum of noise and power spectrum of distortion image [11]: In the formula, () is the Fourier transformation of image (); () is the mean value of sample ().
Regard the circle as target supporting domain field; its diameter is from 8 pixels to 185 pixels.The definition of the changing image noise function   () is In the formula, PSD()   is the power spectral density of the image with constraints of support domain field and PSD()   is the power spectral density of the image without constraints of support domain field.

Improved Expectation Maximization
Algorithm.Reference [12] deduced the point spread function of multiframe degraded images with turbulence and the iteration formula of target image based on maximum likelihood criterion, but using this method leads to the larger solution space.However, we put forward the technique of regularization to modify the maximum likelihood estimation algorithm, introduce constraints in accordance with AO image in physical reality, and adopt the improved EM algorithm on multiframe iteration of AO images to acquire ideal f and { ĥ }.

EM Algorithm by Combining Regularization Technique.
EM algorithm is put forward by Dempster et al. [13] for parameters of maximum likelihood estimate of the iterative algorithm, which is iteratively going on two basic steps: calculation steps E (Expectation) and M (Maximization).In this paper, a new regularization method combined with prior knowledge is applied to the EM algorithm, realizing a method of multiframe AO image restoration.
Based on prior knowledge of target image and PSF joint probability distribution function, defining the maximization function will get target image estimation f and PSF estimation ĥ; that is, In the formula,  is the normalization factor.The cost function is So the problem is converted to minimize cost function (, ℎ); that is, [ f, ĥ] = min ,ℎ ((, ℎ)).The cost function of AO images and the optimization model of its parameter estimation are shown as follows: Symbol explanation: the first item represents the noise cost function, the second item represents the constraints cost function of PSF, and the third item is the cost functions with target edge constraint;  is prior information,  2  is the variance of mixed noise,   is the wavefront phase error, the power spectral density of the point spread function is PSD() ℎ = ⟨| Ĥ() − ()| 2 ⟩, Ĥ() is the Fourier transformation of ĥ(), () is the Fourier transform of average PSF, the estimation of PSD ℎ () can be estimated through the AO system closed-loop control data, and  and (∇) are regularization adjustment parameters, using the target edge constraint theory of Mugnier isotropic [14] as where  is a constant coefficient, (⋅) is a regularized function that has variable regularization coefficient associated with the gradient of each point, (∇) = exp −|∇| 2 /2 2 is based on [15], 0 < (∇) ≤ 1, parameter  is the controlled smoothing coefficient for selection, and ∇ is the gradient amplitude in four directions for the point (, ) (where ,  are pixel coordinates and  is gray value of pixel).This paper will apply all the prior information on the estimation of target and PSF, in order to get the optimal estimation of target and the point spread function (PSF).Combine formulas (2), (8), and (15) to establish the cost function for joint deconvolution with multiframe AO images: where  is the number of frames.

Implement of Image Restoration
Algorithm.This paper is combining with the basic idea of the EM algorithm [13]; through limited times of cycle calculations, find the most suitable PSF estimation; on the purpose of restoring target image drastically, modify the cost function (17) and make it correspond to full data [15]; if { g ( | )} is full data set of  times observation results   () of observed image,   () and the expectation of cost function are as follows: In the iteration process, it is supposed that the  − 1 iteration result is known;   ( multi (, ℎ)) represents the current expectation of the cost function.In order to minimize the expectation of the cost function, formula ( 18) can be differentiated by ℎ   () and   () then make the differentiation equal to zero; that is, The derivation process in detail is omitted, so the final solutions are f+1 ) .

(21)
In this paper, the specific steps of image restoration algorithm are as follows.
(1) According to the method of Section 2.3, constrain the support domain field of observed image and deal with the noise, and so forth, regard the result of one time Wiener filtering as the initial value of target estimation f0 .(2) According to the method of Section 2.2, for the initial PSF estimation ĥ0 , normalized processing should be taken at same time.
(3) Set the Max iteration of extrinsic cycles.
(4) Set Max count of target and PSF estimation.
(5) Set the inner loop counter variable h count by PSF estimation, variable f count of target inner cycle count.
Step 2 (the calculation of parameter values).
(1) According to the formulas of Section 2.3, the variance  2  of mixed noise of each degradation image is estimated.
(2) According to the method and calculation formulas from Section 2.2, calculate the power spectral density PSD ℎ () of PSF and estimate the wavefront phase error   .(3) According to the scheme and formulas of Section 2.4.1, estimate the regularization parameters  and (∇) of target edge.
(1) The inner loop count variable of PSF h count = 0.
(3) The inner loop counter variable of target estimation   = 0.
(5) Judge whether the extrinsic loop is over; if it is  >  , then turn to Step 4.
Step 4 (the algorithm is in the end).The cycle does not stop until the algorithm converges, and the output target image estimation is f that is over.

Results and Analysis of AO Image Restoration Experiments
3.1.The Restoration Experiment on Simulate AO Images.The algorithm proposed in this paper is used to do the simulation experiment on degradation spot images; the evaluation of image restoration quality is judged by using root mean square error (RMSE) and improved ΔSNR.
The evaluation of estimation accuracy is defined as follows: In the formula, 1 and 2 represent the total number of pixels on -axis and -axis on the image, respectively.(, ) represents a pixel point on image, and (, ) represents the gray value of the point (, ) on ideal image.
In the simulation experiments of this paper, the original image is infrared simulation imaging with spots for multiple targets with long distance; 10 frames of simulation degraded image series are generated by image degradation software from the key optical laboratory of Chinese Academy of Sciences, and then add Gaussian white noise and make the signal-to-noise ratio SNR of image for 20 dB.Set parameters the same as adaptive optical telescope of 1.2 m on observatory in Yunnan.The main parameters of telescope imaging system are the atmospheric coherence length  0 = 13 cm, diameter of telescope  = 1.08 m, comprehensive focal length of optical system  = 22.42 m, the wavelength in center  = 700 nm, and the pixel size of CCD is 7.4 m.Experiments on IBD algorithm based on RL iterative [16] (RL-IBD), iterative blind restoration algorithm based on the Wiener filtering [17] (Wiener-IBD), and the proposed algorithm are compared.

Restoration Experiments about Binary AO Images.
Without the reference of ideal image, using the objective quality evaluation with no reference image, the Full Width Half Maximum (abbreviated as FWHM) [18] is used as objective evaluation for image restoration in this paper; FWHM refers to imaging peak pixels is two times of half peak, including the FWHM on  direction and  direction; for the better evaluation in astronomical observation areas, computation formula is shown as follows: The experiment for the binary image restoration without reference images is carried out by the algorithm proposed in this paper.The binary images for experiments were taken by the 1.2 m AO telescope from Chinese Academy of Sciences in Yunnan observatory on December 3, 2006.The size for imaging CCD is 320 × 240 pixels array; the main parameters of imaging system are as follows: the atmospheric coherent length  0 = 13cm, telescope aperture  = 1.03, all imaging observation band is 0.7∼0.9,center wavelength  = 0.72, the CCD pixel size is 6.7 m, and optical imaging focal length  = 20.Recovery experiments based on RL IBD algorithm, Wiener-IBD algorithm, and the proposed algorithm in this paper are compared; the frame selection technology is taken in this paper; 40 frames are picked from 100 frames of degradation images as blind convolution images to be processed.
Figure 3 is the comparison results of multiframe degraded images and restoration images based on three kinds of algorithms.Figures 3(a), 3(b), 3(c), and 3(d) are for the observation of the binary image (to save space, we show only four frames; the rest is omitted); Figure 3(e) is the estimation of the point spread function (PSF) based on the algorithm in this paper; Figure 3(f) is the restored image based on the method proposed in this paper; the variance of wavefront phase error   (⋅) is 6.04 × 10 −3 , ĥ0 = 9.03 × 10 −3 , and  2  = 1.32 × 10 −5 so that the regularization parameter in this experiment is  = 1.15; the parameter of regularization function is  = 1.34, extrinsic iterations for 50 times, the inner loop of PSF for 4 times, and image loop for 2 times, namely.The total times of iteration are 300, FWHM = 6.17 pixels, and the restoration results are very close to the diffraction limit of 1.2 m AO telescope; Figure 3(g) is the restoration result based on RL-IBD algorithm, FWHM = 5.62 pixels; Figure 3(h) is the restoration result based on Wiener-IBD algorithm, FWHM = 4.76 pixels.Figure 4 is the sketch of energy for the observation image and three kinds of restored images based on three restoration algorithms.With the comparison, it shows that the proposed algorithm in this paper is better than Wiener-IBD algorithm and RL-IBD algorithm and needs less time of iterations.

Restoration Experiments about
Stellar AO Images.Also, the restoration experiment of stellar images was carried out using the algorithm in this paper; the stellar images for experiment were taken by 1.2 m AO telescope from the Chinese Academy of Sciences in Yunnan observatory on December 1, 2006 collection.The parameters of the optical imaging system are the same as that of binary images, and the method of parameters selection for experiments in this paper is nearly the same.Figure 5 is the comparison of multiframe stellar observation images and restoration results based on the three algorithms.

Conclusions
In this paper, according to the characteristic of adaptive optics image, we establish a degradation image model for multiframe exposure image series, and the mathematical model of the point spread function (PSF) for adaptive optical system after closed loop correction is derived; the PSF model based on phase error and changing over time is solved, and we use the power spectral density of image and constrain support domain method for AO image denoising processing.Finally, this paper applied the actual AO imaging system parameters as prior knowledge combined with regularization technique to improve the EM blind deconvolution algorithm.Experimental results show that the established blind deconvolution algorithm based on improved EM method can effectively eliminate the atmospheric turbulence and eliminate the noise pollution brought by the AO closed-loop correction, which improves the resolution of the adaptive optical image.The image process of AO observed images shows that the restoration quality is better than the same kind of IBD algorithm.It also shows that the restoration of binary images using the algorithm of this paper has carried on the identification of PSF and completed the AO image after processing; the results of actual AO image restoration have important application value.

Figure 1 :
Figure 1: Works of adaptive optical system for imaging compensation.

Figure 2 :Figure 3 :
Figure 2: Degraded images in experiment and restoration comparison of three algorithms.

Figure 4 :
Figure 4: Energy spectrum of observed image and restoration comparison of three algorithms.

4 Y
(a) Degraded image frame 1 (b) Degraded image frame 2 (c) Degraded image frame 3 (d) Degraded image frame 4 (e) Restored image by Weiner-IBD algorithm (f) Restored image by RL-IBD algorithm (g) Estimated PSF by our algorithm (h) Restored image by our algorithm Energy spectrum of observed image in Figure 5(a) ( p i x e l ) X (p ix el ) Gray (j) Energy spectrum of our algorithm restored image in Figure 5(h) (k) A line plot taken access the horizontal axis of image restored by Weiner-IBD, RL-IBD, and our algorithm

Figure 5 :
Figure 5: The multiframe observed stellar images and restoration comparison of three algorithms.

Figure 2
is the compared results of the original image, simulated multiframe degradation images, and restoration images based on three kinds of algorithm; the image size is 217 × 218, Figure 2(a) is the original target image, and ten frames turbulence degraded image sequences are generated by using image degradation simulation software, as shown in Figures 2(b), 2(c), and 2(d) (to save space, we show only three frames; the rest is omitted); Figure 2(e) is the Wiener-IBD restored image, RMSE = 16.35; Figure 2(f) is the RL-IBD restored image; the number of iterations is 800 times, RMSE = 10.31; Figure 2(g) is restored image used method in this paper, the variance of wavefront phase error   (⋅) is 6.14 × 10 −3 , ĥ0 = 8.26 × 10 −3 , and  2  = 1.1 × 10 −5 that the regularization parameters in the experiments are  = 1.25; the parameter of regularization function is  = 1.34, extrinsic iterations for 100 times, inner loop of the PSF for 4, image Figures 5(a), 5(b), 5(c), and 5(d) are the observation of stellar images (to save space, we show only 4); Figure 5(e) is for the Wiener-IBD restored image; Figure 5(f) is the RL-IBD restored image; Figure 5(g) is the estimation of the point spread function (PSF) based on algorithm in this paper, the iterative initial value ĥ0 = 9.56 × 10 −3 , PSF iterative for six times; Figure 5(h) is the restored image based on algorithm in this paper; Figure 5(i) is the energy of the observed image for stellar; Figure 5(j) is the energy of restoration image based on the algorithm in this paper; Figure 5(k) is the graph of comparison results for three restoration algorithms.

Table 1 :
Comparison of ΔSNR and  estimate values among three algorithms.