# Method of Image Quality Improvement for Atmospheric Turbulence Degradation Sequence Based on Graph Laplacian Filter and Nonrigid Registration

Authors: Jing Wu; Xiuqin Su
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

It is challenging to restore a clear image from an atmospheric degraded sequence. The main reason for the image degradation is geometric distortion and blurring caused by turbulence. In this paper, we present a method to eliminate geometric distortion and blur and to recover a single high-quality image from the degraded sequence images. First, we use optical flow technology to register the sequence images, thereby suppressing the geometric deformation of each frame. Next, sequence images are summed by a temporal filter to obtain a single blurred image. Then, the graph Laplacian matrix is used as the cost function to construct the regularization term. The final clear image and point spread function are obtained by iteratively solving the problem. Experiments show that the method can effectively eliminate the distortion and blur, restore the image details, and significantly improve the image quality.


---

## Body

## 1. Introduction

Atmosphere is an unstable random medium. Its temperature, pressure, density, and water vapor content are constantly changing and moving. Liquid water, dust, and aerosols are also continually changing and moving. Therefore, when light is propagated in the atmosphere, it is affected by molecular absorption, scattering of gas and aerosol, and turbulence disturbances. These effects degrade the quality of the captured video signal or image.In general, the atmosphere can be divided into a discrete turbid medium composed of mote and a “continuous” turbulent medium consisting of thermokinetic molecules. The turbidity medium is generally called aerosol particle. The influence on light propagation is mainly reflected in the scattering of particles, such as fog, haze, and dust clouds, which results in the deviation of light in any direction. Owing to the weak change of temperature, the turbulence media cause the weak fluctuation of refractivity to form turbulence. The influence on the propagation of light is mainly concentrated in the near-axis range of the line propagation, such as undulating hot air on asphalt pavement under the sun.The scattering and absorption effect of the turbid media on light radiation leads to energy attenuation. The quality of the image is degraded and the contrast is reduced, which is equivalent to the deterioration of the optical transmission function (OTF) of the optical system. A typical practical application of light propagation in turbid media is the visibility problem.The turbulent atmospheric medium causes the fluctuation of the refractive index of the atmosphere, which leads to the distortion of the light wavefront. The original wavefront of the light wave is thus changed, the light coherence is destroyed, and the image of the focal plane of the imaging device is seriously blurred. The atmospheric boundary layer, which is particularly near the ground, is the region of atmospheric turbulence and has the most complex characteristics. Turbulence is one of the most significant problems in physics. After several centuries of research on this topic, there remains no real solution. Adaptive optics (AO) [1–6] and image restoration technology [7–11] are used to eliminate atmospheric turbulence and obtain high-quality image. Given that AO is expensive and its processes are complex and AO image may still be blurred, it is necessary to study image restoration.At present, two main types of restoration process can be described for video or sequential images. The first process employs image selection and fusion, which involves the so-called “lucky imaging technique” [12–17]. Those regions of the input frames that have the best quality in the temporal direction are employed. The region of interest is selected and image fusion technology, such as image averaging, pixel-based processing, or region-based processing, is applied. Finally, a deblurring process may be reapplied. The disadvantage of lucky imaging technique is that it can only be applied to short exposure images. Due to the low probability of showing good images, a large number of images are usually required for the quality evaluation. However, some image sequences contain both large area distortion (perceived as waving) and small disturbance (perceived as jittering). This will negatively impact the quality of the final image.The other method employs an image registration technique with deformation estimation [18, 20, 21]. It combines image registration with blind deconvolution to restore images. In such methods, the image degradation process can be modeled using the following equation: y=D(H(x))+n, where y and x denote the degraded image and the latent image, respectively, and n is the noise. We assume that D is the geometric distortion factor and H is the blurring factor. Image registration process attempts to resolve small camera movements and temporal variations due to atmospheric refraction. Subsequently, a deblurring process is applied to the registered image and image details are restored. Gal’s method [20] applies state-of-the-art registration and deblurring algorithms as building blocks within the sum and deblur framework. Accurate preregistration of the input video frames narrows the spatial support of the effective blur kernel. Hence, powerful registration is therefore crucial for successful reconstruction.Motivated by these developments, in this paper, we revisit the turbulence recovery problem and incorporate state-of-the-art registration and deblurring algorithms within the turbulence recovery framework. For registration, we employ optical flow method to estimate the displacement between the two frames and eliminate the geometric distortion caused by atmospheric turbulence. Subsequently, the single blurred image is obtained by the temporal filter of the sequence image after image registration. In this step, the graph Laplacian regularization term is introduced to construct the cost function, suppress noise, and maintain image detail. Finally, a latent sharp image can be obtained by a blind deconvolution and deblurring process. The contributions of this paper can be listed as follows:(1)
An optical flow method is introduced to achieve the image registration between frames, and the pixel displacement between frames is estimated to eliminate the geometric distortion.(2)
Applying the concept of graphs to the image, the normalized graph Laplacian regularization term is constructed, and its filtering characteristics are used to reduce the noise of the image.The remainder of this paper is organized as follows. The proposed scheme for atmospheric degradation is detailed in Section2. The proposed method for image restoration is introduced in Section 3. The performance of the method is evaluated on a set of images and compared with other techniques in Section 4. Finally, Section 5 presents the discussion and conclusion of the paper.


## 2. Atmospheric Degradation Model

It is considered that turbulence is invariable for a certain time. When the working frequency of the detector is greater than the turbulence, we believe that the degraded image is subjected to the same turbulence effect during this period. Therefore, it can be described as the degradation process of the image caused by the spatial and temporal changing. The degraded model of spatial and temporal variation is established as the space-invariant blur kernel and deformation matrix. The space-invariant blur kernel can be regarded as a convolution of the clear image and space-invariant point spread function (PSF). The deformation matrix can calculate the pixel displacement of each frame image by nonrigid image registration. In the near-surface-captured video or image, the degradation process is described through the following linear model:(1)yt=DtHx+ntThe above equation shows that the captured image is blurred due to the influence of optical system and turbulence and degraded by time varying distortions. Here,x is the original static scene, yt is captured image at time t, H is a space and time invariant blur kernel, Dt is the space and time dependent deformation matrix at time t, and nt is white Gaussian noise at time t. In this model, the above variables are a lexicographically ordered vector representation. It is well known that restoring a scene distorted by atmospheric turbulence is a challenging problem in video surveillance. Assuming that the noise of the degraded image is not considered (because the noise can be considerably reduced by applying a simple temporal filter), the inverse process of the degradation can be described as(2)x=Dt-1H-1ytThe deformation matrixDt describes the displacement of each pixel between image yt and image x. Image registration is used to reduce the turbulence-induced distortions. After being processed by the temporal filter, the noise in registered images reduces and images are less blurred, leading to better reconstruction. Therefore, we can estimate Dt by image registration and then perform an inverse mapping on yt so that the registered sequence images are only affected by the blur kernel H. Because image deconvolution is an ill-posed problem, a prior constraint is introduced to transform the solution problem into a posed problem. In addition, the blur kernel H is unknown, its sparse feature is regarded as a regularization term, and the latent sharp image x is obtained by solving the reconstructed cost function.


