Deblurring Turbulent Images via Maximizing L1 Regularization

Atmospheric turbulence significantly degrades image quality. A blind image deblurring algorithm is needed, and a favorable image prior is the key to solving this problem. However, the general sparse priors support blurry images instead of explicit images, so the details of the restored images are lost. The recently developed priors are non-convex, resulting in complex and heuristic optimization. To handle these problems, we first propose a convex image prior; namely, maximizing L1 regularization (ML1). Benefiting from the symmetrybetween ML1 and L1 regularization, the ML1 supports clear images and preserves the image edges better. Then, a novel soft suppression strategy is designed for the deblurring algorithm to inhibit artifacts. A coarse-to-fine scheme and a non-blind algorithm are also constructed. For qualitative comparison, a turbulent blur dataset is built. Experiments on this dataset and real images demonstrate that the proposed method is superior to other state-of-the-art methods in blindly recovering turbulent images.


Introduction
As an essential tool for improving image quality, image deblurring has received considerable attention. Blind image deblurring involves estimating a latent sharp image o and a blur kernel h, with a blurry input image g. Mathematically, the blurring process is modeled as: where * represents the convolution operator and n denotes the additional noise. Involving an infinite number of solutions, this problem is strongly ill-posed.
This paper proposes a novel image deblurring algorithm that achieves competitive results on atmospheric turbulent images. Inspired by the success in previous studies where distinct properties were observed of specific images, such as text images [12], face images [13], and lowlight images [14], our work derived from the effects of turbulence degradation. We first analyze how various types of image priors change with turbulent blur to illustrate the strengths of the proposed prior. We use different sizes and types of 2D signals to verify the validity of the prior. We then present a new deblurring algorithm, including a soft suppression strategy, which effectively inhibits artifacts. We iteratively estimate the latent image and the blur kernel using the alternating minimization (AM) method. We further establish a turbulence blur image dataset to evaluate the effectiveness of the proposed method. The experiments on the dataset and authentic images demonstrate that our method achieves state-of-the-art performance compared to recent blind deblurring algorithms.
The contributions of this paper are as follows: (1) We present the maximizing L1 regularization (ML1) prior for latent sharp images. We theoretically illustrate that this prior is better than other priors for deblurring turbulence blur images. (2) Based on the ML1 prior, we propose a deblurring algorithm with a soft suppression strategy, termed suppressed projected alternating minimization (SPAM). (3) To assess the effectiveness of the proposed method, we build a turbulent blur images dataset. (4) Quantitative and qualitative experiments illustrate that our algorithm is superior to state-of-the-art algorithms.

Related Work
In recent years, blind image deblurring has significantly improved due to the use of various priors on images and blur kernels. Most works are based on statistical priors and gradient sparsity priors.
Fergus et al. [1] adopted a Gaussian mixture model to learn an image gradient prior through variational Bayesian inference. Levin et al. [2] explored the shortcomings of the naïve MAP x,k method and presented a more robust MAP k algorithm. Babacan et al. [6] proposed a general approach with super-Gaussian sparse image priors. However, these methods are computationally expensive.
Efficient approaches based on gradient sparse priors have also been developed. Chan and Wong [3] proposed total variation (TV) regularization, which achieves good performance on image edges. Since then, many different sparse regularizers [4,5,15,16] have been adopted to support the gradient sparsity of natural images. However, Levin et al. [11] showed that these gradient sparsity priors favor blurry images rather than sharp ones.
To support sharp images, various other image priors have been designed. Krishnan et al. [7] proposed the L1/L2 regularization function and developed a fast algorithm to solve this non-convex problem. Xu et al. [8] used the L0 sparse gradient prior for motion deblurring. Then, Pan et al. [12] promoted this prior to deblur text images. Pan et al. [9] noted that the dark channel of blurred images is less sparse than clear images, then they employed L0 regularization on the dark channel to deblur blurry images. Yan et al. [10] further developed extreme channels prior for better robustness. Instead of exploring the statistical properties of entire images, some methods focus on extracting image edges for kernel estimation [17][18][19][20][21]. This type of method is effective for latent images with solid edges. Another kind of method involves obtaining image patch priors for image estimation, such as a local prior for reducing ringing artifacts [22] and patch priors for modeling image edges [23]. Nevertheless, these methods usually involve heuristic or specially designed optimization algorithms, which means they are somewhat computationally complex.
Lately, the neural network has been employed for image restoration. Sun et al. [24] applied a convolutional neural network (CNN) to remove non-uniform motion blur. Zhang et al. [25] employed recurrent neural networks (RNNs) to model the spatially variant deblurring process. Tao et al. [26] proposed a scale-recurrent network as a pyramid strategy to restoring the sharp image. An end-to-end conditional generative adversarial network (GAN) for motion deblurring was designed and promoted [27,28]. Additionally, a dark and bright channel prior was also embedded into a network for dynamic scene deblurring [29]. However, these methods strongly depend on the training and test data, which limit their applications in the complex and variable turbulent blur.

