Research on adaptive optics image restoration algorithm based on improved joint maximum a posteriori method

In this paper, we propose a point spread function (PSF) reconstruction method and joint maximum a posteriori (JMAP) estimation method for the adaptive optics image restoration. Using the JMAP method as the basic principle, we establish the joint log likelihood function of multi-frame adaptive optics (AO) images based on the image Gaussian noise models. To begin with, combining the observed conditions and AO system characteristics, a predicted PSF model for the wavefront phase effect is developed; then, we build up iterative solution formulas of the AO image based on our proposed algorithm, addressing the implementation process of multi-frame AO images joint deconvolution method. We conduct a series of experiments on simulated and real degraded AO images to evaluate our proposed algorithm. Compared with the Wiener iterative blind deconvolution (Wiener-IBD) algorithm and Richardson-Lucy IBD algorithm, our algorithm has better restoration effects including higher peak signal-to-noise ratio (PSNR) and Laplacian sum (LS) value than the others. The research results have a certain application values for actual AO image restoration.


Introduction
In the process of adaptive optics (AO) image, it usually takes post-processing for non-wavefront measurement data of the AO image, in order to get high-definition restoration images. Under the influence of atmospheric turbulence, it is difficult to determine the point spread function (PSF) of the AO image. According to the observed image, the joint estimation of both the PSF and object image should be carried out at the same time, and this process is called the blind image restoration technology [1].
There are many kinds of blind image restoration methods that can be applied to the AO image restoration, and researchers have been making efforts for blind deconvolution problems. Ayers et al. proposed a single frame blind deconvolution method [2], Tian et al. developed an AO image restoration based on the frame selection and multi-frame blind deconvolution method [3], Zhang et al. proposed a multi-frame iterative blind convolution algorithm based on improved expectation maximization for AO image restoration [4], and an adaptive image restoration method based on hierarchical neural Lijuan ZHANG et al.: Research on Adaptive Optics Image Restoration Algorithm Based on Improved Joint Maximum a Posteriori Method 23 network has been proposed by Yap et al. [5]. In this paper, to improve the resolution of the adaptive optics images, we propose to use joint maximum a posteriori (JMAP) estimation method for the blind deconvolution problem by finding the most suitable image data sets, at the same time estimated PSF and object image. This paper is organized as follows: in Section 2, we present the models on the observation, the PSF, and the image; Section 3 describes the JMAP estimation method for the AO image restoration; We present experimental results in Section 4, and conclusions are drawn in Section 5.

AO image model
Blind deconvolution refers to a class of problem for the form as follows: where ( , ) f x y is the original image, ( , ) h x y is the PSF, ( , ) g x y is the observed image, ( , ) n x y is the observation noise, Ω is the area of image, and the operator ⊗ denotes 2-dimensional (2-D) convolution.
During the actual observation process, we can obtain multi-frame AO images of the same object. Define a multi-frame AO images degradation model as m h x y denotes its corresponding PSF, and M is the total number of frames. Therefore, we use the sequence of the degraded multi-frame short exposure images { } 1 M m m g = of the same object to restore the original image f.

PSF reconstruction method
In this section, we propose the PSF reconstruction method that belongs to the model estimation approach, and the PSF estimated method will consider the degraded reasons including atmospheric turbulence, camera shank, and telescope of the adaptive optics system, ect [6].
In order to simplify the notation, we use one-dimensional pattern to describe the AO image. Theoretically, the mathematical model of the PSF after closed-loop correction of the AO system is shown as [6,7] w h e r e 1 within the lens aperture pupil function of the telescope, u is the coordinates for the pupil plane, λ is the central wavelength, l is the focal length of the AO system, and y is the coordinate of the pixel point. Due to the focus error or optical deviation in the AO system, the pupil function should be modified as θ is the wavefront phase error caused by the focus error or optical deviation. The PSF model can be given as follows: From (5), it is known that the PSF is dependent on ( ) u θ . Due to the effect of the atmospheric turbulence, the phase error will change over time, so the PSF model changing with time is given by where θ t (u) is the phase error which is caused by the atmospheric turbulence changing over time.

Proposed algorithm
There is a large number of works on blind deconvolution, originating in good part from astronomy [8]. The blind deconvolution method performs an estimation of both the PSF and the object jointly. We use the joint estimation based on a Bayesian framework for the JMAP estimator, that is  (  , ; ) ; where ( , , ; ) is the joint probability density of the original image g, the object f, and the PSF h.
The noise on the images is mainly photon noise which has a Poisson distribution [8]. However, a strong homogeneous background dominates the AO images. So under this condition, we will assume that the noise is white Gaussian distribution with a variance 2 σ for the short-exposure AO images. For the object, we select a stationary Gaussian distribution with a covariance matrix R f and a mean value f mean , which can be defined as N denotes the number of the pixels in the image, det( ) f R is the determinant of matrix x, and H is the operator performing the convolution by the PSF h.
Therefore, we can define the estimated object and coefficients forf and ĥ that the cost function is given as where the first item is the data-fidelity, and the second item is the regularization term. Then according to (8), we solve (9) by the criterion to be minimized; that is where C is a constant. This minimization is done using a conjugate gradient method. By minimizing (10), we obtain the proposed iterative algorithm as where ζ is a constant. According to (11) and (12), the object image f and the PSF ĥ can be obtained after a certain number of iterations. The following steps describe our proposed algorithm to restore an AO image blurred with an unknown spatially invariant PSF. The pseudo-code for our algorithm is summarized in Algorithm 1.