## 3. Proposed Restoration Algorithm

The proposed restoration algorithm includes two parts: nonrigid image registration and blind deblurring of a single image. In the first place, the nonrigid image registration technique is applied to process the degraded sequenceyn. Then, a single blurred image is obtained through a temporal filter. The second step introduces a graph Laplacian regularization term to construct the cost function, and it solves the minimum solution of the cost function as the final restored image. For senses without ground truth, the reference image is obtained by a temporal filter. For senses with ground truth, the reference image is the ground truth image. The design of the temporal filter is described in detail in [22, 23]. It is also proved by a statistical argument that the median filter is a better choice than the mean filter and provides better results. The specific algorithm and block diagram are presented in Algorithm 1 and Figure 1, respectively.Algorithm 1:Overall algorithm.
Require: Observe all frames of blurry image yn, reference image f, and maximum kernel sizeh, maximum iteration number k. 1. According to the steps of Algorithm2 for image registration, obtain distortion factor Dt of each image frame and the registered sequences rn. 2. Determine the medianrm of each frame in the registered sequence images r n by the temporal filter. 3. Conduct blind estimation of blur matrixH from blurred image rm. Alternate: for i=0 to k-1 do (i) Update latent image x using Algorithm 3. (ii) Update blurring matrix H using the IRLS algorithm. end for Return the latent sharp imagexAlgorithm 2:Nonrigid image registration algorithm.
Require:Observe all frames of degraded image yn and reference image f 1. Minimize the cost function to obtain the displacement mapu. Alternate: (i) Solve equation (11) and update u. (ii) Solve equation (12) and update v. 2. According to the displacement mapu, calculate registered sequences rn. Return registered image sequencern.Algorithm 3:Update sharp image x.
Require: Blur kernel H from the previous H update, Image x from previous x update Require: Regularization parameter η. (i) ComputeK from x using equation (17) (ii) Apply Sinkhorn algorithm in [19] to K to obtain diagonal matrix C. (iii) Compute filtering matrixW and graph Laplace matrix L. (iv) Solve cost function using CG to computex. return update imagexFigure 1
Block diagram of the proposed image restoration method.


### 3.1. Optical Flow Image Registration

Image registration methods include rigid registration and nonrigid registration. Rigid registration is only aimed at shift and rotation transformation of image coordinates, which can not eliminate complex deformation. In this case, nonrigid image registration method is a better choice. Currently, several nonrigid image registration algorithms, such as a warping technique based on diffeomorphisms [24], B-spline image registration [25], optical flow technology [26], and other methods, are widely applied in the image registration. Image registration based on diffeomorphism is being increasingly applied to medical image processing. The limitation of image registration based on B-pline is that nonrigid registration is only effective when the geometric distortion can be adjusted by all control points in the grid. If nonrigid registration has to be used to compensate small disturbance, the number of control points will be huge, making the computation impossible. Because the result of image registration is the initial input to the deblurring method in this paper, fast and accurate image registration technology is needed. Meanwhile, optical flow technology has a high accuracy in calculating the movement displacement of objects. Therefore, we choose the image registration method based on optical flow technology.Optical flow was first proposed by Gibson as early as 1950; however, the truly effective optical flow algorithm was proposed much later by Horn and Schunck in 1981. The latter researchers creatively combined the two-dimensional velocity field with a grayscale. In addition, they introduced the optical flow constraint equation to structure the basic algorithm of the optical flow function [26]. Based on this idea, a large number of improved algorithms were proposed. Nagel adopted a conditional smoothing constraint, which means that the gradient is smoothed by a weighted matrix control [27]. Anandan proposed a piecewise smoothing method for multimotion estimation problem [28]. Fennema and Thompson’s method [29], which is based on the optical flow constraint equation, was introduced in the image discontinuity region under the assumption that velocity is constant. Their method cannot handle rotating targets because it is based on cluster analysis.In this paper, the TV-L1 optical flow algorithm [30] proposed by Zach is used. The algorithm combines the variation regularization with the data fidelity of the L1 norm. This design can maintain the discontinuity of the flow field. Moreover, it has strong robustness to illumination changes, shading, and noise enhancement.Two images frames,I0 and I1:(Ω⊆R2)→R, are given. The objective is to find the disparity map u:Ω→R2, which minimizes an image-based error criterion together with a regularization force. The TV-L1 algorithm focuses on the average intensity difference between pixels as the image similarity score. Hence, we can minimize the following formula to obtain disparity map u:(3)∫ΩλϕI0x-I1x+ux+ψu,∇u,…dxwhere ϕ(I0(x)-I1(x+u(x))) is the data fidelity and ψ(u,∇u,…) is the regularization term of the shape prior. Additionally, λ is the parameter between the fidelity and regularization term. We define ϕ(x)=x and ψ(∇u)=∇u; hence, the energy equation can be expressed as(4)E=∫ΩλI0x-I1x+ux+∇udxSince the equation is not continuously differentiable, a mathematical framework proposed by Chambolle [31] is used to solve the above problem. In one-dimensional space, assuming that the horizontal direction is a nonzero displacement, we can use the scalar u(x) to represent the disparity map u(x) and x+u(x) to x+ux,0T. We linearize the image I1 near x+u0:(5)I1x+u=I1x+u0+u-u0I1xx+u0where I1x is the derivative of the image intensity I1 in the x direction. Fixing u0 and employing the linear approximation for I1, the TV-L1 equation can be written as(6)E=∫ΩλuI1x+I1x+u0-u0I1x-I0+∇udxWe define the current residualI1x+u0+(u-u0)I1x-I0 as ρ(u,u0,x), which is abbreviated as ρ(u). Moreover, we introduce an auxiliary variable v and propose minimizing the convex approximation of the above equation:(7)Eθ=∫Ω∇u+12θu-v2+λρvdxwhere θ is a small constant and v is an approximation of u. This convex minimization problem can be solved by the alternating iterative updating of u and v.Fixingv, we solve(8)minu⁡∫Ω∇u+12θu-v2dxFixingu, we solve(9)minv⁡∫Ω12θu-v2+λρvdxWhen extended to high dimensions, the first-order image residualsρ(u,u0,x) are for two N-dimensional images, I0 and I1. The disparity map u0 is given as I1(x+u0)+∇I1,u-u0-I0(x). Additionally, we write ud for the dth component of ud∈1,2,…N.The equation generalization to more dimensions is the following energy:(10)Eθ=∫Ω∑d∇ud+∑d12θud-vd2+λρvdxSimilar to the one-dimensional case, the method of alternating updating is used to solve the problem.For everyd and fixed vd, we solve(11)minud⁡∫Ω∇ud+12θud-vd2dxForu being fixed, we solve(12)minv⁡∑d12θud-vd2+λρvThe disparity mapu can be obtained by solving (11) and (12), and the specific solution is given in [30]. Finally, the displacement of the reference image can be obtained and the geometric distortion can be eliminated.Through Algorithm2, we obtain the registered sequences rn and remove the geometric distortion factor Dt. The specific experimental results are shown in Figure 2, and the specific information of the image is given in Table 1.Table 1
Degraded sequences.
Type Number Name Color Type Resolution Total Frames Real datasetswith ground truth [10] A1A2 ChimneyBuilding greygrey 237 × 237 237×237 100100  Real datasets without ground- truth One- degreeTurbulence B1B2 B3 BlockhouseChimney Poles greygrey grey 280 × 280 260×260 280×280 100100 100  Real datasets without ground- truth Two- degree Turbulence C1C2 C3 Blockhouse ChimneyPoles greygrey grey 280 × 280 260×260 280×280 100100 100  Real datasets without ground- truth Three- degreeTurbulence D1D2 D3 BlockhouseChimney Poles greygrey grey 280 × 280 260×260 280×280 100100 100Figure 2
Registered result of optical flow for A1-A2 sequences. (a) Degraded sequence A1. (b) Registered sequence A1. (c) Degraded sequence A2. (d) Registered sequence A2.
(a) (b) (c) (d)


