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

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 themethod can effectively eliminate the distortion and blur, restore the image details, and significantly improve the image quality.


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][2][3][4][5][6] and image restoration technology [7][8][9][10][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 socalled "lucky imaging technique" [12][13][14][15][16][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, pixelbased 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:  = (()) + , where  and  denote the degraded image and the latent image, respectively, and  is the noise.We assume that  is the geometric distortion factor and  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-ofthe-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 Section 2. 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.

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: The 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,  is the original static scene,   is captured image at time ,  is a space and time invariant blur kernel,   is the space and time dependent deformation matrix at time , and   is white Gaussian noise at time .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 The deformation matrix   describes the displacement of each pixel between image   and image .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   by image registration and then perform an inverse mapping on   so that the registered sequence images are only affected by the blur kernel .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  is unknown, its sparse feature is regarded as a regularization term, and the latent sharp image  is obtained by solving the reconstructed cost function.

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 sequence {}  .
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.

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.
Require: Observe all frames of degraded image {}  and reference image  1. Minimize the cost function to obtain the displacement map u.
(ii) Solve equation ( 12) and update k.Moreover, it has strong robustness to illumination changes, shading, and noise enhancement.
Two images frames,  0 and  1 : (Ω ⊆ R 2 ) → R, are given.The objective is to find the disparity map u : Ω → R 2 , 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: where ( 0 (x) −  1 (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 () = || and (∇u) = |∇u|; hence, the energy equation can be expressed as Since 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 (x) to represent the disparity map u(x) and  + (x) to x + ((x), 0)  .We linearize the image  1 near x +  0 : where   1 is the derivative of the image intensity  1 in the  direction.Fixing  0 and employing the linear approximation for  1 , the TV-L1 equation can be written as We define the current residual  1 (x +  0 ) + ( −  0 )  1 −  0 as (,  0 , x), which is abbreviated as ().Moreover, we introduce an auxiliary variable V and propose minimizing the convex approximation of the above equation: where  is a small constant and V is an approximation of .This convex minimization problem can be solved by the alternating iterative updating of  and V.
The equation generalization to more dimensions is the following energy: Similar to the one-dimensional case, the method of alternating updating is used to solve the problem.
For every  and fixed V  , we solve For u being fixed, we solve The disparity map u 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 Algorithm 2, we obtain the registered sequences {}  and remove the geometric distortion factor   .The specific experimental results are shown in Figure 2, and the specific information of the image is given in Table 1.i j

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 image  we observe comes from a sharp image  blurred by blur kernel ℎ along with the addition of Gaussian noise , which can be expressed as Our goal is to recover the unknown image  and the blurring kernel ℎ.If , , and  are represented as lexicographically ordered vectors x, y, and n, the model can be rewritten as where  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 Here, we introduce the normalized graph Laplacian filter proposed in [32] to construct the regularization term.Meanwhile, we introduce the sparsity of blur kernel  to reduce the noise on the blur kernel.The constraint (nonnegative, sum of one) of  obeys the physical mechanism of blur formation.Finally, we structure the cost function as follows: where y and x represent the blurred image and latent image, respectively. represents the blur factor, which is the reconstruction of PSF.The first term of the cost function is the data fidelity, where () is a function of the graph Laplacian matrix, and () = ( + ).The second term (x) is the constraint of the latent image x, which is related to the graph Laplacian.The third term () is the constraint of the blur kernel .To reduce noise in the blurring kernel, we use the L1 regularization as its constraint and the expression is () = ‖‖.This constraint on  (sumto-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 (x) will be derived next.

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 Figure 3, any image can be defined as an intensity function on the vertices  of a weighted graph  = (, , ) consisting of a finite set  of vertices (image pixels) and a finite set  ⊂ ×of edges (, ) with the corresponding weights (, ).These weights measure the similarity between vertices (pixels)  and  in the graph.For the image of a  ×  size, the value of its intensity function can be defined by x = [(1), (2), ⋅ ⋅ ⋅ , ()]  , where  =  2 , and the weight of similarity  can be expressed with a symmetric nonnegative matrix of the  ×  size.The graph Laplacian matrix can be derived from .
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][33][34][35].An unnormalized Laplacian  =  −  is proposed in [33], where  is a diagonal matrix whose th diagonal element is the sum of all elements of the th row of matrix .In [34], the author introduces the traditional normalized graph Laplacian  =  −  −1/2  −1/2 .In [35], another form of the graph Laplacian  = − −1  is given.In [32], a new normalized graph Laplacian  =  −  −1/2  −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 (, ) of the similarity matrix in the image is defined as a nonlocal means (NLM) in which x and x are patches around the pixels  and  of image x and ℎ is a smoothing parameter.Owing to the similarity matrix  being a symmetric nonnegative matrix, we only calculate the similarity between each patch and its small neighborhood.Therefore, matrix  is sparse.The sparse structure is appealing from a computational point of view.Then, applying the Sinkhorn algorithm [19] to matrix  will return a diagonal scaling matrix , such that the resulting matrix  =  −1/2  −1/2 is a symmetric nonnegative doubly stochastic filtering matrix.Finally, the normalized graph Laplacian matrix  =  −  =  −  −1/2  −1/2 is given.The specific process is shown in Figure 4.
It is pointed out in [32] that  is symmetric and can be decomposed as  =   , where  is an orthonormal matrix whose columns are the eigenvectors of  and  = diag{ 1 ,  2 , ⋅ ⋅ ⋅ ,   } is a diagonal matrix consisting of eigenvalues of  as its eigenvalues in the range of 1 =  1 >  2 ≥ ⋅⋅⋅ ≥   ≥ 0 [37].It can be observed that the maximum eigenvalue is just equal to one with the corresponding DC eigenvector This means that if matrix  is applied to the signal, the DC component of the signal can be protected.Spectral analysis of the matrix  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  = − 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 of  :  → R of the edge (, ) ∈  in the weighted graph  is given in [32]  (, ) = √  (, ) (  () where (, ) and (, ) are the corresponding  and th diagonal elements of the diagonal matrix  derived from the Sinkhorn algorithm [19].
The weighted gradient vector of  at vector  ∈  can be expressed as Therefore, the Laplacian of  at vector  can be written as As a result, the regularization term (x) can be expressed as 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 Further, the cost function is rewritten in the following form: where  +  =  + ( − ) = Λ  is a symmetric and positive semidefinite matrix.Therefore, the matrix ( + ) 1/2 = Λ 1/2   has a filtering behavior similar to that of  + .By using the eigendecomposition matrix , we can compute matrix Λ, whose th diagonal element is 1 + (1 −   ).Here,   is the th diagonal element of diagonal matrix  decomposed by matrix .Since Laplacian matrix  is a high-pass filter, with  > 0, ( + ) 1/2 corresponds to a sharpening filter for the residual y − x.The same analysis applies to the second term in (23), where both Laplacian matrix  and its square root  1/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  to solve for the minimum of (x, ).Specific means of updating x and  are outlined below.
(A) Update x.In this step, we fix  and obtain x as the minimization subproblem: It is apparent that the x subproblem is a convex function.To minimize (24), the corresponding gradient is equal to zero as This can be written as a symmetric positive definite system of liner equations A fast computational method, such as a conjugate gradient (CG) optimization algorithm, can be used to solve the Require: Blur kernel  from the previous  update, Image x from previous x update Require: Regularization parameter .(i) Compute  from x using equation (17) (ii) Apply Sinkhorn algorithm in [19] to  to obtain diagonal matrix .(iii) Compute filtering matrix  and graph Laplace matrix .(iv) Solve cost function using CG to compute x. return update image x Algorithm 3: Update sharp image x.
above system.Note that  and   are a blur kernel PSF or its reconstructed forms, respectively.The matrix  +  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 .In this step, we fix x and obtain  of the minimization subproblem: Blur kernel  satisfies the condition  ≥ 0, ∑    = 1.In our experiments, the parameter  depends on the user-specified kernel size ℎ, according to the formula:  = (3/13)ℎ.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.

Experiment
In this section, the effectiveness of our approach is verified by two groups of real database experiments.The databases are given in Table 1 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 rankorder 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.
(A) Real Dataset with Ground Truth.The degraded sequences Chimney and Building 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 fullreference 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.
(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.
The correlation coefficients SROCC, KROCC, and PLCC are given in [47][48][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.

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.

Figure 3 :
Figure 3: Graph representation of an image and construct of kernel similarity matrix .(a) Graph representation of the image.(b) Construct of kernel similarity matrix .

Figure 4 :
Figure 4: Diagram of the Laplacian matrix calculation process.
Observe all frames of blurry image {}  , reference image , and maximum kernel size ℎ, maximum iteration number . 1.According to the steps of Algorithm 2 for image registration, obtain distortion factor   of each image frame and the registered sequences {}  .2. Determine the median   of each frame in the registered sequence images {}  by the temporal filter.3. Conduct blind estimation of blur matrix  from blurred image   .
Require:Figure 1: Block diagram of the proposed image restoration method.

Table 2 :
Objective results of restored sequences on one-degree turbulence.

Table 3 :
Objective results of restored sequences on two-degree turbulence.

Table 4 :
Objective results of restored sequences on three-degree turbulence.