Effective Alternating Direction Optimization Methods for Sparsity-Constrained Blind Image Deblurring

Single-image blind deblurring for imaging sensors in the Internet of Things (IoT) is a challenging ill-conditioned inverse problem, which requires regularization techniques to stabilize the image restoration process. The purpose is to recover the underlying blur kernel and latent sharp image from only one blurred image. Under many degraded imaging conditions, the blur kernel could be considered not only spatially sparse, but also piecewise smooth with the support of a continuous curve. By taking advantage of the hybrid sparse properties of the blur kernel, a hybrid regularization method is proposed in this paper to robustly and accurately estimate the blur kernel. The effectiveness of the proposed blur kernel estimation method is enhanced by incorporating both the L1-norm of kernel intensity and the squared L2-norm of the intensity derivative. Once the accurate estimation of the blur kernel is obtained, the original blind deblurring can be simplified to the direct deconvolution of blurred images. To guarantee robust non-blind deconvolution, a variational image restoration model is presented based on the L1-norm data-fidelity term and the total generalized variation (TGV) regularizer of second-order. All non-smooth optimization problems related to blur kernel estimation and non-blind deconvolution are effectively handled by using the alternating direction method of multipliers (ADMM)-based numerical methods. Comprehensive experiments on both synthetic and realistic datasets have been implemented to compare the proposed method with several state-of-the-art methods. The experimental comparisons have illustrated the satisfactory imaging performance of the proposed method in terms of quantitative and qualitative evaluations.


Background and Related Work
Single-image blind deblurring for imaging sensors has recently received increasing attention in modern imaging applications, e.g., the Internet of Things (IoT), astronomical imaging, biomedical imaging, computational photography and microscopy [1][2][3]. It is well known that the image pixel intensity can be determined by the total incoming light sensed by the imaging sensor over the exposure time. As shown in Figure 1, the discrete image degradation model can be written as follows: where B (x, y) is the observed image after camera exposure, L (x, y) denotes the latent sharp image to be restored, ξ (x, y) denotes the additive white Gaussian noise, w m is a weight, which essentially represents the length of exposure time at camera pose m, and H m is a transformation matrix related to the camera rotation or translation at pose m during exposure. In this work, we only consider the case of uniform (i.e., spatially invariant) image deblurring. Thus, the matrix H m only corresponds to the camera translation along both the X and Y axes. For the sake of simplicity, the original image degradation model (1) can be rewritten as a convolution version as follows: where ⊗ is the mathematical operation of convolution, k denotes the blur kernel related to the weight w and transformation matrix H in (1). The purpose of single-image blind deblurring is to recover both k and L from only one blurred image B. It is a challenging ill-conditioned inverse problem, since many different pairs k and L can lead to the same B [4]. The constraints on both the blur kernel and the latent sharp image should be exploited to select the optimal pairs k and L for enhancing imaging performance. To cope with the ill-conditioned nature of blind deblurring, many statistical priors learned from blur kernels and latent sharp images have been developed to regularize the restoration process. In the literature [4,5], current single-image blind deblurring methods are widely divided into two categories: (1) methods that simultaneously estimate both the blur kernel and latent sharp image; (2) methods that first estimate the blur kernel, then recover the latent sharp image. It is well known that the support size of the blur kernel is often extremely smaller compared to the image size. Therefore, the joint maximum a posteriori (MAP) estimation of k and L often fails since the number of unknowns is larger than the number of known variables in B. In contrast, the estimation of the blur kernel can be obtained accurately through the MAP estimation of k alone [4]. To guarantee the high-quality blind deblurring, this paper mainly focuses on the second type of method, i.e., estimating the blur kernel first and then dealing with the corresponding non-blind deconvolution problem.
The pioneering work [6] mainly focused on the estimation of simple and small blur kernels, which are very rare in many practical scenarios. To make blind deblurring more practical, most current state-of-the-art methods were usually proposed by exploiting the prior knowledge from the statistics of blur kernels and sharp images. In 2006, Fergus et al. [7] contributed the original work on practical blind deblurring where the blur kernels were quite large and complex. In particular, the authors proposed a variational Bayesian image deblurring model by combining the mixture-of-Gaussian image prior with the mixture-of-exponential kernel prior. In [8], the blur kernel prior was assumed to follow an exponential distribution. Under the MAP framework, the exponential distribution results in an L 1 -norm constraint on kernel intensity, which has a good interpretation on the sparsity of the blur kernel [9,10]. Under some imaging conditions, the blur kernel can also be assumed as a (piecewise) sufficiently smooth function. As a consequence, many researchers proposed to replace the L 1 -norm of kernel intensity with its squared L 2 -norm version [11][12][13]. Current experiments have shown that both L 1 -and squared L 2 -regularized optimization methods could achieve accurate kernel estimation on the benchmark dataset introduced by [4]. In the case of large blur kernels, however, it is difficult to robustly and accurately estimate the blur kernels using these methods mentioned. Many efforts [4,5] have been devoted to theoretically explain the reason why it is difficult for accurately estimating the blur kernels, especially for the large ones in practice. More recently, hybrid sparsity priors on blur kernels [14,15] have been considered and achieved robust image restoration results.
Once the blur kernel is estimated, the blind deblurring problem (2) essentially becomes a non-blind image deconvolution. During the past several decades, numerous numerical methods have been developed to handle non-blind deconvolution. One of the most popular methods is the Tikhonov regularization [16,17], followed by its various extensions [18,19]. These methods can be easily implemented, but commonly generate over-smoothing effects on the restored images. Other widely-used methods, such as the Richardson-Lucy method [20] and Wiener filter [21], easily suffer from noise amplification and ringing-like artifacts. To overcome the undesirable artifacts, Yuan et al. [22] proposed to develop a progressive inter-scale and intra-scale image deconvolution approach based on the bilateral Richardson-Lucy method. Current research illustrates that images have the properties of sparse gradients. Many efforts [23][24][25] were made to enhance non-blind deconvolution by imposing the total variation (TV) regularizer. From a statistical point of view, the TV regularizer corresponds to an assumption of a Laplacian sparse prior on image gradients. Recently, the extended TV regularizers, such as non-convex TV (NCTV) [11,13,26] and higher-order TV (HOTV) [27,28], have been attracting increasing attention for improving non-blind deconvolution. Both TV and HOTV regularizers have also been combined to overcome the potential disadvantages existing in these two regularizers [29][30][31]. The newly-developed total generalized variation (TGV) regularizer, originally proposed by Bredies et al. [32] in 2010, achieved great success for the restoration of blurred images [33][34][35]. Motivated by the concepts of non-local means (NLM) and graph Laplacian, the non-local TV (NLTV) regularizer has significantly improved the deconvolution quality [36][37][38][39]. The NLTV-regularized variational models can guarantee the highest-quality deconvolution because they take full advantage of the high degree of geometrical self-similarity that is inherent in natural images.