### 3.2. Single Image Deblur

After a nonrigid image registration step, the degraded image eliminates the geometric distortion and only the blur factor exists. Then, a single blurred image is obtained by the temporal filter.The blurry imagey we observe comes from a sharp image x blurred by blur kernel h along with the addition of Gaussian noise n, which can be expressed as(13)y=h⊗x+nOur goal is to recover the unknown imagex and the blurring kernel h. If x, y, and n are represented as lexicographically ordered vectors x, y, and n, the model can be rewritten as(14)y=Hx+nwhere H is a blurring matrix which is constructed from the corresponding PSF and usually has a special structure depending on the type of boundary condition assumptions.Most existing deblurring methods rely on optimizing a cost function as(15)Ex,H=y-Hx2+ηRx+λRHHere, we introduce the normalized graph Laplacian filter proposed in [32] to construct the regularization term. Meanwhile, we introduce the sparsity of blur kernel H to reduce the noise on the blur kernel. The constraint (nonnegative, sum of one) of H obeys the physical mechanism of blur formation. Finally, we structure the cost function as follows:(16)Ex,H=FLy-Hx2+ηRx+λRHwhere y and x represent the blurred image and latent image, respectively. H represents the blur factor, which is the reconstruction of PSF. The first term of the cost function is the data fidelity, where F(L) is a function of the graph Laplacian matrix, and F(L)=(I+βL). The second term R(x) is the constraint of the latent image x, which is related to the graph Laplacian. The third term R(H) is the constraint of the blur kernel H. To reduce noise in the blurring kernel, we use the L1 regularization as its constraint and the expression is RH=H. This constraint on H (sum-to-1 and nonnegativity) follows from physical principles of blur formation. The parameters η and λ affect the correlation between the image regularization and the blur kernel. The form of R(x) will be derived next.


#### 3.2.1. Graph Laplacian Regularization

In this section, we first outline the definition of a graph and then introduce the definition and construction method of the graph Laplacian. As depicted in Figure3, any image can be defined as an intensity function on the vertices V of a weighted graph G=(V,E,K) consisting of a finite set V of vertices (image pixels) and a finite set E⊂V×Vof edges (i,j) with the corresponding weights K(i,j). These weights measure the similarity between vertices (pixels) i and j in the graph. For the image of a n×n size, the value of its intensity function can be defined by x=x1,x2,⋯,xNT, where N=n2, and the weight of similarity K can be expressed with a symmetric nonnegative matrix of the N×N size. The graph Laplacian matrix can be derived from K.Figure 3
Graph representation of an image and construct of kernel similarity matrixK. (a) Graph representation of the image. (b) Construct of kernel similarity matrix K.
(a) (b)There are four different definitions of the graph Laplacian commonly used in the literature in the context of graph signal and image processing, each of which has different spectral characteristics [32–35]. An unnormalized Laplacian L=D-K is proposed in [33], where D is a diagonal matrix whose ith diagonal element is the sum of all elements of the ith row of matrix K. In [34], the author introduces the traditional normalized graph Laplacian L=I-D-1/2KD-1/2. In [35], another form of the graph Laplacian L=I-D-1K is given. In [32], a new normalized graph Laplacian L=I-C-1/2KC-1/2 is proposed and its characteristics are analyzed. Because of its excellent characteristics in mathematics, we use it as a constraint on the image in the cost function of our method. Then, we introduce an approach to obtaining this graph Laplacian matrix and outline its excellent features.The similarity matrix involved in the graph Laplacian is derived from the formula given in [36]. As shown in Figure 2(b), the pixel (i,j) of the similarity matrix in the image is defined as a nonlocal means (NLM)(17)Ki,j=exp⁡-x^i-x^j2h2in which x^i and x^j are patches around the pixels i and j of image x^ and h is a smoothing parameter. Owing to the similarity matrix K being a symmetric nonnegative matrix, we only calculate the similarity between each patch and its small neighborhood. Therefore, matrix K is sparse. The sparse structure is appealing from a computational point of view. Then, applying the Sinkhorn algorithm [19] to matrix K will return a diagonal scaling matrix C, such that the resulting matrix W=C-1/2KC-1/2 is a symmetric nonnegative doubly stochastic filtering matrix. Finally, the normalized graph Laplacian matrix L=I-W=I-C-1/2KC-1/2 is given. The specific process is shown in Figure 4.Figure 4
Diagram of the Laplacian matrix calculation process.It is pointed out in [32] that W is symmetric and can be decomposed as W=VSVT, where V is an orthonormal matrix whose columns are the eigenvectors of W and S=diag⁡λ1,λ2,⋯,λN is a diagonal matrix consisting of eigenvalues of W as its eigenvalues in the range of 1=λ1>λ2≥⋯≥λN≥0 [37]. It can be observed that the maximum eigenvalue is just equal to one with the corresponding DC eigenvector v1=(1/N)1,1,⋯,1T. This means that if matrix W is applied to the signal, the DC component of the signal can be protected. Spectral analysis of the matrix W shows that it is intrinsically low-pass (the largest eigenvalue corresponds to the DC component) [38], which is a very good property for the filter. The normalized Laplacian matrix L=I-W has high-pass filter characteristics (DC eigenvectors correspond to their null eigenvalues). This property is consistent with the expected behavior of Laplacian filters in image processing.The definition of the difference operator ofx:V→R of the edge (i,j)∈E in the weighted graph G is given in [32](18)dxi,j=Ki,jxjCj,j-xiCi,iwhere C(i,i) and C(j,j) are the corresponding j and ith diagonal elements of the diagonal matrix C derived from the Sinkhorn algorithm [19].The weighted gradient vector ofx at vector i∈V can be expressed as(19)∇xi=dxi,j1,⋯,dxi,jmT,∀i,j∈ETherefore, the Laplacian ofx at vector i can be written as(20)Δxi=∇2xi=1Ci,i∑j,j~iKi,jxiCi,i-xjCj,jAs a result, the regularization termR(x) can be expressed as(21)Rx=12∑i=1N∇xi2=12∑i=1N∑j,j~iKi,jxiCi,i-xjCj,j2=xTLx