Comparison of Image Priors
In this section, we analyzed the advantages of the developed image prior compared to previous image priors, also known as image regularization terms.
The L0 regularization [8,31] calculates the number of nonzero values of ∇I. The proposed image regularization term is defined as: where C is a constant that ensures the regularization is loss positive. From the optimization point of view, the developed regularization term is to maximize L1 norm; thus, we named it the maximizing L1 regularization (ML1). To compare the properties of different regularization terms, we simulated atmospheric turbulent blur with a Von Kármán statistic phase screen model [32], as used in [33,34]. To all blurred images, 1% Gaussian noise is added. First, we analyze the effect of image blur on different regularization terms. A natural image, embedded in Figure 1, is used as the original clear image. As shown in Figure 1, the regularization loss of L1 and L2 norms decreases as the image becomes more blurred. Therefore, the L1 and L2 norms are inappropriate image regularization terms because they support blurry images. Levin et al. [11] studied this behavior in detail. The L1/L2 [7] and the L0 [8,31] regularization terms can correctly favor sharp images. However, these priors are non-convex functions, making optimization complicated and unstable. Additionally, compared to the L1 and L2 regularization terms, they reduce the gap between the blurred image and the original image. The ML1 prior behaves correctly for blur and distinguishes the original image and blurred image to a greater extent than the L1/L2 and L0 regularization terms. Furthermore, the ML1 prior is a symmetric function of L1 regularization, which enables it to retain the differentiation and monotonicity of L1 regularization for image blur. Remarkably, ML1 is still a convex function, so the methods to optimize it are relatively straightforward and stable.
where f (I) presents the regular function, and ∇ x and ∇ y are discrete gradient filters. The turbulent blur kernel h ranges from 0 to 50 pixels in size.
Then, we considered the effect of the feature size on the ML1 prior. We simulated three types of classic 2D signals for analysis. As shown in Figure 2, for all types of signals, small-size signals have higher regularization losses under the same blur. This means that the ML1 prior is more sensitive to the ambiguity of small-scale features. To verify the effectiveness of this discovery, we selected an image with rich details and performed similar experiments. In Figure 3, we observe the same phenomenon. The above illustrates the effectiveness of the ML1 prior for preserving small-scale features.  In conclusion, compared to other existing image regularization terms, our proposed regularization term can distinguish between explicit and blurred images and better preserve edge details. Particularly, the presented regularization term is convex, which means many ordinary algorithms are robust to solving this problem. In the next section, the blind image deblurring algorithm using ML1 prior is introduced in detail.

Method
Assuming the image noise is spatially independent Gaussian noise, the general optimization model for the blind deblurring algorithm is: where φ(o) and ϕ(h) denote the regularization terms of the estimated sharp image o and blur kernel h, respectively.