Algorithm 1 Steps for our proposed restoration algorithm
Step 1: Initialize operation. The M frames of images (g 1 , g 2 , …, g M ) are obtained with frame selection technique [3], then the initial object image is 0 Step 2: Obtain the initial estimation of PSF 0 h , according to the algorithm in Section 2.2.
a) The conjugate gradient method is used to optimize (11)

AO image restoration experiment
Because of the absence of the ideal image reference cases, we adopt the objective evaluation criteria peak signal-to-noise ratio (PSNR) [9] and Laplacian sum (LS) [10] in this paper. The PSNR is defined as  (13) where M and N represent the total numbers of pixels along the x-axis and y-axis of the image, respectively. (x, y) denotes a pixel location on the image, f(x, y) represents the gray value of point (x, y) on the ideal image, and ˆ( , ) f x y is the gray value of the image being estimated at point (x, y).
Laplacian sum is defined as follows: A number of experiments have been performed with our proposed algorithm using several synthetically degraded and real AO images, some of which are presented here.

Restoration experiment on simulated images
The test images are selected from [11], and the original image is the "star" of size 256 × 256 pixels; the blur PSF is the matrix [1, 6, 4, 6, 1] T [1, 6, 4, 6, 1]/256, and SNR = 25 dB, corresponding to a noise standard deviation of σ 2 = 8. Set parameters the same as AO telescope of 1.2 m on observatory in Yunnan, those are: the atmospheric coherence length r 0 = 13 cm, diameter of telescope D = 1.08 m, focal length l = 22.42 m, the wavelength in center λ = 700 nm, and the pixel size of charge coupled device (CCD) is 7.4 μm. Experiments on the iterative blind deconvolution (IBD) algorithm based on Richardson-Lucy iterative (RL-IBD) [12], iterative blind restoration algorithm based on the Wiener filtering (Wiener-IBD) [13] and our proposed algorithm are compared. Figure 1 is the original image, simulated degradation images, and the restored results. Figure  1(a) is the original image, and twenty frames degraded images with turbulence and noise are shown in Figs. 1(b) and 1(c) (to save space, only show two frames); Fig. 1(d) is the restored image

Photonic Sensors 26
using the Wiener-IBD algorithm; Fig. 1(e) is the restored image adopting the RL-IBD algorithm; Fig.  1(f) is the restored image using our algorithm, and the variance of wavefront phase error θ t (u) is 6.14 × 10 -3 , , that the parameter ζ = 0.25 was selected experimentally for visually acceptable results, extrinsic iterations for 100 times, inner loop of the PSF for 4, image cycles for 2, namely, and the total times of cycling is 600.
In order to verify our algorithm, we compare our algorithm with two other algorithms of Wiener-IBD and RL-IBD, and the objective evaluation criteria for the experimental results are measured by PSNR, LS, and the computation time, which are shown in Table 1. Compared with the Wiener-IBD and RL-IBD algorithms, we can see that the PSNR measures of our algorithm are increased by 5.42%, and 1.94%, respectively, and the LS values are increased by 15.6%, and 12.8%, respectively. It is shown in Table I that our algorithm can obtain the higher PSNR and LS. The computation load of our algorithm is slightly higher than those of the other two restoration algorithms, and we plan to further improve the performance of our approach in future.

Restoration experiment on AO image
The proposed image restoration method is tested on a set of AO data. A data set corresponds to the observations of the stellar on December 1, 2006 by 1.2 m AO telescope from the Chinese Academy of Sciences in Yunnan observatory. The size for imaging CCD is 320 × 240 pixels array, and the main parameters of the AO imaging system are as follows: the atmospheric coherent length r 0 is 13 cm, the telescope aperture D is 1.03 m, the imaging observation range is 0.7 -0.9μm, the central wavelength λ is 720 nm, the CCD pixel size is 6.7 μm, and the focal length l is 20 m. Figure 2 is the comparison of multi-frame stellar observation AO images and restoration results based on the three algorithms. Figures 2(a), 2(b), and 2(c) are the observation of stellar images (to save space show only 3); Fig. 2(d) is the restored object using the Wiener-IBD method; Fig. 2(e) is the restored object using the RL-IBD method; Fig. 2(f) is the restored object using our algorithm, and the variance of wavefront phase error θ t (u) = 6.04 × 10 -3 , cycles for 2 times, namely. The total iterations are 300; Fig. 2(g) is estimated PSF based on our algorithm. Figure 3 is the graph of comparison results for three restoration algorithms. Table 2 is the PSNR and the LS data, and the computation time results for the three restoration algorithms on the stellar AO images.  From visual and objective evaluations, we can see that our approach can effectively restore AO image, and the quality of the restored AO image has clearly improved.
The restored experiments on the simulated images and a set of real observed AO images are carried out, and the results show that our algorithm has better restoration effects including higher PSNR and higher LS than those of the other two methods. The calculation load of our algorithm is slightly higher than the others. The computational complexity of our proposed method has space per iteration, which are the same as our method. The main advantage of our method over the other two methods is that it incorporates AO devices parameters into the restoration process. This method largely prevents noise amplification from occurring. In the future work, we will optimize our proposed algorithm and consider the use of parallel processing to improve the restored effect of the AO image.

Conclusions
According to the characteristics of the AO image, we propose a novel adaptive optics images restoration algorithm based on the joint maximum a posteriori method. The PSF reconstruction method based on the wavefront phase effect is proposed. The restored experiments on the simulated images and a set of real observed AO images are carried out, and the results demonstrate that our proposed approaches result in high-quality restorations in both simulated and real AO image experiments.