#### 3.2.2. Solve the Cost Function

In this step, we give the process of solving the cost equation. Based on the above deduction, we can finally rewrite the cost function expression as(22)Ex,H=I+βLy-Hx2+ηxTLx+λHFurther, the cost function is rewritten in the following form:(23)Ex,H=I+βL1/2y-Hx2+ηL1/2x2+λHwhere I+βL=I+β(I-W)=VΛVT is a symmetric and positive semidefinite matrix. Therefore, the matrix I+βL1/2=VΛ1/2VT has a filtering behavior similar to that of I+βL. By using the eigendecomposition matrix W, we can compute matrix Λ, whose ith diagonal element is 1+β(1-λi). Here, λi is the ith diagonal element of diagonal matrix S decomposed by matrix W. Since Laplacian matrix L is a high-pass filter, with β>0, I+βL1/2 corresponds to a sharpening filter for the residual y-Hx. The same analysis applies to the second term in (23), where both Laplacian matrix L and its square root L1/2 are adaptive high-pass filters. Therefore, the regularization term in (23) adaptively punishes the high frequency of the final solution. Moreover, it suppresses the noise in the final solution while maintaining the image details.For solving (23), we use the method of alternating estimates of x and H to solve for the minimum of E(x,H). Specific means of updating x and H are outlined below.(A) Update x. In this step, we fix H and obtain x as the minimization subproblem:(24)min⁡I+βL1/2y-Hx2+ηL1/2x2It is apparent that thex subproblem is a convex function. To minimize (24), the corresponding gradient is equal to zero as(25)∂E∂x=-2HTI+βLy-Hx+2ηLx=0This can be written as a symmetric positive definite system of liner equations(26)HTI+βLH+ηLx=HTI+βLyA fast computational method, such as a conjugate gradient (CG) optimization algorithm, can be used to solve the above system. Note thatH and HT are a blur kernel PSF or its reconstructed forms, respectively. The matrix I+βL must have a square root. The parameter β should maintain the above system as positive, thereby satisfying β≥-1. The value of regularization parameters η and β are selected based on the noise variance and blurring scenario. They are fixed at each step of the algorithm. In this paper, the parameter β lies in the range of (0,1) and the parameter η is empirically selected in the range of (0,0.4). The algorithm for updating the latent image x is in Algorithm 3.(B) Update H. In this step, we fix x and obtain H of the minimization subproblem:(27)min⁡I+βL1/2y-Hx2+λHBlur kernelH satisfies the condition H≥0, ∑iHi=1. In our experiments, the parameter λ depends on the user-specified kernel size h, according to the formula: λ=(3/13)h. We use the unconstrained iterated reweighted least squares (IRLS) [39] and then project the results to the constraint conditions (negative elements are set to zero and are renormalized). During the iteration, the weights can be calculated from the last updated kernel. The IRLS internal iteration is solved by CG algorithm.


## 3.1. Optical Flow Image Registration