Motivation and Contributions
In the current literature [7][8][9][10][11][12][13], most existing blur kernel estimation methods were proposed based on the assumptions that the blur kernel was spatially sparse or piecewise smooth within the support of a continuous curve. As a consequence, the proposed methods could not always guarantee high-accuracy estimation under certain degradation conditions. In recent years, more attention has been paid to the sparse image priors for improving the estimation accuracy. To further enhance the estimation quality, in our opinion, it is still necessary to investigate advanced sparsity constraints on the blur kernel. The robust estimation method will be proposed in this paper by taking into account the sparsity and smoothing properties of the blur kernel. In particular, the sparsity property is promoted using the L 1 -norm of kernel intensity; the smoothing property is utilized through the introduction of the squared L 2 -norm of the intensity derivative. By making full use of the advantages of both the L 1 -norm and the squared L 2 -norm on kernel prior representation, the proposed method could potentially generate satisfactory estimation under more different degradation conditions. Essentially, most of the previous works [8,9,[11][12][13] on blur kernel estimation can be considered as a special case of our proposed hybrid regularization method. If we only use the L 1 -norm term, it can take full advantage of the property of spatial sparsity, but the resulting estimated blur kernel easily suffers from the isolated points [15]. If we only use the squared L 2 -norm term, the continuous and smoothing properties of blur kernels under certain imaging conditions could be well preserved. However, the potential spatial sparsity property may be ignored, leading to inaccurate estimation of the blur kernel. Therefore, to guarantee the accuracy of the estimated blur kernel, it is necessary to combine the L 1 -norm of kernel intensity with the squared L 2 -norm of the intensity derivative in our work. If images are degraded by Gaussian, average or pillbox (disc) blur kernels, which have the weak spatial sparsity properties, but high smoothing properties, the proposed hybrid regularization method could theoretically generate higher estimation accuracy compared with traditional single regularization methods. It is worth mentioning that the hybrid blur kernel prior proposed in this work is extremely different from the current hybrid versions [14,15]. Current work [13,[40][41][42] has illustrated that the L 0 quasi-norm has a good natural interpretation of the sparsity property of the image gradient and benefits for image detail enhancement. In particular, it performs well in penalizing small gradient magnitudes and encouraging large ones to preserve fine details. To improve the accuracy of blur kernel estimation, the L 0 quasi-norm of the image gradient is also incorporated into our blur kernel estimation method. Owing to the non-convex nature of the L 0 quasi-norm and the non-smooth nature of the L 1 -norm, the commonly-used numerical methods could not be effectively adopted to solve the blur kernel estimation problem. To guarantee a feasible solution, the resulting non-convex non-smooth optimization problem will be effectively dealt with by developing an alternating direction method of multipliers (ADMM)-based numerical method [43]. The preliminary results on blur kernel estimation can be found in our previous short-version conference paper [44].
Existing work has illustrated that the TV regularizer, first proposed by Rudin et al. [23] in 1992, has the capacity of preserving edges and smoothing flat regions. TV-regularized variational image restoration models with the L 1 -norm [11] or the squared L 2 -norm [24] data-fidelity terms have gained considerable attention. However, the image quality could be degraded because the results often suffer from undesirable staircase-like artifacts in regions with gradual intensity variations [45]. The reason behind this phenomenon is that the TV regularizer favors solutions that are piecewise constant. To effectively suppress the artifacts, many extensions of TV [11,13,[26][27][28][36][37][38][39] could be used to improve the image quality. For example, the patch-based NLTV regularizer has the capacity of guaranteeing the highest-quality image restoration. However, the NLTV-regularized variational model is practically limited due to the high computational cost. To make it easier to implement blind deblurring in practice, it is necessary to balance the trade-off between computational cost and imaging performance. Motivated by the success of the TGV regularizer, we tend to propose an effective non-blind deconvolution method based on the TGV regularizer of second-order (i.e., TGV 2 ) [32]. The TGV 2 regularizer is able to suppress the undesirable artifacts, while preserving the image edges since it favors piecewise polynomial intensities [34]. The quality of restored images could be correspondingly enhanced. From an optimization point of view, the resulting image deconvolution model could not be directly solved using traditional numerical methods because of the non-smooth nature of the TGV 2 regularizer. To achieve a robust and effective solution, an ADMM-based optimization method will be developed to solve the resulting non-smooth minimization problem. In particular, the original complex minimization problem can be decomposed into several simple subproblems by introducing some auxiliary variables. Each of these subproblems has a closed-form solution or can be efficiently solved using the current numerical method. The effectiveness of the proposed method will be demonstrated using comprehensive experiments on both synthetic and realistic blurred images.
In conclusion, the main contributions of this paper, given the state-of-the-art research work, are mainly summarized by the following three aspects: • To accurately estimate the blur kernel, a hybrid regularization method was proposed by combining the L 1 -norm of kernel intensity with the squared L 2 -norm of the intensity derivative. An alternating direction method was presented to effectively solve the resulting blur kernel estimation problem. • The TGV 2 -regularized variational model with an L 1 -norm data-fidelity term was proposed for enhancing the non-blind deconvolution result. To guarantee the stability and effectiveness of the solution, an ADMM-based numerical method was developed to solve the resulting non-smooth optimization problem. • The satisfactory blind deblurring performance of the proposed method has been illustrated using comprehensive experiments on both synthetic and realistic blurred images (with large blur kernels). The proposed method has also been successfully exploited for single-image deblurring in the field of ocean engineering.
The main benefit of the proposed method is that it takes full advantage of the hybrid constraints for blur kernel estimation and the TGV 2 regularizer for non-blind deconvolution. Therefore, it can accurately estimate the blur kernel and guarantee high-quality image deconvolution. Experiments using synthetic, as well as realistic blurred images will be implemented to verify the effectiveness of our proposed method in practical applications.

Hybrid Regularized Blur Kernel Estimation
As discussed in Section 1.2, our robust two-step framework for single-image blind deblurring is illustrated in Figure 2. This section mainly focuses on the first blur kernel estimation step, which is separated into the following two aspects: (1) sharp edges restoration; and (2) blur kernel estimation. In order to enhance the deblurring performance, we exploit the following statistical priors for blur kernel estimation: an L 0 -sparsity prior on the latent gradient image x and a hybrid sparsity prior on the blur kernel k. Under these sparsity-constrained priors, blur kernel estimation in this paper is equivalent to solving the following minimization problem: where γ, η 1 , η 2 are predefined positive regularization parameters, the L 0 quasi-norm • 0 counts the number of nonzero elements, x denotes ∇L = (∂ h L, ∂ v L) T and y denotes ∇B = (∂ h B, ∂ v B) T with ∂ h and ∂ v being the finite differences along the horizontal and vertical directions, respectively. The proposed blur kernel estimation model (3) is mainly composed of four terms: The first term, called the squared L 2 -norm data-fidelity term, denotes a measure of the distance between the restored data and the observed version. The second term is the L 0 quasi-norm regularization term, which can preserve the sparsity of natural image gradients. The third and fourth terms are, respectively, the L 1 -norm and the squared L 2 -norm constraints on the blur kernel, which can stabilize the final estimation result. The accurate estimation of the blur kernel is beneficial for generating high-quality non-blind image deconvolution.

Sharp Edge Restoration
Since x and k are independent in (3), in the first step, the latent sharp edges x at the (m + 1)-th outer iteration can be recovered by solving the following minimization problem: for m = 0, 1, · · · , M max with M max denoting the maximum number of outer iterations. As discussed in [11,13], strong edges are not always beneficial for accurate estimation of the blur kernel. To select the informative edges, an effective method used to measure the usefulness of gradients is given by: where B is the blurred image and N h (p) denotes an h × h window centered at pixel p ∈ Ω (image domain). The measure metric (5), first proposed by Xu and Jia [11], enables accurate estimation of the blur kernel by removing some narrow strips. To incorporate the measure metric into our kernel estimation framework, the problem (4) can be reformulated as: where • represents the pointwise product and κ (p) = exp(− |r (p)| 0.8 ) for p ∈ Ω with the variable r (p) being defined in (5). For the sake of simplicity, the convolution version in (6) is expressed in a matrix-vector multiplication form. In (6), K m is a block Toeplitz matrix with Toeplitz blocks transformed from the blur kernel k at the m-th outer iteration. x and y represent the vector version of x and y, respectively. It is well known that the model (6) is difficult to solve directly because of the non-smooth and non-convex natures of the L 0 quasi-norm κ • x 0 . To guarantee solution efficiency and stability, we propose to develop an ADMM-based numerical method [43,46,47] to solve the unconstrained optimization problem (6). To apply ADMM, we first replace x by v and then transform (6) into the following constrained optimization problem: Note that the updates of v and x are independent of each other. Let L A (v, x; ϕ x ) represent the augmented Lagrangian function of (7), which is defined as follows: where β 1 is a pre-defined penalty parameter and ϕ x denotes the Lagrangian multiplier. In particular, ADMM solves problem (6) by minimizing L A (v, x; ϕ x ) with respect to v and x alternatively given the other fixed, followed by an update of the Lagrangian multiplier ϕ x , i.e., ) for i = 0, 1, · · · , I max with I max denoting the maximum number of inner iterations. Here, τ ∈ 0, ( √ 5 + 1)/2 denotes the step length. In order to accelerate the convergence of numerical solution, during each iteration, the penalty parameter β 1 can be updated as β 1 ← ρβ 1 with ρ being a positive step length. The explicit solution of the v-subproblem in (9) can be directly obtained using the element-wise hard thresholding operator formulated in [48] where H a,b (·) is defined as: where both a and b are intermediate variables.
Essentially, the x-subproblem in (9) is a least-squares optimization problem; the corresponding normal equation can be readily obtained as follows: where superscript denotes the transpose operator for real matrices or vectors and I denotes an identity matrix. Under the periodic boundary condition, K m K m is a block circulant matrix with circulant blocks. It can be diagonalized using the two-dimensional discrete Fourier transform. Let F denote the forward fast Fourier transform (FFT) operator. Applying by F on both sides of (11) and yields: To decrease the computational cost, both F (K m )F (K m ) + β 1 F (I) and F (K m )F (y) can be computed only once at the beginning of the iterative algorithm. Thus, solving (12) is straightforward, which means that it is relatively easy to achieve the solution of (11). As a consequence, the solution of the least-squares optimization problem (11) is given by: where F −1 (·) denotes the inverse FFT operator and F (·) represents the complex conjugate operator. The minimization process (9) is implemented alternately until the solution converges to the optimal one. Finally, the recovered sharp edge x m+1 = x m,I max is achieved to enhance the accuracy of blur kernel estimation in the next step. The whole optimization procedure of ADMM for Subproblem (6) is summarized in Algorithm 1.

Blur Kernel Estimation
In the blur kernel estimation step, given the recovered sharp edge x m+1 , the blur kernel k in (3) at the (m + 1)-th outer iteration can be estimated by considering the following minimization problem: where the optimal parameters η 1 and η 2 are manually selected by extensive experiments. It is obvious that Model (14) is essentially a convex optimization problem. Analogous to the optimization of problem (6), ADMM can also be adopted to efficiently solve (14) in our experiments. We first replace k by h and then obtain the corresponding augmented Lagrangian function L A h, k; ϕ k as follows: where β 2 is a pre-defined penalty parameter and ϕ k denotes the Lagrangian multiplier. Given the fixed k m,j , the minimization of L A h, k; ϕ k with respect to h could be easily handled through the widely-used shrinkage operator [49][50][51], which operates pointwise on scalars or matrices. The solution h j+1 at the (j + 1)-th inner iteration is given by: with k m,0 = k m . Here, the sign function sign (·) is defined as: Given the fixed h j+1 , the minimization of L A h, k; ϕ k with respect to k is equivalent to solving a least-squares optimization problem. Analogous to solving the x-subproblem in (9), the corresponding solution k m,j+1 is obtained as follows: It is tractable to obtain the efficient solution k m,j+1 in (17) using one forward and one inverse FFT. At each iteration, the Lagrangian multiplier ϕ k can be updated as follows: for j = 0, 1, · · · , J max with J max denoting the maximum number of inner iterations. The estimated blur kernel k m+1 = k m,J max can be achieved for sharp edge restoration in the next step. Note that the convergence of our proposed numerical algorithm can be guaranteed according to the existing convergence results for ADMM in the literature [43,52,53]. Finally, our proposed hybrid regularized variational model for blur kernel estimation is summarized in Algorithm 2.

Robust Non-Blind Deconvolution
This section mainly focuses on developing a high-order variational model for robust non-blind deconvolution. For the sake of better writing, the original image degradation model B = L ⊗ k + ξ in (2) can be rewritten as follows: Once the blur kernel K (i.e., k in Algorithm 2) is estimated accurately, blind deblurring can be simplified to the non-blind deconvolution problem. As a consequence, this problem could be handled through the commonly-used regularization methods. One of the most famous methods is the TV-regularized deconvolution method [23][24][25]. However, the undesirable staircase-like artifacts generated in restored images often lead to significant degradation of visual image quality. In order to enhance the imaging performance, a robust non-blind image deconvolution method will be proposed based on the second-order regularizer TGV 2 [32]. Recently, TGV 2 has been successfully utilized as a regularization scheme in various practical applications [33,[54][55][56] and outperforms the popular TV regularizer. In particular, TGV 2 is capable of preserving image edges and suppressing undesirable artifacts.
Inspired by the work in [11], the L 1 -norm data-fidelity term will be incorporated into our TGV 2 -regularized non-blind deconvolution to suppress the potential outliers. The used assumption behind the squared L 2 -norm data-fidelity term is that the data-fidelity costs follow a Gaussian distribution. This assumption often fails because the squared L 2 -norm could make the restored images vulnerable to undesirable outliers. In contrast, the L 1 -norm introduced in this work is much more robust to the presence of outliers compared with the squared L 2 -norm. The proposed non-blind deconvolution model L 1 -TGV 2 is given by: where λ > 0 denotes a predefined regularization parameter. For a scalar field L ∈ L 1 (Ω), the discretized TGV 2 α (L) [32] is defined as follows: where α 1 and α 0 are positive tuning parameters, C 2 c Ω, R 2 denotes the space of the vector field and E (V) = 1 2 ∇V + ∇V T with V = [V 1 V 2 ] T being the symmetrized gradient of a complex-valued vector field V. Due to the non-smooth nature of the TGV 2 regularizer, in this paper, we propose to develop an ADMM-based numerical algorithm to effectively solve the non-smooth optimization problem (20). Three auxiliary variables W, Y and Z are first introduced, and (20) is then transformed into the following constrained minimization problem: min W,Y,Z,L,V It is obvious that L and V are coupled together. For the fixed values of L and V, the updates of W, Y and Z are independent of each other. Thus, the variables W, Y, Z, L and V can be decomposed into two blocks, i.e., (W, Y, Z) and (L, V). Let L A (W, Y, Z, L, V; ξ, ζ, η) denote the augmented Lagrangian function of (22), which can be defined as follows: where ξ ∈ R mn , ζ ∈ R 2mn and η ∈ R 4mn denote the Lagrange multipliers and ρ 1 , ρ 2 and ρ 3 represent the positive penalty parameters that control the weights of penalty terms. It is numerically intractable to directly obtain the solutions of (22) through commonly-used methods [57]. In order to guarantee a stable solution, it is necessary to alternatively solve the W-, Y-, Z-, L-and V-subproblems and then update the Lagrange multipliers (i.e., ξ, ζ and η) until the obtained solution meets the predefined threshold. In particular, each of these subproblems has a closed-form solution or can be efficiently solved using the existing simple numerical method.

{W, Y, Z}-Subproblems
Given the fixed values of variables L t , V t , ξ t , ζ t and η t , the W-subproblem, the Y-subproblem and the Z-subproblem in Equation (23) can be efficiently solved by considering the following L 1 -regularized least-squares minimization problems: Note that the unknown variables W, Y and Z are componentwise separable in the {W, Y, Z}-subproblems (24)- (26). These subproblem can be effectively dealt with through the commonly-used shrinkage operator [49,50]. This operator is fast and easy to implement in practice. The solutions W t+1 , Y t+1 and Z t+1 are obtained as follows:

(L, V)-Subproblems
The minimization with respect to (L, V) in (23) is essentially a least-squares optimization problem. However, it is impossible to directly obtain the solutions L t+1 and V t+1 through the forward and inverse FFT operators because the updates of L and V are coupled to each other. To guarantee the solution stability, the minimizations with respect to both L t+1 and V t+1 should be simultaneously implemented. Given the fixed values of W t+1 , Y t+1 , Z t+1 , ξ t , ζ t and η t , the coupled (L, V)-subproblem in the augmented Lagrangian function (23) is quadratic, resulting in the following system of linear equations: Instead of directly solving the system of linear equations (30), we tend to solve the corresponding first-order necessary optimality conditions as follows: For the sake of better reading, the original system of linear equations (31) can be rewritten as follows: with: and: Let F denote the discrete Fourier transform operator for real (complex) matrices or vectors. To efficiently solve the system of linear equations (32), we multiply both sides of (32) by F , such that the coefficient matrix will be blockwise diagonal, i.e., Essentially, (33) is a linear system with three equations and three variables. We propose to directly use Cramer's rule to effectively yield the closed-form solutions as follows: where F −1 (•) represents the inverse Fourier transform operator. In particular, we have the determinants det L = S T 2 T 3 * , det V 1 = T 1 S T 3 * , det V 2 = T 1 T 2 S * and det T = . The definition of determinant | · | * in this work is briefly introduced as follows: r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33 * To shorten the computational time, all matrix elements in det T for (34) should be calculated before the execution of our ADMM-based numerical algorithm. At each iteration, D 1 , D 2 and D 3 are first can be easily calculated correspondingly. The final solutions L t+1 , V t+1 1 and V t+1 2 can be naturally obtained using Cramer's rule (34).

Update the Lagrange Multipliers
At each iteration of our proposed numerical method, the Lagrange multipliers (ξ, ζ, η) should be updated as follows: where the step length τ = 1.618 is used throughout this paper. In conclusion, an ADMM-based numerical method was proposed to decompose the original complex optimization problem (20) into several simpler subproblems. Each subproblem has a closed-form solution or can be efficiently solved using the existing numerical method. In particular, the {W, Y, Z}-subproblems (24)(25)(26) could be easily solved using the shrinkage operator. The solutions of L and V were simultaneously obtained through Cramer's rule (34). The optimization procedure of our proposed method for non-blind image deconvolution is summarized in Algorithm 3. Compute W t+1 according to W t+1 = shrinkage KL t − B + ξ t /ρ 1 , 1/ρ 1 .

Experimental Results and Discussion
Comprehensive blind deblurring experiments on both synthetic and realistic blurred images will be performed to verify the effectiveness of our proposed method in this section.

Experimental Settings
The proposed blind deconvolution framework was evaluated on a synthetic blur-image dataset [4] and realistic blurred images. The synthetic dataset has been widely exploited as a benchmark dataset to evaluate the performance of blur kernel estimation. Our numerical experiments were implemented using MATLAB R2011a (The MathWorks, Natick, Inc., MA) on a machine with a 3.30-GHz Intel(R) Pentium(R) G3260 CPU and 4 GB RAM. For both synthetic and realistic datasets, the parameter values for blur kernel estimation in Section 2 were set as follows: τ = 1.618, ρ = 3, I max = 5, γ = 5 × 10 −2 , η 1 = η 2 = 10 −3 and M max = 15. The resulting optimal parameters for non-blind deconvolution in Section 3 were set empirically, i.e., ρ 1 = 50, ρ 2 = 0.5, ρ 3 = 5, α 1 = 1, α 2 = 1.5 and T max = 10. The deblurring results have illustrated the satisfactory performance of the manually-selected parameters in our experiments. For the sake of better comparison, the competing blind deblurring methods yield the restoration results with the input parameters manually optimized by the authors. To further improve deblurring performance, there is a great potential to develop automatic estimation methods to adaptively select the optimal parameters in our future work. Similar to [4], the sum of squared differences (SSD) and SSD ratio were used simultaneously to quantitatively evaluate the performance of blur kernel estimation in Section 4.2. In particular, the SSD ratio is measured between the deconvolution error with the estimated blur kernel and the deconvolution error with the ground truth kernel [4].

Experiments on Synthetically-Blurred Images
Numerous experiments are implemented in this subsection to evaluate the performance of our proposed method on one widely-used synthetic blur image dataset [4], which can be downloaded from the link: www.wisdom.weizmann.ac.il/~levina/papers/LevinEtalCVPR2011Code.zip. In Figure 3, it could be found that the dataset is composed of four grayscale images of the size 256 × 256 and eight different uniform blur kernels, resulting in a total of 32 synthetic blurred images in our experiments. The proposed method will be compared with several state-of-the-art blind deblurring methods [7,[11][12][13]58] in terms of SSD and the SSD ratio. In order to guarantee an unbiased comparison, the final deblurred results are all generated using the sparse non-blind deconvolution method proposed in [58]. Furthermore, to enhance the robustness of blur kernel estimation, a widely-used multiscale scheme [7] was introduced in Algorithm 2. Experimental results on the SSD ratio for different kernel estimation methods are summarized in Figure 4. It can be found that our proposed method is able to generate more robust estimation results on this synthetic dataset under consideration in most of the cases. In contrast, the accuracy of blur kernel estimation for other competing methods is limited due to the simple assumption of the blur kernel prior.  The objective function Equation (3), comprising the L 0 quasi-norm and the L 1 -norm regularization terms, is both non-convex and non-smooth. Thus, it is much more difficult to estimate the complexity of the proposed method from a theoretical point of view. It is well known that computational cost is highly dependent on the algorithm complexity. For the sake of simplicity, we only tend to compare the computational time for different blur kernel estimation methods under the same imaging conditions. The methods of Xu and Jia [11] and Cho and Lee [12] are efficiently implemented using C code in our experiments. In contrast, the other competing methods are performed in MATLAB. Since all test images in Figure 3 are mainly composed of fine details of different textures, we take images Im02 and Im04 as examples to evaluate the computational efficiency. The computational time of different competing blur kernel estimation methods is summarized in Table 1. It could be found that the methods of Xu and Jia [11] and Cho and Lee [12] generate the lowest computational cost due to the C implementation. Our proposed method yields significantly faster computational speeds compared with Fergus et al. [7] and Levin et al. [58]. However, the method of Pan and Su [13] achieves the highest computational efficiency due to the fast alternating direction optimization method. As shown in Table 1, the proposed method yields the best evaluation results in terms of the SSD metric under consideration in most of the cases. The satisfactory performance of our proposed method benefits from the combination of the L 1 -norm of kernel intensity with the squared L 2 -norm of the intensity derivative. It could be further observed that both computational cost and image quality highly depend on the size of the blur kernel. As the size becomes larger, the computational cost obviously becomes higher, and the image quality could becomes lower for all competing methods under different imaging conditions. In addition, there is no significant difference in blind deblurring performance between different test images for the same blurring degradation. For a better comparison, Figure 5 visually displays a test image from the synthetic blur-image dataset [4] recovered by different deblurring methods. As can be observed, our proposed method estimates a more accurate blur kernel. The proposed method generates the highest quality of deblurred image, since it achieves a more natural-looking appearance. We can conclude that our proposed method has superior performance compared with other competing methods.  [4]. From top-left to bottom-right: (a) latent sharp image, deblurred versions with (b) the truth blur kernel, estimated blur kernels generated by (c) Fergus et al. [7], (d) Xu and Jia [11], (e) Cho and Lee [12], (f) Levin et al. [58], (g) Pan and Su [13] and (h) our proposed method, respectively. The estimated blur kernels are illustrated in the upper-left panels.

Experiments on a Large Blur Kernel
Single-image blind deconvolution with a large blur kernel is an extremely challenging problem in practical application. It is necessary to investigate whether our proposed method is able to deal with the large blur kernel. In this subsection, our proposed method will be compared with three state-of-the-art deblurring methods, i.e., Xu and Jia [11] in ECCV-2010, Krishnan et al. [9] in CVPR-2011 and Pan et al. [59] in CVPR-2016. Once the blur kernel is estimated, the method of Xu and Jia [11] implements a robust non-blind deconvolution by integrating the L 1 -norm data-fidelity term and the TV regularizer. Krishnan et al. [9] directly uses the hyper-Laplacian prior-based fast non-blind deconvolution method [26] to generate the final recovered image. The non-blind deconvolution results in our proposed method will be achieved using the combination of the L 1 -norm data-fidelity term and the TGV 2 regularizer summarized in Algorithm 3. In contrast, the method in Pan et al. [59] simultaneously estimates the latent sharp image and blur kernel by introducing the assumption of the dark channel prior.
The blind deblurring results with large blur kernels are visually displayed in Figures 6 and 7. As shown in Figure 6, the method of Krishnan et al. [9] is unable to yield the satisfactory estimation of the large blur kernel of the size 159 × 159. The resulting deblurred image suffers from the loss of important geometrical structures, which significantly degrades the visual quality. In contrast, the methods of Xu and Jia [11], Pan et al. [59] and our proposed method have the capacity of guaranteeing the accurate estimation of the blur kernel. The main geometrical structures in the deblurred images could be reconstructed correspondingly. However, an excess smoothing could be observed in the methods of Xu and Jia [11] and Pan et al. [59], which leads to the loss of small structural features. Owing to the second-order regularizer TGV 2 , our non-blind deconvolution method produces a much more natural-looking result than Xu and Jia [11] and Pan et al. [59]. As shown by the arrows in Figure 6, more details could be preserved in our proposed method; thus, the resulting deblurring performance outperforms other comparative methods. More blind deblurring results with large blur kernels are visually illustrated in Figure 7. The sizes of the estimated blur kernels in our experiments are 101 × 101, 95 × 95 and 91 × 91, respectively. It could be found that the method of Krishnan et al. [9] fails to yield the accurate estimation of the blur kernel and generates unsatisfactory deblurring performance. In contrast, the proposed method is able to generate comparable results to the state-of-the-art blind deblurring methods, i.e., Xu and Jia [11] and Pan et al. [59]. The final high-quality recovered images could be achieved with more structures and details preserved. Therefore, there is a huge potential to use our proposed method to restore the blurred images with large blur kernels in practice.   [11], (c) Krishnan et al. [9], (d) Pan et al. [59] and (e) our proposed method, respectively. The sizes of the estimated blur kernels from top to bottom are 101 × 101, 95 × 95 and 91 × 91, respectively.

Experiments on Ocean Engineering
In the field of ocean engineering, computer vision-assisted automatic detection and tracking systems with airborne and shipborne imaging sensors have been widely used to improve maritime control, safety and rescue operations. However, the resulting imaging performance sometimes suffers from motion blur, noise, haze and sensor nonlinearities, which could significantly degrade the visual image quality under poor weather conditions. In this paper, we mainly focus our attention on the restoration of blurred images, since this degradation condition is more common than other conditions in ocean engineering. The experimental images captured with airborne and shipborne imaging sensor systems, as well as the corresponding blur kernel estimation and image deconvolution results are visually displayed in Figures 8 and 9. The sizes of the estimated blur kernels are 35 × 35 and 95 × 95, respectively. As shown in Figure 8, the method of Krishnan et al. [9] is unable to achieve high-quality blur kernel estimation resulting in unsatisfactory deblurring performance with ringing-like artifacts. In contrast, Pan et al. [59] and our proposed method have the capacity of producing accurate estimation of the blur kernel for this example. The final high-quality restored images could be guaranteed using non-blind deconvolution methods. However, due to the low-contrast structure shown in Figure 9, the latest dark channel prior-based method [59] could not accurately estimate the blur kernel in the case of shipborne imaging. The reason behind this phenomenon may be that the statistical properties between images captured by shipborne cameras and natural images are essentially different. The assumption of the dark channel prior is not always valid under different imaging conditions. The method of Xu and Jia [11] also fails to estimate the blur kernel and generates a low-quality restored image. Figure 9 visually illustrates that our proposed method is still able to guarantee accurate kernel estimation and image deconvolution. More geometric structures and fine details could be preserved in our recovered images, beneficial for detecting and tracking moving vessels in practice. The maritime control, safety and rescue operations could be correspondingly improved in the field of maritime management and ocean engineering.

Experiments on More Realistic Blurred Images
To better evaluate our proposed method, this subsection is concluded by testing blind deblurring on more realistic human and nature images. Our experimental results will be compared with the recovered results generated by the three above-mentioned methods, i.e., Xu and Jia [11], Krishnan et al. [9] and Pan et al. [59]. The estimation of the blur kernel for each method is implemented directly using the codes and parameter settings provided by the authors. As shown in Figure 10, the recovered result generated by Krishnan et al. [9] suffers from the over-smoothing of detailed texture structures due to the inaccurate estimation of the blur kernel. The loss of geometrical structures easily makes the deblurred image look less natural, resulting in significant visual quality degradation. The local magnification views shown in Figure 10 visually illustrate that the proposed method yields a more "natural-looking" performance compared with Xu and Jia [11] and Pan et al. [59]. In particular, more geometrical details in the face and hand regions could be preserved in our proposed method. The sharp edges, slightly over-smoothed by Xu and Jia [11] and Pan et al. [59], have been reconstructed accurately using our proposed method. Its superior performance benefits from the hybrid blur kernel constraints and edge-preserving TGV 2 regularizer. The excellent deblurring performance of our proposed method could also be visually found in Figure 11. As shown in the local magnification views, the inaccurate estimation of the blur kernel for Xu and Jia [11] causes the degraded visual quality with a significant loss of fine details. Krishnan et al. [9] is able to guarantee the quality of blur kernel estimation in this case. However, the final deblurring result tends to be unsatisfactory because the restored "number" could not be preserved correctly, shown in the local magnification views. The main geometrical structures and fine details in the recovered images are preserved by Pan et al. [59] and our proposed method. The ringing-like artifacts could be visually found near the ear region in the deblurring result by Pan et al. [59]. Our proposed method is able to overcome this limitation, but still generates slight ringing-like artifacts in the jaw region. The reason may be that our non-blind deconvolution method summarized in Algorithm 3 is performed with a constant regularization parameter λ. To further enhance the image quality, the regularization parameter should be selected spatially variant to suppress the ringing-like artifacts. More realistic deblurring results on human images are visually displayed in Figure 12. It could be observed that our proposed method yields deblurring results that are visually comparable with the current state-of-the-art methods. Figure 13 illustrates the realistic deblurring results on five different natural images. The sizes of the estimated blur kernels from top to bottom are 35 × 35, 55 × 55, 41 × 41, 35 × 35 and 55 × 55, respectively. Since these realistic images contain sufficient textures and geometrical structures, all competing methods have the capacity of accurately estimating the blur kernels in these cases. Therefore, our experimental results are visually comparable to others. The quality of the deblurred images could be correspondingly guaranteed in practical applications. Figure 11. Restoration of a blurred image with a blur kernel of the size 27 × 27 (top) and the corresponding local magnification views (middle and bottom). From left to right: (a) input blurred images, deblurred versions generated by (b) Xu and Jia [11], (c) Krishnan et al. [9], (d) Pan et al. [59] and (e) our proposed method, respectively.

Conclusions and Future Work
The major contributions of this work are mainly two-fold. First, a hybrid regularization method was developed to robustly estimate the blur kernel by incorporating both the L 1 -norm of kernel intensity and the squared L 2 -norm of the intensity derivative. The underlying assumption behind the proposed method was that the blur kernel was not only spatially sparse, but also piecewise smooth within the support of a continuous curve. An alternating direction algorithm was then proposed to effectively solve the resulting problem of blur kernel estimation. Second, to guarantee high-quality non-blind deconvolution, the TGV 2 -regularized variational model with an L 1 -norm data-fidelity term was presented to enhance the final image quality. The resulting optimization problem was then effectively solved using an ADMM-based numerical method. Comprehensive experiments implemented on both synthetic and realistic blurred images have illustrated the effectiveness of the proposed method. Given the recent progress in image deblurring, the proposed deblurring framework has several potential limitations in its current version. To further improve the blind deblurring performance, there is a huge potential to extend our future work along the following directions: • The constant parameters (i.e., η 1 and η 2 ) for both the L 1 -norm of kernel intensity and the squared L 2 -norm of intensity derivative in (3) are manually selected in our current work. Essentially, it is necessary to automatically and adaptively select the parameters according to the statistical properties of the blur kernel. For instance, if the blur kernel can be better sparsely represented in the spatial domain, η 1 should be larger; whereas η 2 plays a more important role if the blur kernel has a significant piecewise smooth structure. In our future work, an automatic estimation method should be developed to adaptively select the weighting parameters η 1 and η 2 in (3) to enhance the accuracy of blur kernel estimation. • The single-image blind deblurring method proposed in this work is performed based on a common assumption that the blur kernel is uniform (i.e., spatially invariant) across the image plane. Recent work in the literature [2,[60][61][62][63][64][65] has illustrated that the uniform simple assumption does not always hold in practice. To further enhance image quality, the assumption of the non-uniform (i.e., spatially variant) blur kernel has gained increasing attention in modern imaging sciences. In our opinion, the proposed hybrid regularized blur kernel estimation method discussed in Section 2 can be naturally extended to the case of non-uniform deblurring in future work.
As discussed beforehand, our proposed method suffers from some potential limitations (i.e., constant weighting parameters and uniform blurring assumption). Numerous experiments implemented on both synthetic and realistic blurred images have demonstrated its satisfactory deblurring performance. Therefore, it is still worthy of consideration since it could guarantee reliable performance compared to current state-of-the-art uniform blind deblurring methods. Recent research [9,41,59] has indirectly shown that our proposed method could be easily extended to the case of the non-uniform scenario. We believe there is a great potential for restoring blurred images using the proposed method in practical applications.