Turbulent image deblurring using a deblurred blur kernel

In the context of addressing a noisy turbulence-degraded image, it is common to use a denoising low-pass filter before implementing a deblurring algorithm. However, this filter not only suppresses noise but also induces a certain degree of blur into the degraded image. This blur effect causes a blurred estimate of the true blur kernel and ultimately leads to a distorted estimate of the latent clear image. To tackle this issue, this paper presents an innovative single-image deblurring method. It integrates a dedicated blur kernel deblurring step to mitigate the effects of the denoising filter. The L 0 norm and L 2 norm serve as the respective constraints for latent clear image and blur kernel. Experimental results on both synthetic and real-world turbulence-degraded images demonstrate the effectiveness and efficiency of the proposed method.


Introduction
Image deblurring is a crucial tool to improve atmospheric turbulence-degraded image quality.The aim of blind image deblurring is to estimate the latent clear image o and the turbulent blur kernel h with only the degraded image g.The degradation process is typically modeled as Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where * denotes the 2D convolution operator and n represents the additional noise.The maximum-a-posterior (MAP) framework is widely used to solve this problem [1][2][3][4][5], where the key lies in accurately estimating blur kernels [6,7].These methods involve two steps.In the initial step, they aim to find o and h by maximizing the posterior p(o, h|g).Assuming independent Gaussian noise, the problem is formulated as where λ and γ are weight parameters; ϕ(o) and φ(h) represent the prior terms (also known as constraint terms) for the latent clear image and blur kernel, respectively.In the subsequent step, a non-blind image restoration algorithm is implemented to recover the ultimate latent clear image, using the blur kernel estimated by solving (2).To mitigate noise in the degraded image, a Gaussian filter is commonly applied.Thus the problem of equation ( 2) is transformed into where G denotes the discrete Gaussian filter.Zhong et al [8] have pointed out that the use of the Gaussian filter leads to a biased kernel estimation, subsequently affecting the accuracy of the final latent clear image estimation.When neglecting the prior term of the kernel, the optimization of (3) results in the following kernel estimation ≈ arg min where h represents the actual blur kernel for the original degraded image g.To solve this issue, Zhong et al [8] utilized a series of directional low-pass filters and inverse radon transforms to estimate the blur kernel.However, this method is more suitable for motion deblurring.In addition, this approach involves multiple image deconvolutions and inverse radon transforms, which significantly increases the overall algorithm time.
In this paper, we introduce a new strategy for accurately estimating turbulent blur kernels in turbulence-degraded image deblurring.In contrast to [8], our approach only necessitates a denoising filter and a pseudo Wiener filter for deblurring the blur kernel, effectively avoiding the need for baseline image restoration across a series of directional low-pass filtered blurry images.Instead of the conventional L α gradient prior [2,8,9], we adopt an unnatural L 0 gradient prior [10,11] to encourage sparser and more localized representations.In comparison to other recently proposed regularization priors [4,[12][13][14][15][16], our simplified L 0 gradient prior excels in promoting high sparsity while ensuring rapid acquisition of solutions in our turbulent image restoration algorithm.To regularize the kernel, we utilize the L 2 norm prior [10,13], enabling the fast computation of involved subproblems through the use of the fast Fourier transform (FFT).Additionally, we incorporate a common multi-scale strategy [2,4,8,10,12,13,17] into our approach.The experimental results show that the proposed approach can effectively and efficiently improve the quality of the estimated turbulent blur kernel and latent clear image.
The main contributions of this paper include the following: (1) A novel strategy aimed at addressing the problem of biased kernels induced by denoising filters is presented; (2) The utilization of a pseudo Wiener filter is proposed to estimate a clear kernel from its blur version; (3) The image and kernel are constrained using L 0 gradient norm and L 2 norm, respectively; (4) A multi-scale restoration strategy is also deployed to achieve robust kernel estimation; (5) Building upon the aforementioned contributions, we present a novel approach for deblurring turbulencedegraded images.

Related work
Blind image deblurring is an ill-posed problem, which has infinite paired solutions.In order to constrain the problem to obtain a clear solution, a variety of image priors have been proposed.The initial prior gained popularity was the total variation prior [18].This prior aimed at promoting sparsity in image gradients, thus improving the overall image quality.Subsequently, researchers, such as Krishnan and Fergus [19] and Babacan et al [9] adopted the hyper-Laplacian function as an image prior due to its strong ability to enforce sparsity constraints.Others explored diverse functions as image priors to approximate the heavy-tailed distribution of natural image gradients, including a Gaussian function [20], a mixture of Gaussians [21,22] and a piecewise-defined function [23].
Particularly, a fractional-order total variation is proposed to more precisely and delicately capture image textures in image deblurring [24].However, Levin et al [25] stated that these priors tend to favor blur solutions over clear ones.Afterward, various novel image priors were developed to favor clear solutions, including L 1 /L 2 prior [1], sound L 0 prior [10,26], nonzero constraint prior [27], nonconvex L 1 − αL 2 prior [15,28], and saturation-value geo-metric spatial-feature prior [16].Leveraging the sparsity advantage of L 0 prior and an efficient optimization framework [10], a series of L 0 + X style priors have been introduced [17], including dark channel prior [4], extreme channels prior [12], and local minimal intensity prior [13], and enhanced sparse prior [14].However, the incorporation of X prior, motivated differently from L 0 , tends to markedly diminish the algorithm's efficiency [17].An alternative strategy to avoid blur solutions is extracting edge images from the original degraded image and utilizing them for blur kernel estimation [29,30].The efficacy of these methods is highly contingent on the performance of the edge filter.
Recently, deep learning-based methods have achieved remarkable strides in visual processing tasks.Diverse network architectures have been also employed for image deblurring, encompassing convolutional neural networks (CNNs) [31,32], recurrent neural networks [33], end-to-end conditional generative adversarial networks [34,35], multi-scale residual networks [36], and U-net like models [37][38][39].Particularly, Cai et al [32] embedded the dark and bright channels prior into a multi-scale CNN for motion image deblurring.These models excel in capturing intricate structures and details within an image, potentially acquiring a more sophisticated representation of prior knowledge related to clear images, surpassing conventional methods.While these methods demonstrate outstanding performance on their respective test datasets, their heavy reliance on training data features makes them more suitable for particular scenarios, consequently restricting their overall generality.

Proposed method
The L 0 norm counts the number of nonzero pixel values in the image gradient domain.Formulated mathematically to encourage high sparsity in regularization, it effectively promotes both sparsity and rapid convergence.Figure 1 illustrates the horizontal and vertical gradient histograms of a clear image and a corresponding turbulence-degraded counterpart.The density of nonzero values in turbulent image gradients surpasses that of clear image gradients, and the frequencies of nonzero values in the turbulent image are higher than those in the clear image.Thus we utilized the L 0 norm as the prior of the image gradients.
We employ L 2 norm as the kernel constraint term and impose value constraints for the image and kernel following the methods of [2,18].With these considerations, our objective function becomes We utilize the alternating minimization (AM) [18,40] to optimize this problem.Given the blur kernel, the optimization problem of the latent clear image becomes and the optimization problem of the turbulent blur kernel becomes where ∇ denotes the gradient operator.Assuming g(p, q) represents the grayscale value of a pixel at coordinates (p, q) in the degraded image, and G(p, q) represents the value after applying Gaussian filtering G at the same coordinates.The Gaussian filter is defined as follows where k represents the radius of the filter, and x and y denote the horizontal and vertical distances from the point (p, q), respectively.The parameter σ x and σ y control the filtering strength along the respective horizontal and vertical directions.The filter size is determined by truncating the cumulative distribution function of the Gaussian kernel at a threshold of 0.05.As analyzed in equation ( 4), the optimization of equation ( 7) will produce a blurred version of the actual kernel.We employ a pseudo Wiener filter to offset this effect to recover a clear kernel from the blurred one.
In accordance with [2,18], we first optimize the minimization problems of equations ( 6) and ( 7) without incorporating value constraints and then enforce the value constraints in a subsequent step.The half-quadratic splitting L 0 minimization method [41] is employed to optimize the problem outlined in equation (6).Introducing an auxiliary variable u corresponding to the gradient of o, we can rewrite equation (6) as where β is a weight parameter.This equation is equal to equation ( 6) without the value constraint when β gets close to ∞.It can be solved by alternatively minimizing o and u in each iteration.Fixing u, the problem of o is This is a least squares minimization problem, and its solution is computed via where F(•) and F −1 (•) represent the FFT and inverse FFT operators, respectively; F(•) denotes the complex conjugate operator; , where ∇ x and ∇ y are the horizontal and vertical differential operators, respectively.
Fixing o, the problem of u is This is a pixel-wise minimization problem, and the solution is given by Subsequently, we impose the value constraint of o via The main steps for solving the latent clear image o are outlined in algorithm 1. Gaussian filter G, noise-to-signal-power ratio δ, and parameter γ.Solve for h using equation (16).Update h sequentially with equations ( 17)- (19).Output: Turbulent blur kernel h.
Without the value constraints, the problem of equation ( 7) is also formulated as a least squares minimization problem.Similar to the approach presented in [4,26], we estimate the blur kernel in the gradient domain, i.e.
The solution is computed by FFTs [42] via Then, we enforce the value constraint of h g via As analyzed in equation ( 4), this solution is a blurred version of the actual turbulent kernel.To obtain a clear estimation and mitigate the effect of the Gaussian denoising filter, we apply a pseudo Wiener filter, i.e.
where G is the discrete presence of the Gaussian filter, and δ is a parameter representing the noise-to-signal-power ratio in the kernel.The ultimate estimate of the turbulent kernel is updated after enforcing the value constraint again with Equation ( 18) improves the accuracy of the estimated blur kernel.The value constraints of equation ( 19) ensure the validity and reasonableness of the final estimated blur kernel.The main steps for solving equation (15) are outlined in algorithm 2.
With the estimated turbulent kernel, a non-blind image restoration method from [26] is applied to retrieve the ultimate clear image.Similar to other existing methods, we also incorporate a multiscale strategy into our deblurring algorithm.The main steps of the overall algorithm are summarized in algorithm 3.

Algorithm 3. Multiscale turbulent image deblurring.
Input: Turbulence-degraded image g with size m × n, kernel size s × s, weight parameters λ, γ and βmax, noise-to-signal-power ratio δ, parameter of Gaussian filter r, and number of iteration T.

end if if l == 1 then
Initialize h with an s l × s l zeros matrix adding [0.5,0.5] in center.
Initialize o with g. else Upsample h to size of s l × s l using bilinear interpolation method.end if for t = 1 to T do Solve for o using algorithm 1. Solve for h using algorithm 2. Update weight parameter λ as λ = max(λ/1.1, 1 * 10 −4 ).end for end for Solve for o using non-blind image restoration method in [26].Output: Deblurred image o and turbulent blur kernel h.
Several hyperparameter settings are involved in applying our algorithm to specific images.However, we have empirically found that the following parameter settings consistently yield satisfactory recovery results for different images: λ = 4 * 10 −3 , γ = 2, β max = 1 * 10 5 , δ = 0.1, r = 1 and T = 5.The parameter settings for the image and kernel prior constraints λ, γ, and β max refer to the settings in [26].It is generally recommended to set δ to 0.1, with larger values suggested for very complex turbulent blurring situations.The Gaussian denoising parameter r is typically set in the range of 1 to 4, with larger values recommended for noisier images.Figure 2 illustrates the objective loss function curve with respect to iterations at the image scale in the multiscale strategy approach.The algorithm converges around the 5th iteration.Therefore, we recommend setting T = 5 as a compromise between achieving sufficient accuracy and ensuring computational efficiency.

Experimental results
We performed evaluations of our proposed and four stateof-the-art MAP deblurring algorithms [1,2,4,13] on both synthetic and real-world turbulence-degraded images for comparison.The algorithms were implemented in Matlab R2014a on a PC with an Intel Core i5 CPU and 8GB of RAM.

Synthetic data
We employed the Von Kármán statistic phase screen model [43] to generate turbulent blur kernels.The random phase screen function was generated by inverting the Von Kármán power spectrum model [44].By adjusting the values of D/r 0 around 4, where D represents the aperture diameter and r 0 is the coherence diameter of the atmosphere, we generated three turbulent blur kernels, each corresponding to a different level of turbulence.Given that the evaluated algorithms required the size of the blur kernel as a prior, we conducted boundary cropping on each kernel to achieve an appropriate size.In detail, we removed boundaries smaller than 1/255 of the maximum value within each corresponding blur kernel.Subsequently, we added zero rows or columns to ensure that each kernel image had equal and odd length and width.Eventually, the kernels were normalized to 1. Figure 3 presents the three synthetic turbulent blur kernels.
The turbulence-degraded images were simulated by convolving two clear images with the three synthetic blur kernels and then adding Gaussian noise with a mean of zero and a standard deviation of 0.01.Before applying convolution, all images were normalized to the range of [0, 1].The simulated turbulence-degraded images are depicted in figure 4.
We performed deblurring experiments on these six degraded images.The visual results are presented in figure 5.For our algorithm, we set the parameters as follows: λ = 4 * 10 −3 , γ = 2, β max = 1 * 10 5 , δ = 0.1, r = 1, and T = 5.The other four methods were fine-tuned according to the reference parameters provided by the authors to achieve the best results.To ensure a fair comparison, all algorithms adopt the non-blind image deblurring method from [26], where the weights of the Laplacian prior and L 0 prior are set to 5 * 10 −3 and 5 * 10 −4 , respectively.As shown in figure 5, Our method can estimate more accurate turbulent blur kernels and eventually recover clearer images.These results demonstrate the effectiveness of our blur kernel deblurring strategy.
We employed normalized mean square error (NMSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM) [45] to evaluate the results.The NMSE was used to measure the accuracy of the estimated blur kernels.The PSNR and SSIM were used to evaluate the quality of the estimated images.The evaluation results are presented in tables 1 and 2. Of all the methods, our results for each image were consistently among the highest.Moreover, our method achieves the best average performance in all three metrics.This evidences the superiority of our method in estimating turbulent blur kernel and improving turbulent image quality.Figure 6 illustrates the intermediate estimates of our kernel (figure 5(d)) at each scale.Our kernel deblurring step can efficiently refine the estimation of blur kernels.Visual results on synthetic turbulence-degraded images.Left to right: deblurred images using the method of Krishnan et al [1], method of Perrone and Favaro [2], method of Pan et al [4], method of Wen et al [13], and our proposed method.
Table 1.Quantitative results on the blur kernels: NMSE.

Methods
Figure 5  The bold values highlight the best results for each image.
In the above experiments, the kernel size is assumed to be a known prior.In practical situations, this prior is often unavailable.Many blind deblurring methods opt to use a large kernel size to ensure the support domain.However, using oversized kernels can introduce errors in kernel estimation, ultimately resulting in pronounced artifacts in the deblurring results [6].We repeated the deblurring experiments on the synthetic images, this time using larger blur kernel sizes.Precisely, all kernel sizes were set to 3/2 times that of their ground truths.The deblurring results are shown in figure 7.All results are degraded compared to figure 5, but our results show the least degradation and achieve the highest visual performance.
We also applied PSNR and SSIM to assess the quality of the deblurred results.The averaged results across all images are depicted in table 3. Our method outperforms all other methods.This demonstrates the robustness of our method in dealing with oversized kernels.The bold values highlight the best results for each image in each evaluation criterion.

Real data
The real-world turbulence-degraded test data includes three images of a pavilion, sunspots, and an international space station (ISS).The pavilion image and sunspots image were acquired by our teams, while the ISS image was sourced from www.tracking-station.de/images/images.html.We conducted deblurring experiments on these three images to evaluate the performance of different deblurring algorithms in real-world scenarios.The results are presented in figures 8-10.Due to the high noise level in figures 8 and 10, we increased the Gaussian filter parameter r to 2. For figure 10, we utilized the default non-blind deconvolution algorithm.For figure 9, we employed the non-blind deconvolution algorithm for partially saturated images from [46].The results demonstrate that our deblurred images contain more high-frequency details and fewer artifacts than the other four competing methods.
We recorded the time taken by each method to process the three images.The details are given in table 4. The time required for executing the non-blind deconvolution algorithms is not included in the measurements.Our method surpasses all the others in terms of speed, and this advantage becomes even more significant when handling large images.This demonstrates the exceptional efficiency of our method.Results on synthetic turbulence-degraded images with larger kernel size settings.Left to right: deblurred images using the method of Krishnan et al [1], method of Perrone and Favaro [2], method of Pan et al [4], method of Wen et al [13], and our proposed method.
Table 3. Quantitative evaluation on the synthetic turbulence-degraded images with larger kernel size settings: average PSNR, average SSIM.

Methods
Average PSNR Average SSIM Krishnan et al [1] 22.16 0.7673 Perrone and Favaro [2] 18.94 0.5891 Pan et al [4] 13.93 0.4247 Wen et al [13] 23.81 0.8218 Ours 24.53 0.8295 Results on a real-world image of a pavilion.Left to right: degraded image, results using the method of Krishnan et al [1], method of Perrone and Favaro [2], method of Pan et al [4], method of Wen et al [13], and our proposed method.Results on a real-world image of sunspots.Left-top to right-bottom: degraded image, results using the method of Krishnan et al [1], method of Perrone and Favaro [2], method of Pan et al [4], method of Wen et al [13], and our proposed method.Results on a real-world image of an international space station (ISS).Left to right: degraded image, results using the method of Krishnan et al [1], method of Perrone and Favaro [2], method of Pan et al [4], method of Wen et al [13], and our proposed method.

Conclusion
In this study, we developed a multiscale deblurring method for turbulence-degraded images.Benefiting from the proposed kernel deblurring step, our method is able to estimate more accurate kernels and clearer images.The experimental results on synthetic and real-world degraded images demonstrate the superiority of our method in terms of effectiveness and efficiency.The proposed kernel deblurring strategy holds potential applications on various iterative deblurring algorithms.

Figure 1 .
Figure 1.Histograms of a clear image and a corresponding turbulence-degraded image.In the top row, from left to right: the clear image, the histogram of horizontal gradients of the clear image, and the histogram of vertical gradients of the clear image.In the bottom row, from left to right: the turbulence-degraded image, the histogram of horizontal gradients of the turbulence-degraded image, and the histogram of vertical gradients of the turbulence-degraded image.

Figure 2 .
Figure 2. Objective loss function curve with iterations.The blue and red curves on the left correspond to the two examples depicted on the right, respectively.The red line represents the loss adjusted by subtracting a constant value to improve visibility.The algorithm converges at the 5th iteration.

Figure 4 .
Figure 4. Synthetic turbulence-degraded images: (a) the original image of a satellite; (e) the original image of the cameraman; (b) and (f) degraded images using blur kernel 1; (c) and (g) degraded images using blur kernel 2; (d) and (h) degraded images using blur kernel 3.

Figure 5 .
Figure 5.Visual results on synthetic turbulence-degraded images.Left to right: deblurred images using the method of Krishnan et al[1], method of Perrone and Favaro[2], method of Pan et al[4], method of Wen et al[13], and our proposed method.

Figure 6 .
Figure 6.Estimated kernels of our multiscale deblurring algorithm at different scales.Top: estimated kernels before applying kernel deblurring step.Bottom: estimated kernels after applying kernel deblurring step.

Figure 7 .
Figure 7.Results on synthetic turbulence-degraded images with larger kernel size settings.Left to right: deblurred images using the method of Krishnan et al[1], method of Perrone and Favaro[2], method of Pan et al[4], method of Wen et al[13], and our proposed method.

Figure 8 .
Figure 8. Results on a real-world image of a pavilion.Left to right: degraded image, results using the method of Krishnan et al[1], method of Perrone and Favaro[2], method of Pan et al[4], method of Wen et al[13], and our proposed method.

Figure 9 .
Figure 9.Results on a real-world image of sunspots.Left-top to right-bottom: degraded image, results using the method of Krishnan et al[1], method of Perrone and Favaro[2], method of Pan et al[4], method of Wen et al[13], and our proposed method.

Figure 10 .
Figure10.Results on a real-world image of an international space station (ISS).Left to right: degraded image, results using the method of Krishnan et al[1], method of Perrone and Favaro[2], method of Pan et al[4], method of Wen et al[13], and our proposed method.

Table 4 .
Processing time comparison on real-world turbulence-degraded images.