Image registration methods include rigid registration and nonrigid registration. Rigid registration is only aimed at shift and rotation transformation of image coordinates, which can not eliminate complex deformation. In this case, nonrigid image registration method is a better choice. Currently, several nonrigid image registration algorithms, such as a warping technique based on diffeomorphisms [24], B-spline image registration [25], optical flow technology [26], and other methods, are widely applied in the image registration. Image registration based on diffeomorphism is being increasingly applied to medical image processing. The limitation of image registration based on B-pline is that nonrigid registration is only effective when the geometric distortion can be adjusted by all control points in the grid. If nonrigid registration has to be used to compensate small disturbance, the number of control points will be huge, making the computation impossible. Because the result of image registration is the initial input to the deblurring method in this paper, fast and accurate image registration technology is needed. Meanwhile, optical flow technology has a high accuracy in calculating the movement displacement of objects. Therefore, we choose the image registration method based on optical flow technology.Optical flow was first proposed by Gibson as early as 1950; however, the truly effective optical flow algorithm was proposed much later by Horn and Schunck in 1981. The latter researchers creatively combined the two-dimensional velocity field with a grayscale. In addition, they introduced the optical flow constraint equation to structure the basic algorithm of the optical flow function [26]. Based on this idea, a large number of improved algorithms were proposed. Nagel adopted a conditional smoothing constraint, which means that the gradient is smoothed by a weighted matrix control [27]. Anandan proposed a piecewise smoothing method for multimotion estimation problem [28]. Fennema and Thompson’s method [29], which is based on the optical flow constraint equation, was introduced in the image discontinuity region under the assumption that velocity is constant. Their method cannot handle rotating targets because it is based on cluster analysis.In this paper, the TV-L1 optical flow algorithm [30] proposed by Zach is used. The algorithm combines the variation regularization with the data fidelity of the L1 norm. This design can maintain the discontinuity of the flow field. Moreover, it has strong robustness to illumination changes, shading, and noise enhancement.Two images frames,I0 and I1:(Ω⊆R2)→R, are given. The objective is to find the disparity map u:Ω→R2, which minimizes an image-based error criterion together with a regularization force. The TV-L1 algorithm focuses on the average intensity difference between pixels as the image similarity score. Hence, we can minimize the following formula to obtain disparity map u:(3)∫ΩλϕI0x-I1x+ux+ψu,∇u,…dxwhere ϕ(I0(x)-I1(x+u(x))) is the data fidelity and ψ(u,∇u,…) is the regularization term of the shape prior. Additionally, λ is the parameter between the fidelity and regularization term. We define ϕ(x)=x and ψ(∇u)=∇u; hence, the energy equation can be expressed as(4)E=∫ΩλI0x-I1x+ux+∇udxSince the equation is not continuously differentiable, a mathematical framework proposed by Chambolle [31] is used to solve the above problem. In one-dimensional space, assuming that the horizontal direction is a nonzero displacement, we can use the scalar u(x) to represent the disparity map u(x) and x+u(x) to x+ux,0T. We linearize the image I1 near x+u0:(5)I1x+u=I1x+u0+u-u0I1xx+u0where I1x is the derivative of the image intensity I1 in the x direction. Fixing u0 and employing the linear approximation for I1, the TV-L1 equation can be written as(6)E=∫ΩλuI1x+I1x+u0-u0I1x-I0+∇udxWe define the current residualI1x+u0+(u-u0)I1x-I0 as ρ(u,u0,x), which is abbreviated as ρ(u). Moreover, we introduce an auxiliary variable v and propose minimizing the convex approximation of the above equation:(7)Eθ=∫Ω∇u+12θu-v2+λρvdxwhere θ is a small constant and v is an approximation of u. This convex minimization problem can be solved by the alternating iterative updating of u and v.Fixingv, we solve(8)minu⁡∫Ω∇u+12θu-v2dxFixingu, we solve(9)minv⁡∫Ω12θu-v2+λρvdxWhen extended to high dimensions, the first-order image residualsρ(u,u0,x) are for two N-dimensional images, I0 and I1. The disparity map u0 is given as I1(x+u0)+∇I1,u-u0-I0(x). Additionally, we write ud for the dth component of ud∈1,2,…N.The equation generalization to more dimensions is the following energy:(10)Eθ=∫Ω∑d∇ud+∑d12θud-vd2+λρvdxSimilar to the one-dimensional case, the method of alternating updating is used to solve the problem.For everyd and fixed vd, we solve(11)minud⁡∫Ω∇ud+12θud-vd2dxForu being fixed, we solve(12)minv⁡∑d12θud-vd2+λρvThe disparity mapu can be obtained by solving (11) and (12), and the specific solution is given in [30]. Finally, the displacement of the reference image can be obtained and the geometric distortion can be eliminated.Through Algorithm2, we obtain the registered sequences rn and remove the geometric distortion factor Dt. The specific experimental results are shown in Figure 2, and the specific information of the image is given in Table 1.Table 1
Degraded sequences.
Type Number Name Color Type Resolution Total Frames Real datasetswith ground truth [10] A1A2 ChimneyBuilding greygrey 237 × 237 237×237 100100  Real datasets without ground- truth One- degreeTurbulence B1B2 B3 BlockhouseChimney Poles greygrey grey 280 × 280 260×260 280×280 100100 100  Real datasets without ground- truth Two- degree Turbulence C1C2 C3 Blockhouse ChimneyPoles greygrey grey 280 × 280 260×260 280×280 100100 100  Real datasets without ground- truth Three- degreeTurbulence D1D2 D3 BlockhouseChimney Poles greygrey grey 280 × 280 260×260 280×280 100100 100Figure 2
Registered result of optical flow for A1-A2 sequences. (a) Degraded sequence A1. (b) Registered sequence A1. (c) Degraded sequence A2. (d) Registered sequence A2.
(a) (b) (c) (d)


## 3.2. Single Image Deblur

After a nonrigid image registration step, the degraded image eliminates the geometric distortion and only the blur factor exists. Then, a single blurred image is obtained by the temporal filter.The blurry imagey we observe comes from a sharp image x blurred by blur kernel h along with the addition of Gaussian noise n, which can be expressed as(13)y=h⊗x+nOur goal is to recover the unknown imagex and the blurring kernel h. If x, y, and n are represented as lexicographically ordered vectors x, y, and n, the model can be rewritten as(14)y=Hx+nwhere H is a blurring matrix which is constructed from the corresponding PSF and usually has a special structure depending on the type of boundary condition assumptions.Most existing deblurring methods rely on optimizing a cost function as(15)Ex,H=y-Hx2+ηRx+λRHHere, we introduce the normalized graph Laplacian filter proposed in [32] to construct the regularization term. Meanwhile, we introduce the sparsity of blur kernel H to reduce the noise on the blur kernel. The constraint (nonnegative, sum of one) of H obeys the physical mechanism of blur formation. Finally, we structure the cost function as follows:(16)Ex,H=FLy-Hx2+ηRx+λRHwhere y and x represent the blurred image and latent image, respectively. H represents the blur factor, which is the reconstruction of PSF. The first term of the cost function is the data fidelity, where F(L) is a function of the graph Laplacian matrix, and F(L)=(I+βL). The second term R(x) is the constraint of the latent image x, which is related to the graph Laplacian. The third term R(H) is the constraint of the blur kernel H. To reduce noise in the blurring kernel, we use the L1 regularization as its constraint and the expression is RH=H. This constraint on H (sum-to-1 and nonnegativity) follows from physical principles of blur formation. The parameters η and λ affect the correlation between the image regularization and the blur kernel. The form of R(x) will be derived next.


### 3.2.1. Graph Laplacian Regularization

