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

Authors: Lijuan Zhang; Dongming Li; Wei Su; Jinhua Yang; Yutong Jiang
Journal: Abstract and Applied Analysis
 (2014)

---

## 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.


---

## Body

## 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. Section2 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; Figure1 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.Figure 1
Works of adaptive optical system for imaging compensation.


### 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:
(1)g(x)=(f*h)(x)+n(x)=∫Ωf(y)h(x-y)dy+n(x),x,y∈Ω,
where f(x) is the original image for energy concentration, h(x) is the point spread function (abbreviated as PSF) or impulse response for image system, g(x) is output observation images for the AO system, n(x) is the additive noise, and Ω is the image area.In application, we use multiple frames of short exposure images{gm}m=1M for one target to restore the target image f; m is the number of frames. The model of multiframe short exposure image with degradation can be written as
(2)gm(x)=f(x)*hm(x)+n(x),1<m≤M,
where gm(x) and hm(x) represent the mth 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:
(3)h(y)=|∫P(u)e-j·(2π/λl)y·udu|2.
In the formula,
(4)P(u)={1,withinthelensaperture,0,others
is the pupil function of telescope; u is the coordinates for the pupil plane; λ is the central wavelength of imaging beam; l is the focal length of the optical system; y 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:(5)p(u)=P(u)ejθ(u),
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(6)h(y;θ)=|∫p(u)e-j·(2π/λl)·y·udu|2=|∫P(u)ejθ(u)e-j·(2π/λl)·y·udu|2.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:(7)h(y;θt)=|∫P(u)ejθt(u)e-j·(2π/λl)·y·udu|2.
In the formula, θt(·) 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(8)g(y;θtm,f)=∫h(y-x;θtm)f(x)dx.
In the formula, θtm is the phase error of turbulence for the mth 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 g(x), and gc(x) is constrained by field of support domain:
(9)gc(x)=g(x)w(x).
In the formula, w(x) is weighting function in support field, the image spectrum Gc(u) is:
(10)Gc(u)=∫G(v)W(u-v)du,
where u is the spatial frequency of image, G(u)=S(u)*N(u), S(u) is the signal frequency spectrum function, N(u) is the noise spectrum function, and bandwidth function W(u)=G(u)*δ. In the limited support domain field, W(u) is very narrow, which nearly equals the δ times of the width of G(u).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]:
(11)PSD(u)nn=〈|Re{G(u)}-Re{G-(u)}|2〉,PSD(u)ff=〈|Im⁡{G(u)}-Im⁡{G-(u)}|2〉.
In the formula, G(u) is the Fourier transformation of image g(x); G-(u) is the mean value of sample G(u).Regard the circle as target supporting domain field; its diameter is from 8 pixels to 185 pixels. The definition of the changing image noise functionFn(G) is
(12)Fn(G)=∫(PSD(u)nnc-PSD(u)nnc)du∫PSD(u)nncdu.
In the formula, PSD(u)nnc is the power spectral density of the image with constraints of support domain field and PSD(u)nnc 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 f^ and {h^k}.