Deblurring Model
Using the ML1 prior for image and the delayed normalization in [15] for the blur kernel, our objective function becomes: where ∇ is the gradient operator. Inspired by the method of non-maximum suppression (NMS) [35,36], we propose a soft suppression strategy for inhibiting artifacts in the blind image deblurring algorithm, which is where γ is a regulator. The value of s ij is regarded as the possibility that an image point o ij is an edge feature point. Applied this function to adjust the regularization weights, the optimization formulation (4) becomes: where represents the Hadamard product. This strategy prevents our algorithm from overprotecting the edge.
The alternating minimization (AM) [3,37] is a commonly used solution in blind image deblurring due to its fast speed. We performed a delayed weighting in the AM algorithm to implement the soft suppression strategy. Precisely, we first calculated s based on the previous iteration image and then updated the image with fixed s. The details of the algorithm are described in the next section.

Optimization
The proposed ML1 prior is convex, so we use the gradient descent method as the specific optimization method for each step. With an initialization of o and h, we conduct one gradient descent step on o and h. The main steps are shown in Algorithm 1.

Algorithm 1 Image deblurring.
Input: The degraded image g, initialization h, and parameters γ and λ o for t = 1 to T do compute s using Equation (9) update o using Equations (10) and (11) update h using Equations (13) and (14) We first compute and fix the s at the tth iteration. With the estimated image obtained from the initialization or the previous iteration, s is updated by: Then, update o at the tth iteration with: for image step ε o > 0. Lastly, we add a non-negative projection on the image, which is to project o by: This step guarantees that all intermediate images conform to the natural rule that the image values are non-negative.

Updating h
In this step, we update h. Given the estimated image, the loss function for kernel h is: Using gradient descent algorithm, h update witĥ for kernel step ε h > 0. Then, project h by Modified from projected alternating minimization (PAM) [15], this iterative algorithm is called the suppressed projected alternating minimization (SPAM).

Pyramid Scheme
In turbulent conditions, the blur kernel is often large. Then, a coarse-to-fine recovery strategy is urgently needed. The pyramid scheme can substantially reduce the computational complexity caused by directly recovering on the blind image size domain and speed up the algorithm [38]. We used a pyramid scheme with a similar method as in [1,7]. We downsampled the blurred image with a size ratio of 2 until the corresponding kernel size is 3 × 3. The blur kernel is initialized as a uniform blur. We start the recovery algorithm at the lowest size and then upsample the recovered sharp image and the estimated blur kernel to initialize the original image and the blur kernel at the next size. At each scale, the number of iterations is not less than 200. The parameters λ o decrease correspondingly with the increase in scale by multiplying by the attenuation factor α. We use bicubic interpolation to complete all the resizing operations. The overall pyramid scheme is as Algorithm 2.

Non-Blind Algorithm
Given the estimated kernel at the finest level, we execute a non-blind image deblurring algorithm similar to the existing approaches [7,9,39]. In this paper, the final clear image is obtained by applying the non-blind image deblurring algorithm in [7].

Experiments and Discussion
We conducted experiments on the turbulent synthetic dataset and real turbulencedegraded images. Our testing environment was a PC running MS Windows 1064 bit version with an Intel Core i5 CPU. We used MATLAB R2014a to test all the MATLAB codes. We compared our method with the previous eight blind image restoration approaches, including the latest [9,15,16,39].