In this section, we first outline the definition of a graph and then introduce the definition and construction method of the graph Laplacian. As depicted in Figure3, any image can be defined as an intensity function on the vertices V of a weighted graph G=(V,E,K) consisting of a finite set V of vertices (image pixels) and a finite set E⊂V×Vof edges (i,j) with the corresponding weights K(i,j). These weights measure the similarity between vertices (pixels) i and j in the graph. For the image of a n×n size, the value of its intensity function can be defined by x=x1,x2,⋯,xNT, where N=n2, and the weight of similarity K can be expressed with a symmetric nonnegative matrix of the N×N size. The graph Laplacian matrix can be derived from K.Figure 3
Graph representation of an image and construct of kernel similarity matrixK. (a) Graph representation of the image. (b) Construct of kernel similarity matrix K.
(a) (b)There are four different definitions of the graph Laplacian commonly used in the literature in the context of graph signal and image processing, each of which has different spectral characteristics [32–35]. An unnormalized Laplacian L=D-K is proposed in [33], where D is a diagonal matrix whose ith diagonal element is the sum of all elements of the ith row of matrix K. In [34], the author introduces the traditional normalized graph Laplacian L=I-D-1/2KD-1/2. In [35], another form of the graph Laplacian L=I-D-1K is given. In [32], a new normalized graph Laplacian L=I-C-1/2KC-1/2 is proposed and its characteristics are analyzed. Because of its excellent characteristics in mathematics, we use it as a constraint on the image in the cost function of our method. Then, we introduce an approach to obtaining this graph Laplacian matrix and outline its excellent features.The similarity matrix involved in the graph Laplacian is derived from the formula given in [36]. As shown in Figure 2(b), the pixel (i,j) of the similarity matrix in the image is defined as a nonlocal means (NLM)(17)Ki,j=exp⁡-x^i-x^j2h2in which x^i and x^j are patches around the pixels i and j of image x^ and h is a smoothing parameter. Owing to the similarity matrix K being a symmetric nonnegative matrix, we only calculate the similarity between each patch and its small neighborhood. Therefore, matrix K is sparse. The sparse structure is appealing from a computational point of view. Then, applying the Sinkhorn algorithm [19] to matrix K will return a diagonal scaling matrix C, such that the resulting matrix W=C-1/2KC-1/2 is a symmetric nonnegative doubly stochastic filtering matrix. Finally, the normalized graph Laplacian matrix L=I-W=I-C-1/2KC-1/2 is given. The specific process is shown in Figure 4.Figure 4
Diagram of the Laplacian matrix calculation process.It is pointed out in [32] that W is symmetric and can be decomposed as W=VSVT, where V is an orthonormal matrix whose columns are the eigenvectors of W and S=diag⁡λ1,λ2,⋯,λN is a diagonal matrix consisting of eigenvalues of W as its eigenvalues in the range of 1=λ1>λ2≥⋯≥λN≥0 [37]. It can be observed that the maximum eigenvalue is just equal to one with the corresponding DC eigenvector v1=(1/N)1,1,⋯,1T. This means that if matrix W is applied to the signal, the DC component of the signal can be protected. Spectral analysis of the matrix W shows that it is intrinsically low-pass (the largest eigenvalue corresponds to the DC component) [38], which is a very good property for the filter. The normalized Laplacian matrix L=I-W has high-pass filter characteristics (DC eigenvectors correspond to their null eigenvalues). This property is consistent with the expected behavior of Laplacian filters in image processing.The definition of the difference operator ofx:V→R of the edge (i,j)∈E in the weighted graph G is given in [32](18)dxi,j=Ki,jxjCj,j-xiCi,iwhere C(i,i) and C(j,j) are the corresponding j and ith diagonal elements of the diagonal matrix C derived from the Sinkhorn algorithm [19].The weighted gradient vector ofx at vector i∈V can be expressed as(19)∇xi=dxi,j1,⋯,dxi,jmT,∀i,j∈ETherefore, the Laplacian ofx at vector i can be written as(20)Δxi=∇2xi=1Ci,i∑j,j~iKi,jxiCi,i-xjCj,jAs a result, the regularization termR(x) can be expressed as(21)Rx=12∑i=1N∇xi2=12∑i=1N∑j,j~iKi,jxiCi,i-xjCj,j2=xTLx


### 3.2.2. Solve the Cost Function

In this step, we give the process of solving the cost equation. Based on the above deduction, we can finally rewrite the cost function expression as(22)Ex,H=I+βLy-Hx2+ηxTLx+λHFurther, the cost function is rewritten in the following form:(23)Ex,H=I+βL1/2y-Hx2+ηL1/2x2+λHwhere I+βL=I+β(I-W)=VΛVT is a symmetric and positive semidefinite matrix. Therefore, the matrix I+βL1/2=VΛ1/2VT has a filtering behavior similar to that of I+βL. By using the eigendecomposition matrix W, we can compute matrix Λ, whose ith diagonal element is 1+β(1-λi). Here, λi is the ith diagonal element of diagonal matrix S decomposed by matrix W. Since Laplacian matrix L is a high-pass filter, with β>0, I+βL1/2 corresponds to a sharpening filter for the residual y-Hx. The same analysis applies to the second term in (23), where both Laplacian matrix L and its square root L1/2 are adaptive high-pass filters. Therefore, the regularization term in (23) adaptively punishes the high frequency of the final solution. Moreover, it suppresses the noise in the final solution while maintaining the image details.For solving (23), we use the method of alternating estimates of x and H to solve for the minimum of E(x,H). Specific means of updating x and H are outlined below.(A) Update x. In this step, we fix H and obtain x as the minimization subproblem:(24)min⁡I+βL1/2y-Hx2+ηL1/2x2It is apparent that thex subproblem is a convex function. To minimize (24), the corresponding gradient is equal to zero as(25)∂E∂x=-2HTI+βLy-Hx+2ηLx=0This can be written as a symmetric positive definite system of liner equations(26)HTI+βLH+ηLx=HTI+βLyA fast computational method, such as a conjugate gradient (CG) optimization algorithm, can be used to solve the above system. Note thatH and HT are a blur kernel PSF or its reconstructed forms, respectively. The matrix I+βL must have a square root. The parameter β should maintain the above system as positive, thereby satisfying β≥-1. The value of regularization parameters η and β are selected based on the noise variance and blurring scenario. They are fixed at each step of the algorithm. In this paper, the parameter β lies in the range of (0,1) and the parameter η is empirically selected in the range of (0,0.4). The algorithm for updating the latent image x is in Algorithm 3.(B) Update H. In this step, we fix x and obtain H of the minimization subproblem:(27)min⁡I+βL1/2y-Hx2+λHBlur kernelH satisfies the condition H≥0, ∑iHi=1. In our experiments, the parameter λ depends on the user-specified kernel size h, according to the formula: λ=(3/13)h. We use the unconstrained iterated reweighted least squares (IRLS) [39] and then project the results to the constraint conditions (negative elements are set to zero and are renormalized). During the iteration, the weights can be calculated from the last updated kernel. The IRLS internal iteration is solved by CG algorithm.


## 3.2.1. Graph Laplacian Regularization