#### 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 estimationf^ and PSF estimation h^; that is,
(13)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
In the formula, z is the normalization factor. The cost function is
(14)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
So the problem is converted to minimize cost function J(f,h); that is, [f^,h^]=min⁡f,h(J(f,h)). The cost function of AO images and the optimization model of its parameter estimation are shown as follows:
(15)J(f,h)=Jn(g∣f,h,a)+Jh(h∣a)+Jf(f∣a)s.t.Jn(g∣f,h,a)=∑i,j1σn2[g(i,j)-(h*f)(i,j)]2s.t.Jh(h∣a)=θt2∑u|H^(u)-H-(u)|2PSDh(u)s.t.Jf(f∣a)=α∑i,jΦ(r(f^,γ(i,j))).
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; a is prior information, σn2 is the variance of mixed noise, θt is the wavefront phase error, the power spectral density of the point spread function is PSD(u)h=〈|H^(u)-H-(u)|2〉, H^(u) is the Fourier transformation of h^(x), H-(u) is the Fourier transform of average PSF, the estimation of PSDh(u) can be estimated through the AO system closed-loop control data, and α and γ(∇x) are regularization adjustment parameters, using the target edge constraint theory of Mugnier isotropic [14] as
(16)Φ(p)=p-ln⁡(1+p),p(f^,γ(∇x))=∥∇f^(i,j)∥γ(∇x),∥∇f^(i,j)∥=((∇fx(i,j))2+(∇fy(i,j))2)1/2,
where α is a constant coefficient, γ(·) is a regularized function that has variable regularization coefficient associated with the gradient of each point, γ(∇x)=exp-|∇x|2/2ξ2 is based on [15], 0<γ(∇x)≤1, parameter ξ is the controlled smoothing coefficient for selection, and ∇x is the gradient amplitude in four directions for the point x(i,j) (where i,j are pixel coordinates and x 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:
(17)Jmulti(f,h)=Jn(gm∣f,hm,a)+Jh(hm∣a)+Jf(f∣a),
where m 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 {g~k(y∣x)} is full data set of k times observation results gk(y) of observed image, gk(y) and the expectation of cost function are as follows:
(18)E(gk(y))=E(∑xg~k(y∣x))=∑xhk(y-x)f(x)E(Jmulti(f,h))=E(Jn(∑xg~k∣f,hm,a))+E(Jh(hm∣a))+E(Jf(f∣a)).In the iteration process, it is supposed that then-1 iteration result is known; En(Jmulti(f,h)) represents the current expectation of the cost function. In order to minimize the expectation of the cost function, formula (18) can be differentiated by hkn(x) and fn(x) then make the differentiation equal to zero; that is,
(19)∂En(Jmulti(f,h))∂hkn(x)=0,∂En(Jmulti(f,h))∂fn(x)=0.
The derivation process in detail is omitted, so the final solutions are
(20)f^in+1(x)=f^in(x)+θtα*Φ(r(f^in(x),γ(x)))×∑y∈Y∑mh^in(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2),(21)h^in+1(x)=h^in(x)+θt2∑u|H^i(u)-H-i(u)|2PSDh_i(u)×∑y∈Yf^in+1(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2).In this paper, the specific steps of image restoration algorithm are as follows.Step 1 (initialization operation).
(1)
According to the method of Section2.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 estimationf^0.
(2)
According to the method of Section2.2, for the initial PSF estimation h^0, normalized processing should be taken at same time.
(3)
Set theMax_iteration of extrinsic cycles.
(4)
SetMax_count of target and PSF estimation.
(5)
Set the inner loop counter variableh_count by PSF estimation, variablef_count of target inner cycle count.Step 2 (the calculation of parameter values).
(1)
According to the formulas of Section2.3, the variance σn2 of mixed noise of each degradation image is estimated.
(2)
According to the method and calculation formulas from Section2.2, calculate the power spectral density PSDh(u) of PSF and estimate the wavefront phase error θt.
(3)
According to the scheme and formulas of Section2.4.1, estimate the regularization parameters α and γ(∇x) of target edge.Step 3 (iterate calculation,i=0,1,…,Max_iteration).
(1)
The inner loop count variable of PSFh_count= 0.
(2)
The iteration process of PSF estimationj=0,1,…,Max_count.
(a)
Complete the PSF estimationh^(i) with the conjugate gradient method using optimized formula (20).
(b)
Considerh_ count++;j++.
(c)
Judge the value of the loop variablej; if j<Max_count, continue; otherwise, turn to (3).
(3)
The inner loop counter variable of target estimationf_count=0.
(4)
The iteration process of target estimationk=0,1,…,Max_count.
(a)
The conjugate gradient method was used to optimize formula (21) to complete target image estimation f^(i);
(b)
considerf_count++; k++;
(c)
judge loop variablek; if k<Max_count, continue; otherwise, turn to (5).
(5)
Judge whether the extrinsic loop is over; if it isi>max_iteration, then turn to Step 4.
(6)
i
=
i
+
1, return to (1) 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 isf^ that is over.


## 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:
(1)g(x)=(f*h)(x)+n(x)=∫Ωf(y)h(x-y)dy+n(x),x,y∈Ω,
where f(x) is the original image for energy concentration, h(x) is the point spread function (abbreviated as PSF) or impulse response for image system, g(x) is output observation images for the AO system, n(x) is the additive noise, and Ω is the image area.In application, we use multiple frames of short exposure images{gm}m=1M for one target to restore the target image f; m is the number of frames. The model of multiframe short exposure image with degradation can be written as
(2)gm(x)=f(x)*hm(x)+n(x),1<m≤M,
where gm(x) and hm(x) represent the mth 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:
(3)h(y)=|∫P(u)e-j·(2π/λl)y·udu|2.
In the formula,
(4)P(u)={1,withinthelensaperture,0,others
is the pupil function of telescope; u is the coordinates for the pupil plane; λ is the central wavelength of imaging beam; l is the focal length of the optical system; y 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:(5)p(u)=P(u)ejθ(u),
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(6)h(y;θ)=|∫p(u)e-j·(2π/λl)·y·udu|2=|∫P(u)ejθ(u)e-j·(2π/λl)·y·udu|2.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:(7)h(y;θt)=|∫P(u)ejθt(u)e-j·(2π/λl)·y·udu|2.
In the formula, θt(·) 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(8)g(y;θtm,f)=∫h(y-x;θtm)f(x)dx.
In the formula, θtm is the phase error of turbulence for the mth 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 g(x), and gc(x) is constrained by field of support domain:
(9)gc(x)=g(x)w(x).
In the formula, w(x) is weighting function in support field, the image spectrum Gc(u) is:
(10)Gc(u)=∫G(v)W(u-v)du,
where u is the spatial frequency of image, G(u)=S(u)*N(u), S(u) is the signal frequency spectrum function, N(u) is the noise spectrum function, and bandwidth function W(u)=G(u)*δ. In the limited support domain field, W(u) is very narrow, which nearly equals the δ times of the width of G(u).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]:
(11)PSD(u)nn=〈|Re{G(u)}-Re{G-(u)}|2〉,PSD(u)ff=〈|Im⁡{G(u)}-Im⁡{G-(u)}|2〉.
In the formula, G(u) is the Fourier transformation of image g(x); G-(u) is the mean value of sample G(u).Regard the circle as target supporting domain field; its diameter is from 8 pixels to 185 pixels. The definition of the changing image noise functionFn(G) is
(12)Fn(G)=∫(PSD(u)nnc-PSD(u)nnc)du∫PSD(u)nncdu.
In the formula, PSD(u)nnc is the power spectral density of the image with constraints of support domain field and PSD(u)nnc 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 f^ and {h^k}.


### 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 estimationf^ and PSF estimation h^; that is,
(13)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
In the formula, z is the normalization factor. The cost function is
(14)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
So the problem is converted to minimize cost function J(f,h); that is, [f^,h^]=min⁡f,h(J(f,h)). The cost function of AO images and the optimization model of its parameter estimation are shown as follows:
(15)J(f,h)=Jn(g∣f,h,a)+Jh(h∣a)+Jf(f∣a)s.t.Jn(g∣f,h,a)=∑i,j1σn2[g(i,j)-(h*f)(i,j)]2s.t.Jh(h∣a)=θt2∑u|H^(u)-H-(u)|2PSDh(u)s.t.Jf(f∣a)=α∑i,jΦ(r(f^,γ(i,j))).
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; a is prior information, σn2 is the variance of mixed noise, θt is the wavefront phase error, the power spectral density of the point spread function is PSD(u)h=〈|H^(u)-H-(u)|2〉, H^(u) is the Fourier transformation of h^(x), H-(u) is the Fourier transform of average PSF, the estimation of PSDh(u) can be estimated through the AO system closed-loop control data, and α and γ(∇x) are regularization adjustment parameters, using the target edge constraint theory of Mugnier isotropic [14] as
(16)Φ(p)=p-ln⁡(1+p),p(f^,γ(∇x))=∥∇f^(i,j)∥γ(∇x),∥∇f^(i,j)∥=((∇fx(i,j))2+(∇fy(i,j))2)1/2,
where α is a constant coefficient, γ(·) is a regularized function that has variable regularization coefficient associated with the gradient of each point, γ(∇x)=exp-|∇x|2/2ξ2 is based on [15], 0<γ(∇x)≤1, parameter ξ is the controlled smoothing coefficient for selection, and ∇x is the gradient amplitude in four directions for the point x(i,j) (where i,j are pixel coordinates and x 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:
(17)Jmulti(f,h)=Jn(gm∣f,hm,a)+Jh(hm∣a)+Jf(f∣a),
where m 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 {g~k(y∣x)} is full data set of k times observation results gk(y) of observed image, gk(y) and the expectation of cost function are as follows:
(18)E(gk(y))=E(∑xg~k(y∣x))=∑xhk(y-x)f(x)E(Jmulti(f,h))=E(Jn(∑xg~k∣f,hm,a))+E(Jh(hm∣a))+E(Jf(f∣a)).In the iteration process, it is supposed that then-1 iteration result is known; En(Jmulti(f,h)) represents the current expectation of the cost function. In order to minimize the expectation of the cost function, formula (18) can be differentiated by hkn(x) and fn(x) then make the differentiation equal to zero; that is,
(19)∂En(Jmulti(f,h))∂hkn(x)=0,∂En(Jmulti(f,h))∂fn(x)=0.
The derivation process in detail is omitted, so the final solutions are
(20)f^in+1(x)=f^in(x)+θtα*Φ(r(f^in(x),γ(x)))×∑y∈Y∑mh^in(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2),(21)h^in+1(x)=h^in(x)+θt2∑u|H^i(u)-H-i(u)|2PSDh_i(u)×∑y∈Yf^in+1(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2).In this paper, the specific steps of image restoration algorithm are as follows.Step 1 (initialization operation).
(1)
According to the method of Section2.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 estimationf^0.
(2)
According to the method of Section2.2, for the initial PSF estimation h^0, normalized processing should be taken at same time.
(3)
Set theMax_iteration of extrinsic cycles.
(4)
SetMax_count of target and PSF estimation.
(5)
Set the inner loop counter variableh_count by PSF estimation, variablef_count of target inner cycle count.Step 2 (the calculation of parameter values).
(1)
According to the formulas of Section2.3, the variance σn2 of mixed noise of each degradation image is estimated.
(2)
According to the method and calculation formulas from Section2.2, calculate the power spectral density PSDh(u) of PSF and estimate the wavefront phase error θt.
(3)
According to the scheme and formulas of Section2.4.1, estimate the regularization parameters α and γ(∇x) of target edge.Step 3 (iterate calculation,i=0,1,…,Max_iteration).
(1)
The inner loop count variable of PSFh_count= 0.
(2)
The iteration process of PSF estimationj=0,1,…,Max_count.
(a)
Complete the PSF estimationh^(i) with the conjugate gradient method using optimized formula (20).
(b)
Considerh_ count++;j++.
(c)
Judge the value of the loop variablej; if j<Max_count, continue; otherwise, turn to (3).
(3)
The inner loop counter variable of target estimationf_count=0.
(4)
The iteration process of target estimationk=0,1,…,Max_count.
(a)
The conjugate gradient method was used to optimize formula (21) to complete target image estimation f^(i);
(b)
considerf_count++; k++;
(c)
judge loop variablek; if k<Max_count, continue; otherwise, turn to (5).
(5)
Judge whether the extrinsic loop is over; if it isi>max_iteration, then turn to Step 4.
(6)
i
=
i
+
1, return to (1) 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 isf^ that is over.


## 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 estimationf^ and PSF estimation h^; that is,
(13)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
In the formula, z is the normalization factor. The cost function is
(14)[f^,h^]=max⁡f,h(p(g∣f,h,a))=max⁡f,h(exp⁡(-J(f,h))z).
So the problem is converted to minimize cost function J(f,h); that is, [f^,h^]=min⁡f,h(J(f,h)). The cost function of AO images and the optimization model of its parameter estimation are shown as follows:
(15)J(f,h)=Jn(g∣f,h,a)+Jh(h∣a)+Jf(f∣a)s.t.Jn(g∣f,h,a)=∑i,j1σn2[g(i,j)-(h*f)(i,j)]2s.t.Jh(h∣a)=θt2∑u|H^(u)-H-(u)|2PSDh(u)s.t.Jf(f∣a)=α∑i,jΦ(r(f^,γ(i,j))).
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; a is prior information, σn2 is the variance of mixed noise, θt is the wavefront phase error, the power spectral density of the point spread function is PSD(u)h=〈|H^(u)-H-(u)|2〉, H^(u) is the Fourier transformation of h^(x), H-(u) is the Fourier transform of average PSF, the estimation of PSDh(u) can be estimated through the AO system closed-loop control data, and α and γ(∇x) are regularization adjustment parameters, using the target edge constraint theory of Mugnier isotropic [14] as
(16)Φ(p)=p-ln⁡(1+p),p(f^,γ(∇x))=∥∇f^(i,j)∥γ(∇x),∥∇f^(i,j)∥=((∇fx(i,j))2+(∇fy(i,j))2)1/2,
where α is a constant coefficient, γ(·) is a regularized function that has variable regularization coefficient associated with the gradient of each point, γ(∇x)=exp-|∇x|2/2ξ2 is based on [15], 0<γ(∇x)≤1, parameter ξ is the controlled smoothing coefficient for selection, and ∇x is the gradient amplitude in four directions for the point x(i,j) (where i,j are pixel coordinates and x 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:
(17)Jmulti(f,h)=Jn(gm∣f,hm,a)+Jh(hm∣a)+Jf(f∣a),
where m 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 {g~k(y∣x)} is full data set of k times observation results gk(y) of observed image, gk(y) and the expectation of cost function are as follows:
(18)E(gk(y))=E(∑xg~k(y∣x))=∑xhk(y-x)f(x)E(Jmulti(f,h))=E(Jn(∑xg~k∣f,hm,a))+E(Jh(hm∣a))+E(Jf(f∣a)).In the iteration process, it is supposed that then-1 iteration result is known; En(Jmulti(f,h)) represents the current expectation of the cost function. In order to minimize the expectation of the cost function, formula (18) can be differentiated by hkn(x) and fn(x) then make the differentiation equal to zero; that is,
(19)∂En(Jmulti(f,h))∂hkn(x)=0,∂En(Jmulti(f,h))∂fn(x)=0.
The derivation process in detail is omitted, so the final solutions are
(20)f^in+1(x)=f^in(x)+θtα*Φ(r(f^in(x),γ(x)))×∑y∈Y∑mh^in(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2),(21)h^in+1(x)=h^in(x)+θt2∑u|H^i(u)-H-i(u)|2PSDh_i(u)×∑y∈Yf^in+1(y-x)×(gm-∑z∈Xh^in(y-z)f^in(z)σi2).In this paper, the specific steps of image restoration algorithm are as follows.Step 1 (initialization operation).
(1)
According to the method of Section2.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 estimationf^0.
(2)
According to the method of Section2.2, for the initial PSF estimation h^0, normalized processing should be taken at same time.
(3)
Set theMax_iteration of extrinsic cycles.
(4)
SetMax_count of target and PSF estimation.
(5)
Set the inner loop counter variableh_count by PSF estimation, variablef_count of target inner cycle count.Step 2 (the calculation of parameter values).
(1)
According to the formulas of Section2.3, the variance σn2 of mixed noise of each degradation image is estimated.
(2)
According to the method and calculation formulas from Section2.2, calculate the power spectral density PSDh(u) of PSF and estimate the wavefront phase error θt.
(3)
According to the scheme and formulas of Section2.4.1, estimate the regularization parameters α and γ(∇x) of target edge.Step 3 (iterate calculation,i=0,1,…,Max_iteration).
(1)
The inner loop count variable of PSFh_count= 0.
(2)
The iteration process of PSF estimationj=0,1,…,Max_count.
(a)
Complete the PSF estimationh^(i) with the conjugate gradient method using optimized formula (20).
(b)
Considerh_ count++;j++.
(c)
Judge the value of the loop variablej; if j<Max_count, continue; otherwise, turn to (3).
(3)
The inner loop counter variable of target estimationf_count=0.
(4)
The iteration process of target estimationk=0,1,…,Max_count.
(a)
The conjugate gradient method was used to optimize formula (21) to complete target image estimation f^(i);
(b)
considerf_count++; k++;
(c)
judge loop variablek; if k<Max_count, continue; otherwise, turn to (5).
(5)
Judge whether the extrinsic loop is over; if it isi>max_iteration, then turn to Step 4.
(6)
i
=
i
+
1, return to (1) 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 isf^ 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ΔSNR. Consider
(22)RMSE=1N1N2∑i,j(f(i,j)-f^(i,j))2,ΔSNR=100log10∥g-f∥2∥f^-f∥2.The evaluation of estimation accuracy is defined as follows:(23)εestimate=∑i,j∈Sh|h(i,j)-h^(i,j)|2∑i,j∈Sh|h(i,j)|2.
In the formula, N1 and N2 represent the total number of pixels on x-axis and y-axis on the image, respectively. (i,j) represents a pixel point on image, and f(i,j) represents the gray value of the point (i,j) 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 lengthr0=13 cm, diameter of telescope D=1.08 m, comprehensive focal length of optical system l=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.Figure2 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 θt(·) is 6.14×10-3, h^0=8.26×10-3, and σn2=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 cycles for 2, namely, the total times of cycling is 600, RMSE = 10.05. Comparing with ΔSNR 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.Table 1
Comparison ofΔSNR and εestimate values among three algorithms.
Algorithm name
ΔSNR (dB)
ε
estimate
Iteration
Time-consuming
Wiener-IBD
1.982
3.62%
800
20.645 s
RL-IBD
3.046
4.49%
800
17.236 s
Our algorithm
4.164
5.10%
600
10.923 sDegraded images in experiment and restoration comparison of three algorithms.
(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 x direction and y direction; for the better evaluation in astronomical observation areas, computation formula is shown as follows:
(24)FWHM=FWHM_x2+FWHM_y2.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 is320×240 pixels array; the main parameters of imaging system are as follows: the atmospheric coherent length r0=13 cm, telescope aperture D=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 l=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.Figure3 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 θt(·) is 6.04×10-3, h^0=9.03×10-3, and σn2=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.Multiframe original images and restoration comparison of three algorithms.
(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 algorithmEnergy spectrum of observed image and restoration comparison of three algorithms.
(a)
Energy spectrum of observed image in Figure3(a)
(b)
Energy spectrum of our algorithm restored image in Figure3(f)
(c)
Energy spectrum of RL-IBD algorithm restored image in Figure3(g)
(d)
Energy spectrum of Wiener-IBD algorithm restored image in Figure3(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. Figure5 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 h^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.The multiframe observed stellar images and restoration comparison of three 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 Figure5(a)
(j)
Energy spectrum of our algorithm restored image in Figure5(h)
(k)
A line plot taken access the horizontal axis of image restored by Weiner-IBD, RL-IBD, and our algorithm


## 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. Consider
(22)RMSE=1N1N2∑i,j(f(i,j)-f^(i,j))2,ΔSNR=100log10∥g-f∥2∥f^-f∥2.The evaluation of estimation accuracy is defined as follows:(23)εestimate=∑i,j∈Sh|h(i,j)-h^(i,j)|2∑i,j∈Sh|h(i,j)|2.
In the formula, N1 and N2 represent the total number of pixels on x-axis and y-axis on the image, respectively. (i,j) represents a pixel point on image, and f(i,j) represents the gray value of the point (i,j) 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 lengthr0=13 cm, diameter of telescope D=1.08 m, comprehensive focal length of optical system l=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.Figure2 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 θt(·) is 6.14×10-3, h^0=8.26×10-3, and σn2=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 cycles for 2, namely, the total times of cycling is 600, RMSE = 10.05. Comparing with ΔSNR 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.Table 1
Comparison ofΔSNR and εestimate values among three algorithms.
Algorithm name
ΔSNR (dB)
ε
estimate
Iteration
Time-consuming
Wiener-IBD
1.982
3.62%
800
20.645 s
RL-IBD
3.046
4.49%
800
17.236 s
Our algorithm
4.164
5.10%
600
10.923 sDegraded images in experiment and restoration comparison of three algorithms.
(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 x direction and y direction; for the better evaluation in astronomical observation areas, computation formula is shown as follows:
(24)FWHM=FWHM_x2+FWHM_y2.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 is320×240 pixels array; the main parameters of imaging system are as follows: the atmospheric coherent length r0=13 cm, telescope aperture D=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 l=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.Figure3 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 θt(·) is 6.04×10-3, h^0=9.03×10-3, and σn2=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.Multiframe original images and restoration comparison of three algorithms.
(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 algorithmEnergy spectrum of observed image and restoration comparison of three algorithms.
(a)
Energy spectrum of observed image in Figure3(a)
(b)
Energy spectrum of our algorithm restored image in Figure3(f)
(c)
Energy spectrum of RL-IBD algorithm restored image in Figure3(g)
(d)
Energy spectrum of Wiener-IBD algorithm restored image in Figure3(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. Figure5 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 h^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.The multiframe observed stellar images and restoration comparison of three 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 Figure5(a)
(j)
Energy spectrum of our algorithm restored image in Figure5(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.



---

*Source: 781607-2014-07-06.xml*