Abstract
To improve the effect of adaptive optics images’ restoration, we put forward a deconvolution algorithm improved by the EM algorithm which joints multiframe adaptive optics images based on expectation-maximization theory. Firstly, we need to make a mathematical model for the degenerate multiframe adaptive optics images. The function model is deduced for the points that spread with time based on phase error. The AO images are denoised using the image power spectral density and support constraint. Secondly, the EM algorithm is improved by combining the AO imaging system parameters and regularization technique. A cost function for the joint-deconvolution multiframe AO images is given, and the optimization model for their parameter estimations is built. Lastly, the image-restoration experiments on both analog images and the real AO are performed to verify the recovery effect of our algorithm. The experimental results show that comparing with the Wiener-IBD or RL-IBD algorithm, our iterations decrease 14.3% and well improve the estimation accuracy. The model distinguishes the PSF of the AO images and recovers the observed target images clearly.
1. 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 real-time. 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.
2. 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.
2.1. 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 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.
2.2. 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, 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 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
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.
2.3. 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, is the power spectral density of the image with constraints of support domain field and is the power spectral density of the image without constraints of support domain field.
2.4. 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 and .
2.4.1. 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 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, . 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, is the variance of mixed noise, is the wavefront phase error, the power spectral density of the point spread function is , is the Fourier transformation of , is the Fourier transform of average PSF, the estimation of 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, is based on [15], , 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.
2.4.2. 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 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 iteration result is known; 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
In this paper, the specific steps of image restoration algorithm are as follows.
Step 1 (initialization operation). (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.(2)According to the method of Section 2.2, for the initial PSF estimation , 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 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 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.
Step 3 (iterate calculation, ). (1)The inner loop count variable of PSF h_count = 0.(2)The iteration process of PSF estimation .(a)Complete the PSF estimation with the conjugate gradient method using optimized formula (20).(b)Consider h_ count++; j++.(c)Judge the value of the loop variable ; if , continue; otherwise, turn to .(3)The inner loop counter variable of target estimation .(4)The iteration process of target estimation .(a)The conjugate gradient method was used to optimize formula (21) to complete target image estimation ;(b)consider f_count; ;(c)judge loop variable ; if , continue; otherwise, turn to .(5)Judge whether the extrinsic loop is over; if it is , then turn to Step 4.(6), return to and continue to loop.
Step 4 (the algorithm is in the end). The cycle does not stop until the algorithm converges, and the output target image estimation is that is over.
3. 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 . Consider
The evaluation of estimation accuracy is defined as follows: In the formula, and 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 cm, diameter of telescope m, comprehensive focal length of optical system m, the wavelength in center 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.
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 , 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, ; Figure 2(f) is the RL-IBD restored image; the number of iterations is 800 times, ; Figure 2(g) is restored image used method in this paper, the variance of wavefront phase error is , , and that the regularization parameters in the experiments are ; the parameter of regularization function is , extrinsic iterations for 100 times, inner loop of the PSF for 4, image cycles for 2, namely, the total times of cycling is 600, RMSE = 10.05. Comparing with and precision from Figures 2(e), 2(f), and 2(g) are shown in Table 1. Comparison results show that the proposed algorithm in this paper is better than Wiener-IBD algorithm and RL-IBD algorithm and needs less time of iterations.
(a) Original image
(b) Degraded image frame 1
(c) Degraded image frame 2
(d) Degraded image frame 3
(e) Restored image by Wiener-IBD
(f) Restored image by RL-IBD
(g) Restored image by our algorithm
3.2. 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 pixels array; the main parameters of imaging system are as follows: the atmospheric coherent length cm, telescope aperture , all imaging observation band is 0.7~0.9, center wavelength , the CCD pixel size is 6.7 m, and optical imaging focal length . 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 , , and so that the regularization parameter in this experiment is ; the parameter of regularization function is , 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, 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.
(a) Degraded image frame 1
(b) Degraded image frame 2
(c) Degraded image frame 3
(d) Degraded image frame 4
(e) Estimated PSF by our algorithm
(f) Restored image by our algorithm
(g) Restored image by RL-IBD algorithm
(h) Restored image by Weiner-IBD algorithm
(a) Energy spectrum of observed image in Figure 3(a)
(b) Energy spectrum of our algorithm restored image in Figure 3(f)
(c) Energy spectrum of RL-IBD algorithm restored image in Figure 3(g)
(d) Energy spectrum of Wiener-IBD algorithm restored image in Figure 3(h)
3.3. 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. 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 , 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.
(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
(i) Energy spectrum of observed image in Figure 5(a)
(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
4. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project is supported by the Scientific and Technological Research Project of Jinlin province, Department of Education, China (no. 2013145 and no. 201363), and the basis study project is supported by Jinlin Provincial Science & Technology Department, China (no. 201105055).