In this section, we first outline the definition of a graph and then introduce the definition and construction method of the graph Laplacian. As depicted in Figure3, any image can be defined as an intensity function on the vertices V of a weighted graph G=(V,E,K) consisting of a finite set V of vertices (image pixels) and a finite set E⊂V×Vof edges (i,j) with the corresponding weights K(i,j). These weights measure the similarity between vertices (pixels) i and j in the graph. For the image of a n×n size, the value of its intensity function can be defined by x=x1,x2,⋯,xNT, where N=n2, and the weight of similarity K can be expressed with a symmetric nonnegative matrix of the N×N size. The graph Laplacian matrix can be derived from K.Figure 3
Graph representation of an image and construct of kernel similarity matrixK. (a) Graph representation of the image. (b) Construct of kernel similarity matrix K.
(a) (b)There are four different definitions of the graph Laplacian commonly used in the literature in the context of graph signal and image processing, each of which has different spectral characteristics [32–35]. An unnormalized Laplacian L=D-K is proposed in [33], where D is a diagonal matrix whose ith diagonal element is the sum of all elements of the ith row of matrix K. In [34], the author introduces the traditional normalized graph Laplacian L=I-D-1/2KD-1/2. In [35], another form of the graph Laplacian L=I-D-1K is given. In [32], a new normalized graph Laplacian L=I-C-1/2KC-1/2 is proposed and its characteristics are analyzed. Because of its excellent characteristics in mathematics, we use it as a constraint on the image in the cost function of our method. Then, we introduce an approach to obtaining this graph Laplacian matrix and outline its excellent features.The similarity matrix involved in the graph Laplacian is derived from the formula given in [36]. As shown in Figure 2(b), the pixel (i,j) of the similarity matrix in the image is defined as a nonlocal means (NLM)(17)Ki,j=exp⁡-x^i-x^j2h2in which x^i and x^j are patches around the pixels i and j of image x^ and h is a smoothing parameter. Owing to the similarity matrix K being a symmetric nonnegative matrix, we only calculate the similarity between each patch and its small neighborhood. Therefore, matrix K is sparse. The sparse structure is appealing from a computational point of view. Then, applying the Sinkhorn algorithm [19] to matrix K will return a diagonal scaling matrix C, such that the resulting matrix W=C-1/2KC-1/2 is a symmetric nonnegative doubly stochastic filtering matrix. Finally, the normalized graph Laplacian matrix L=I-W=I-C-1/2KC-1/2 is given. The specific process is shown in Figure 4.Figure 4
Diagram of the Laplacian matrix calculation process.It is pointed out in [32] that W is symmetric and can be decomposed as W=VSVT, where V is an orthonormal matrix whose columns are the eigenvectors of W and S=diag⁡λ1,λ2,⋯,λN is a diagonal matrix consisting of eigenvalues of W as its eigenvalues in the range of 1=λ1>λ2≥⋯≥λN≥0 [37]. It can be observed that the maximum eigenvalue is just equal to one with the corresponding DC eigenvector v1=(1/N)1,1,⋯,1T. This means that if matrix W is applied to the signal, the DC component of the signal can be protected. Spectral analysis of the matrix W shows that it is intrinsically low-pass (the largest eigenvalue corresponds to the DC component) [38], which is a very good property for the filter. The normalized Laplacian matrix L=I-W has high-pass filter characteristics (DC eigenvectors correspond to their null eigenvalues). This property is consistent with the expected behavior of Laplacian filters in image processing.The definition of the difference operator ofx:V→R of the edge (i,j)∈E in the weighted graph G is given in [32](18)dxi,j=Ki,jxjCj,j-xiCi,iwhere C(i,i) and C(j,j) are the corresponding j and ith diagonal elements of the diagonal matrix C derived from the Sinkhorn algorithm [19].The weighted gradient vector ofx at vector i∈V can be expressed as(19)∇xi=dxi,j1,⋯,dxi,jmT,∀i,j∈ETherefore, the Laplacian ofx at vector i can be written as(20)Δxi=∇2xi=1Ci,i∑j,j~iKi,jxiCi,i-xjCj,jAs a result, the regularization termR(x) can be expressed as(21)Rx=12∑i=1N∇xi2=12∑i=1N∑j,j~iKi,jxiCi,i-xjCj,j2=xTLx


## 3.2.2. Solve the Cost Function

In this step, we give the process of solving the cost equation. Based on the above deduction, we can finally rewrite the cost function expression as(22)Ex,H=I+βLy-Hx2+ηxTLx+λHFurther, the cost function is rewritten in the following form:(23)Ex,H=I+βL1/2y-Hx2+ηL1/2x2+λHwhere I+βL=I+β(I-W)=VΛVT is a symmetric and positive semidefinite matrix. Therefore, the matrix I+βL1/2=VΛ1/2VT has a filtering behavior similar to that of I+βL. By using the eigendecomposition matrix W, we can compute matrix Λ, whose ith diagonal element is 1+β(1-λi). Here, λi is the ith diagonal element of diagonal matrix S decomposed by matrix W. Since Laplacian matrix L is a high-pass filter, with β>0, I+βL1/2 corresponds to a sharpening filter for the residual y-Hx. The same analysis applies to the second term in (23), where both Laplacian matrix L and its square root L1/2 are adaptive high-pass filters. Therefore, the regularization term in (23) adaptively punishes the high frequency of the final solution. Moreover, it suppresses the noise in the final solution while maintaining the image details.For solving (23), we use the method of alternating estimates of x and H to solve for the minimum of E(x,H). Specific means of updating x and H are outlined below.(A) Update x. In this step, we fix H and obtain x as the minimization subproblem:(24)min⁡I+βL1/2y-Hx2+ηL1/2x2It is apparent that thex subproblem is a convex function. To minimize (24), the corresponding gradient is equal to zero as(25)∂E∂x=-2HTI+βLy-Hx+2ηLx=0This can be written as a symmetric positive definite system of liner equations(26)HTI+βLH+ηLx=HTI+βLyA fast computational method, such as a conjugate gradient (CG) optimization algorithm, can be used to solve the above system. Note thatH and HT are a blur kernel PSF or its reconstructed forms, respectively. The matrix I+βL must have a square root. The parameter β should maintain the above system as positive, thereby satisfying β≥-1. The value of regularization parameters η and β are selected based on the noise variance and blurring scenario. They are fixed at each step of the algorithm. In this paper, the parameter β lies in the range of (0,1) and the parameter η is empirically selected in the range of (0,0.4). The algorithm for updating the latent image x is in Algorithm 3.(B) Update H. In this step, we fix x and obtain H of the minimization subproblem:(27)min⁡I+βL1/2y-Hx2+λHBlur kernelH satisfies the condition H≥0, ∑iHi=1. In our experiments, the parameter λ depends on the user-specified kernel size h, according to the formula: λ=(3/13)h. We use the unconstrained iterated reweighted least squares (IRLS) [39] and then project the results to the constraint conditions (negative elements are set to zero and are renormalized). During the iteration, the weights can be calculated from the last updated kernel. The IRLS internal iteration is solved by CG algorithm.


## 4. Experiment