Synthetic Data
To evaluate the performance of our method, we established a turbulence-degraded dataset containing 32 images. The dataset was generated by conducting convolutions between four images and eight simulated atmosphere turbulence blur kernels, plus 1% white noise. The ground truth images were obtained from the dataset in [11].
These eight blur kernels were generated with the Von Kármán statistic phase screen model [32]. The blur kernel of atmospheric turbulence can be modeled as: (15) where h(x, y) is the atmospheric turbulence blur kernel function, A(u, v) is the pupil function of the imaging system, ω(u, v) is the random phase screen function, F −1 presents the inverse Fourier transform, and j is an imaginary unit. The random phase screen function was generated by the inversion of the Von Kármán power spectrum model [40]. Choosing the simulated turbulent images similar to images recorded by a D = 1.50 m telescope through an atmospheric turbulence of r 0 to 0.045-0.055 m, the simulated blur kernel nearly occupied 35 × 35 pixels in the 256 × 256 pixels image pane. Setting the parameters as above, we obtain edmany blur kernels by generating multiple times under the same parameters and setting different parameters. Then, we cropped the edges of the simulated blur kernel image with all pixel values less than 1/225 of its maximum value and resized the kernel image to an odd number by adding zeros to the edges. Finally, we selected these eight blur kernels with less feature repetition. Their sizes range from 29 × 29 to 39 × 39. Figure 4 shows these eight turbulent blur kernels.
In Figures 5 and 6, we depict some visual results recovered by eight competing algorithms and ours on this synthetic turbulent dataset. We downloaded the codes from the authors from their websites and adjusted the parameters to those performing best on our dataset. As for our algorithm, we set γ equal to 3, α equal to 0.9, and final λ o equal to 0.0001. As shown in Figures 5 and 6, our method performed better than the other algorithms visually, especially in the details of the restored images. In Figure 5, our method recovers clearer fences, fine-grained patterns on clothes, and thin branches of bushes. In Figure 6, our restored images contain more obvious textures and hair. In general, the image restored by our method has sharper edges compared to [1,2,9,15,16,39] and less artifacts compared to [2,7,22].   [16], (e) Pan et al. [9], and (f) our methods.  [22] (c) Krishnan et al. [7], (d) Perrone and Favaro [15], (e) Wen et al. [39], and (f) our methods.
Additionally, we evaluated the restored results using the sum of squared differences (SSD) [11], error ratio (ER) [11], peak signal-to-noise ratio (PSNR) [41], and structural similarity (SSIM) [42]. In this study, the sparse prior algorithm in [4] was used as the non-blind image deblurring algorithm for calculating the error ratios. Figure 7 shows the cumulative performance of SSD error ratios on the synthetic turbulent dataset. Our curve reaches the highest cumulative distribution of all ERs. Figure 8 presents the cumulative histogram of SSD results per image.   In Figure 9, we plot the gradient distributions of an image restored by different methods. The two best-performing methods were selected for comparison. As shown in Figure 9, our curve fits the original curve best. This confirms the effectiveness of our method.

Real Data
In this part, we compare our method with five competing methods on real turbulent images. We blindly deblurred an image with different scenes, and our method restored more explicit structures and less artifacts than others. One part of the real images is blurred images that we actually selected, with scenes involving houses, trees and moving targets. The other part of the real images is real blurred images from the sky, with scenes involving planets, moons, and nebulae. The deblurring results of some images are shown here. In Figure 10, compared to [15,16], our image has fewer artifacts in the house part and, compared to [9,39], more details in the bush part. In Figure 11, our method recovers the detailed structures with less artifacts, which can be observed from the circled part. In Figure 12, although the image noise is more serious, our method is also competitive compared to the others. In Figure 13, our image has more natural edges than other methods, which can be clearly observed from the enlarged parts. Notably, the shape of the blur kernel estimated by our method is also more in line with the turbulent kernel, which has multiple peaks.   Entropy was also used for quantitative evaluation, which is calculated as: where p v denotes the probability of gray level v in the image. Table 2 shows the entropy values of the deblurred images. The entropy value is positively associated with the clarity of images. The entropy values obtained by our proposed method are commonly higher than those of previous studies [7,9,15,16,39]. The proposed method achieves the highest average entropy value, which verifies its stability.

Limitation
In our approach, we model image noise with the Gaussian function. If the blurred image contains other scattered noises, such as speckle noise, the model may misjudge the noise as the texture of the image. In this case, the SPAM algorithm may recover a noisy image. However, to the best of our knowledge, this is also challenging for other blind image deblurring algorithms. How to balance the removal of outliers and the preservation of image details in the blind deblurring algorithm is still a topic requiring future research.

Conclusions
In this paper, we proposed a maximizing L1 regularization prior, which is motivated by effectively supporting sharp latent images instead of turbulence blurry images. To efficiently restore the image with this prior, we developed a suppression projected alternating minimization algorithm to estimate latent sharp images and blur kernels. Then, we simulated an atmospheric turbulent blur dataset for quantitative evaluation. Experiments on this dataset and other turbulent images indicated that the output images deblurred by the proposed method have superior visual performance. The proposed approach can be extended to other types of degraded images. In future work, we plan to develop our model to handle more severe noise.