In this section, the effectiveness of our approach is verified by two groups of real database experiments. The databases are given in Table1 and Figure 5. The first databases (A1-A2) are degraded image sequences with ground truth. The second databases (B1-B3, C1-C3, and D1-D3) are degraded image sequences without ground truth, which are affected by different degrees of turbulence. These databases are image sequences that were collected at different times of the day. The focus of the image acquisition device was 2,000 mm, in which the exposure time of B1-B2, C1-C2, and D1-D2 was 4 ms, and the exposure time of B3, C3, and D3 was 8 ms. In image quality evaluation, SROCC (Spearman rank-order correlation coefficient), KROCC (Kendall rank-order correlation coefficient), and PLCC (Pearson linear correlation coefficient) are used to measure the correlation between the objective value and the subjective value. The higher the correlation between the objective value and the subjective value, the closer the value of aforementioned three correlation coefficients to 1. The closer the value is to 1, the more accurate the algorithm is.Figure 5
Degraded sequences. (a) Real dataset with ground truth. (b)-(d) Real dataset without ground truth.
(a) (b) (c) (d)(A) Real Dataset with Ground Truth. The degraded sequencesChimney andBuilding with ground truth have been made available with their ground truth by Hirsch [10]. The sequences were also used in Zhu’s method [18] and N. Anantrasirichai’s method [17]. We employed these images to verify the effectiveness of algorithm. The experimental results and objective evaluation are given in Figures 6 and 7, respectively, and they are compared with the algorithm [10, 17, 18]. Lin Zhang’s method [40] compares several full-reference image quality evaluation (FRQE) algorithms. The correlation coefficients SROCC, KROCC, and PLCC between the objective value and the subjective value of each algorithm are given. The correlation coefficients of FSIM [41], IW-SSIM [42], RFSIM [43], and MS-SSIM [44] are closer to 1, which indicates that the algorithm has high accuracy. In addition, we also used PSIM [45] and ADD-GSIM [46] to evaluate the restoration results. The PSNR and MS-SSIM of the algorithm [10, 17, 18] in Figure 7 are derived from algorithm [17]. PSNR, MS-SSIM, IW-SSIM, RFSIM, and FSIM values reveal that our method outperforms other methods and the performance of PSIM is similar to ADD-GSIM. Interestingly, the subjective results shown in Figure 6 reveal that our results are the most similar to the ground truth.Figure 6
Restored result of A1-A3. (a) Blurred image. (b) Ground truth. (c) Algorithm [10] results. (d) Algorithm [18] results. (e) Algorithm [17] results. (f) Proposed algorithm results.
(a) (b) (c) (d) (e) (f)Figure 7
Comparison of objective results using FRQE methods. (a) A1. (b) A2.
(a) (b)(B) Real Datasets without Ground Truth. Nine sequences, Figure 5 B1-D3, that exhibit significant turbulence distortions are employed in this trial. B1-B3, C1-C3, and D1-D3 are three groups of degraded sequences affected by varying degrees of turbulence. These degraded sequences were restored with our proposed method and algorithm [17], and the specific results are shown in Figures 8, 9, and 10, respectively.Figure 8
Restoration results of one-degree turbulence-degraded sequence image. (a) Degraded sequence. (b) Algorithm [17] results. (c) Proposed Algorithm results.
(a) (b) (c)Figure 9
Restoration results of two-degree turbulence-degraded sequence images. (a) Degraded sequence. (b) Algorithm [17] results. (c) Proposed algorithm results.
(a) (b) (c)Figure 10
Restoration results of three-degree turbulence-degraded sequence images. (a) Degraded sequence. (b) Algorithm [17] results. (c) Proposed Algorithm results.
(a) (b) (c)The correlation coefficients SROCC, KROCC, and PLCC are given in [47–49]. They indicate the correlation between objective results BRISQURE [47], NIQE [48], and SSEQ [49] and the subjective values. In addition, we also used NFERM [50] and BIQI [51] to evaluate the restoration results. Tables 2, 3, and 4 give objective evaluation results of three groups of restored images, respectively. These data indicate that the BRISQURE, NIQE, and SSEQ values of this algorithm are higher than the value of the compared algorithm [17]. The NFERM values of B3-D3 reveal that our method outperforms algorithm [17]; however, the NFERM value of B1-B2 is highest for algorithm [17]. The results of algorithm [17] show that noise still exists. The noise forms a strong boundary between the target and the background.Table 2
Objective results of restored sequences on one-degree turbulence.
Method BIQI BRISQURE NIQE SSEQ NFERM B1 [17] 85.9318 7.4147 7.2710 17.3282 48.7388 Our 29.7874 34.1766 8.7277 35.4349 16.1062  B2 [17] 91.1420 5.4071 7.6322 3.2639 42.3259 Our 66.4207 31.2737 9.1349 21.4994 14.0716  B3 [17] 88.5376 4.0857 7.8358 1.3358 23.4485 Our 46.8050 50.3921 8.9871 28.6650 40.5583Table 3
Objective results of restored sequences on two-degree turbulence.
Method BIQI BRISQURE NIQE SSEQ NFERM C1 [17] 92.6088 10.7997 6.6298 16.6803 35.8362 Our 40.4595 49.2980 8.3991 42.4510 58.5245  C2 [17] 71.3531 2.0506 7.1338 0.5063 35.6659 Our 78.8264 25.4969 9.3054 16.4931 39.0322  C3 [17] 78.1963 2.3538 8.0516 1.6752 27.5387 Our 59.2097 32.9010 10.5940 1.9928 51.1725Table 4
Objective results of restored sequences on three-degree turbulence.
Method BIQI BRISQURE NIQE SSEQ NFERM D1 [17] 78.8978 1.0422 9.5785 10.4626 36.2913 Our 73.5125 76.7459 12.9764 53.6843 72.5698  D2 [17] 77.7949 3.1257 10.4056 3.6847 30.4773 Our 64.9152 58.5710 11.5842 41.7042 57.0911  D3 [17] 76.7288 5.6604 9.9089 2.8860 32.7111 Our 52.2458 83.8020 12.6648 46.2560 61.0262


## 5. Discussion and Conclusion

It remains a difficult problem to restore a latent image from a turbulence-degraded sequence. In this paper, a method is presented to restore a high-quality image from a degraded sequence image caused by atmospheric turbulence. An optical flow technique is used to register sequence images for eliminating distortion of degraded images. The image can be constrained by constructing the graph Laplacian regularization term so that the restored image can maintain a considerable amount of detail and texture. Compared with the previous method, our method demonstrates both subjective and objective improvement. However, our algorithm uses the optical flow method to register the image and the time cost is thus slightly higher. Further research can be considered for reducing the time complexity.



---

*Source: 4970907-2018-07-19.xml*