A Simple Local Minimal Intensity Prior and An Improved Algorithm for Blind Image Deblurring

TL;DR

This paper introduces a simplified sparsity prior called patch-wise minimal pixels (PMP) for blind image deblurring, along with an efficient algorithm that outperforms existing methods in stability, quality, and computational speed.

Contribution

The work proposes a novel sparsity prior (PMP) and an improved algorithm that avoids approximations, enhancing efficiency and effectiveness in blind image deblurring.

Findings

The proposed algorithm achieves better stability than state-of-the-art methods.

It demonstrates superior deblurring quality and robustness.

It is more computationally efficient than existing approaches.

Abstract

Blind image deblurring is a long standing challenging problem in image processing and low-level vision. Recently, sophisticated priors such as dark channel prior, extreme channel prior, and local maximum gradient prior, have shown promising effectiveness. However, these methods are computationally expensive. Meanwhile, since these priors involved subproblems cannot be solved explicitly, approximate solution is commonly used, which limits the best exploitation of their capability. To address these problems, this work firstly proposes a simplified sparsity prior of local minimal pixels, namely patch-wise minimal pixels (PMP). The PMP of clear images is much more sparse than that of blurred ones, and hence is very effective in discriminating between clear and blurred images. Then, a novel algorithm is designed to efficiently exploit the sparsity of PMP in deblurring. The new algorithm flexibly imposes sparsity inducing on the PMP under the MAP framework rather than directly uses the half quadratic splitting algorithm. By this, it avoids non-rigorous approximation solution in existing algorithms, while being much more computationally efficient. Extensive experiments demonstrate that the proposed algorithm can achieve better practical stability compared with state-of-the-arts. In terms of deblurring quality, robustness and computational efficiency, the new algorithm is superior to state-of-the-arts. Code for reproducing the results of the new method is available at https://github.com/FWen/deblur-pmp.git.

Figures

Tables

Algorithm 1: Latent image estimation
Input: Blurred image $B$, interim kernel estimation $k^i$.
$\beta \leftarrow {\beta }_0$, $I^0\leftarrow B$.
While $\beta \le {\beta }_{max}$ do ($t=0,1,2,\cdots$)
$I^{t+1,0}\leftarrow I^t$.
For $j=0:J-1$ do
Compute the mask $M^{t+1,j}$ via (13) based on $I^{t+1,j}$.
Obtain ${\tilde{I}}_s^{t+1,j}$ via (12) and further update ${\tilde{I}}^{t+1,j}$
via (14).
Compute gradient thresholding to obtain $G^{t+1,j+1}$
via (16).
Update $I^{t+1,j+1}$ via (18).
End for
$I^{t+1}\leftarrow I^{t+1,J}$.
$\beta \leftarrow a\beta$.
End while
$I^{i+1}\leftarrow I^{t+1}$.
Output: Intermediate latent image estimation $I^{i+1}$.

Algorithm 2: Blind blur kernel estimation
Input: Blurred image $B$, kernel initialization $k^0$ from the estimation in the last coarser-scale.
For $i=1:max\_iter$ do
Estimate $I^i$ via Algorithm 1 using $k^{i-1}$.
Estimate $k^i$ via (20).
End for
$\hat{k}\leftarrow k^i$, $\hat{I}\leftarrow I^i$.
Output: Kernel estimation $\hat{k}$, intermediate image $\hat{I}$.

Algorithm 3: Latent image estimation via solving (22)
Input: Blurred image $B$, interim kernel estimation $k^i$.
$\rho \leftarrow {\rho }_0$, $I^0\leftarrow B$.
While $\rho \le {\rho }_{max}$ do ($t=0,1,2,\cdots$)
Compute $Z^{t+1}$ via (23) using $I^t$.
$\beta \leftarrow {\beta }_0$, $I^{t+1,0}\leftarrow I^t$.
While $\beta \le {\beta }_{max}$ do ($j=0,1,2,\cdots$)
Obtain $G^{t+1,j+1}$ via (24) using $I^{t+1,j}$.
Solve (25) to update $I^{t+1,j+1}$.
$\beta \leftarrow a\beta$.
End while
$I^{t+1}\leftarrow I^{t+1,J}$.
$\rho \leftarrow a\rho$.
End while
$I^{i+1}\leftarrow I^{t+1}$.
Output: Intermediate latent image estimation $I^{i+1}$.

Table 4. TABLE I: Quantitative evaluation of the new method versus patch size on the dataset [ 50 ] (average PSNR and SSIM).

Patch size	PSNR (dB)	SSIM
$0.015\cdot mean(m,n)$	29.5891	0.8754
$0.020\cdot mean(m,n)$	29.6397	0.8778
$0.025\cdot mean(m,n)$	29.9764	0.8944
$0.030\cdot mean(m,n)$	29.7355	0.8837
$0.035\cdot mean(m,n)$	29.7088	0.8822

Table 5. TABLE II: Average PSNR and SSIM on deblurring 12 samples of the dataset [ 50 ] associated with the kernels 8, 9 and 10.

Kernel size	111	121	131	141	151
Pan et al. (PSNR)	22.40	22.85	21.71	23.13	21.69
Ours (PSNR)	23.05	23.35	23.42	23.39	22.49
Pan et al. (SSIM)	0.7503	0.7651	0.7087	0.7612	0.7141
Ours (SSIM)	0.7540	0.7708	0.7806	0.7619	0.7148

Table 6. TABLE III: Quantitative results on the dataset of Kohler et al. [ 50 ] , including the average PSNR and average SSIM.

Method	PSNR (dB)	SSIM
Cho and Lee	28.9831	0.8746
Xu and Jia	29.5373	0.8851
Shan et al.	25.8912	0.7748
Fergus et al.	22.7303	0.7048
Krishnan et al.	25.7246	0.7608
Whyte et al.	27.8441	0.8116
Hirsch et al.	26.8388	0.8095
Pan et al.	29.9513	0.8853
Ours	29.9764	0.8944

Table 7. TABLE IV: Quantitative results on the dataset of Levin et al. [ 2 ] , including the average PSNR and average SSIM.

Method	PSNR (dB)	SSIM
Levin et al.	31.1372	0.8960
Fergus et al.	29.4629	0.8451
Cho and Lee	30.7927	0.8837
Xu and Jia	31.3604	0.9083
Pan et al.	31.7297	0.9148
Ours	32.4450	0.9344

Table 8. TABLE V: Quantitative results on the dataset [ 7 ] , including the average PSNR and average SSIM.

Method	PSNR (dB)	SSIM
Cho and Lee	26.7548	0.8224
Cho et al.	17.9599	0.4922
Krishnan et al.	23.4366	0.7571
Levin et al.	25.4989	0.8079
Xu and Jia	29.1466	0.8553
Pan et al.	30.2502	0.8587
Ours	30.5660	0.8593

Table 9. TABLE VI: Runtime comparison in seconds (the kernel size is fixed at 51 × 51 51 51 51\times 51 for each algorithm).

Method	$256\times 256$	$512\times 512$	$800\times 800$
Xu et al. (C++)	1.05	2.43	5.35
Levin et al. (Matlab)	155.9	657.8	1598.6
Pan et al. (Matlab)	162.5	548.6	1261.3
HQS (Alg. 3) (Matlab)	107.3	401.2	938.8
Ours w/o PMP (Matlab)	5.59	14.51	35.03
Ours ($r=4$) (Matlab)	38.88	132.7	295.3
Ours ($r=8$) (Matlab)	16.39	55.73	134.8
Ours ($r=16$) (Matlab)	11.02	41.69	94.74
Ours ($r=0.025m$) (Matlab)	18.61	44.36	95.55

Equations

$$B=k\otimes I+n,$$

$$\mathcal{P}(I)(i)=\underset{(x,y)\in {\Omega }_i}{\min }\left(\underset{c\in \{r,g,b\}}{\min }I(x,y,c)\right),$$

$$\mathcal{P}(B)\ge \mathcal{P}(I).$$

$$I_p:={\mathcal{P}}^T(\mathcal{P}(I))=I\circ M,$$

$$\mathcal{D}(I)(i,j)=\underset{(x,y)\in {\Omega }_{i,j}}{\min }\left(\underset{c\in \{r,g,b\}}{\min }I(x,y,c)\right),$$

$$\underset{k,I}{\min }\mathcal{L}(k\otimes I,B)+\gamma G(k)+\mu R(I),$$

$$\underset{k,I}{\min }{\Vert k\otimes I-B\Vert }_2^2+\gamma {\Vert k\Vert }_2^2+\mu {\Vert \nabla I\Vert }_0.$$

$$\begin{array}{rl}& \underset{k,I}{\min }{\Vert k\otimes I-B\Vert }_2^2+\gamma {\Vert k\Vert }_2^2+\mu {\Vert \nabla I\Vert }_0\\ & \text{subject to }\mathcal{P}(I)(i)\sim p(x),\mathrm{for}i\in \{1,\cdots ,P\}.\end{array}$$

$$\underset{k,I}{\min }{\Vert k\otimes I-B\Vert }_2^2+\gamma {\Vert k\Vert }_2^2+\mu {\Vert \nabla I\Vert }_0+\alpha {\Vert \mathcal{P}(I)\Vert }_0,$$

$$\begin{array}{rl}& \underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\mu {\Vert \nabla I\Vert }_0\\ & \text{subject to }\mathcal{P}(I)(i)\sim p(x),\mathrm{for}i\in \{1,\cdots ,P\}.\end{array}$$

$$\begin{array}{rl}& \underset{I,G}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\beta {\Vert \nabla I-G\Vert }_2^2+\mu {\Vert G\Vert }_0\\ & \text{subject to }\mathcal{P}(I)(i)\sim p(x),\mathrm{for}i\in \{1,\cdots ,P\},\end{array}$$

$$\begin{split}&\tilde{I}_{s}^{t+1,j}(i)=\left\{\!{\begin{array}[]{*{20}{l}}\!{0,}&{|I_{s}^{t+1,j}(i)|<\lambda}\!\\ \!{I_{s}^{t+1,j}(i),}&{{\rm{otherwise}}}\!\end{array}}\!\right.\!,\\ &~{}~{}{\rm{for}}~{}i\in\{1,\cdots,P\}.\end{split}$$

$${M^{t+1,j}}(i,j)=\left\{\!{\begin{array}[]{*{20}{l}}\!{1,}&{{\rm{if}}~{}(i,j)\in{\Omega^{t+1,j}}{\rm{}}}\!\\ \!{0,}&{{\rm{otherwise}}}\!\end{array}}\!\right.\!.$$

$${\tilde{I}}^{t+1,j}=I^{t+1,j}\circ (1-M^{t+1,j})+{\mathcal{P}}^T({\tilde{I}}_s^{t+1,j}),$$

$$G^{t+1,j+1}=\arg \underset{G}{\min }\beta {\Vert \nabla {\tilde{I}}^{t+1,j}-G\Vert }_2^2+\mu {\Vert G\Vert }_0,$$

$$\begin{split}&{G^{t+1,j+1}}(i,j)=\left\{\!{\begin{array}[]{*{20}{l}}\!{0,}&\!{{{(T(i,j))}^{2}}<{\mu}/{\beta}}\!\\ \!{T(i,j),}\!&\!{{\rm{otherwise}}}\!\end{array}}\!\right.\!,\\ &~{}~{}{\rm{with}}~{}~{}T=\nabla{\tilde{I}^{t+1,j}}.\end{split}$$

$$I^{t+1,j+1}=\arg \underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\beta {\Vert \nabla I-G^{t+1,j+1}\Vert }_2^2,$$

$$I^{t+1,j+1}={\mathcal{F}}^{-1}\left(\frac{\overline{\mathcal{F}(k^i)}\circ \mathcal{F}(B)+\beta \left(\overline{\mathcal{F}({\nabla }_h)}\circ \mathcal{F}(G_h^{t+1,j+1})+\overline{\mathcal{F}({\nabla }_v)}\circ \mathcal{F}(G_v^{t+1,j+1})\right)}{\overline{\mathcal{F}(k^i)}\circ \mathcal{F}(k^i)+\beta \left(\overline{\mathcal{F}({\nabla }_h)}\circ \mathcal{F}({\nabla }_h)+\overline{\mathcal{F}({\nabla }_v)}\circ \mathcal{F}({\nabla }_v)\right)}\right).$$

$$k^{i+1}=\arg \underset{k}{\min }{\Vert k\otimes (\nabla I^i)-\nabla B\Vert }_2^2+\gamma {\Vert k\Vert }_2^2.$$

$$\begin{array}{rl}& k^{i+1}=\\ & {\mathcal{F}}^{-1}\left(\frac{\overline{\mathcal{F}({\nabla }_h{I}^i)}\circ \mathcal{F}({\nabla }_{h}B)+\overline{\mathcal{F}({\nabla }_v{I}^i)}\circ \mathcal{F}({\nabla }_{v}B)}{\overline{\mathcal{F}({\nabla }_h{I}^i)}\circ \mathcal{F}({\nabla }_h{I}^i)+\overline{\mathcal{F}({\nabla }_v{I}^i)}\circ \mathcal{F}({\nabla }_v{I}^i)+\gamma }\right).\end{array}$$

$$\underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\mu {\Vert \nabla I\Vert }_0+\alpha {\Vert \mathcal{P}(I)\Vert }_0.$$

$$\begin{array}{rl}\underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\mu {\Vert G\Vert }_0+\alpha {\Vert Z\Vert }_0& +\beta {\Vert \nabla I-G\Vert }_2^2\\ & +\rho {\Vert \mathcal{P}(I)-Z\Vert }_2^2,\end{array}$$

$$\begin{split}&{Z^{t+1}}(i,j)=\!\left\{\!{\begin{array}[]{*{20}{l}}\!{0,}&\!{{{(Y(i,j))}^{2}}<{\alpha}/{\rho}}\!\\ \!{Y(i,j),}&\!{{\rm{otherwise}}}\!\end{array}}\!\right.\!,\\ &~{}~{}{\rm{with}}~{}~{}Y={\cal P}({I^{t}}),\end{split}$$

$$\begin{split}&{G^{t+1}}(i,j)=\!\left\{\!{\begin{array}[]{*{20}{l}}\!{0,}&\!{{{(T(i,j))}^{2}}<{\mu}/{\beta}}\!\\ \!{T(i,j),}&\!{{\rm{otherwise}}}\!\end{array}}\!\right.\!,\\ &~{}~{}{\rm{with}}~{}~{}T=\nabla{I^{t}}.\end{split}$$

$$\underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2+\beta {\Vert \nabla I-G^{t+1}\Vert }_2^2+\rho {\Vert \mathcal{P}(I)-Z^{t+1}\Vert }_2^2.$$

$$\begin{array}{rl}\underset{I}{\min }{\Vert k^i\otimes I-B\Vert }_2^2& +\beta {\Vert \nabla I-G^{t+1}\Vert }_2^2\\ & +\rho {\Vert I_p-{\mathcal{P}}^T(Z^{t+1})\Vert }_2^2.\end{array}$$

$$\begin{array}{rl}\underset{I_p}{\min }{\Vert k^i\otimes I_p-(B-k^i\otimes {\breve{I}}_p)\Vert }_2^2+& \beta {\Vert \nabla I_p-(G^{t+1}-\nabla {\breve{I}}_p)\Vert }_2^2\\ +& \rho {\Vert I_p-{\mathcal{P}}^T(Z^{t+1})\Vert }_2^2,\end{array}$$

$$\underset{{\breve{I}}_p}{\min }{\Vert k^i\otimes {\breve{I}}_p-(B-k^i\otimes I_p)\Vert }_2^2+\beta {\Vert \nabla {\breve{I}}_p-(G^{t+1}-\nabla I_p)\Vert }_2^2.$$

$${\Vert {\mathcal{D}}_{I^t}(I)-u\Vert }_2^2\approx {\Vert I-\left({\mathcal{D}}_{I^t}^T(u)+{\breve{I}}_d^t\right)\Vert }_2^2,$$

$$\begin{array}{rl}& {\Vert {\mathcal{D}}_{I^t}^T({\mathcal{D}}_{I^t}(I))-{\mathcal{D}}_{I^t}^T(u)\Vert }_2^2\\ & ={\Vert I_d-{\mathcal{D}}_{I^t}^T(u)\Vert }_2^2={\Vert I-\left({\mathcal{D}}_{I^t}^T(u)+{\breve{I}}_d\right)\Vert }_2^2,\end{array}$$

Full text

I Introduction

Blind image deblurring, also known as blind deconvolution, aims to recover a sharp latent image from its blurred observation when the blur kernel is unknown. It is a fundamental problem in image processing and low level vision, which has been extensively studied and is still a very active research topic in image processing and computer vision.

For single image deblurring and under the assumption that the blur is uniform and spatially invariant, the blur process can be modeled as a convolution operation [1], given by

[TABLE]

where $B$, $k$, $I$, $n$, and $\otimes$ denote the blurred image, blur kernel, latent (clear) image, additive noise, and convolution operator, respectively. In the blind deblurring problem, only the blurred image $B$ is known a prior. The objective is to recover the kernel $k$ and the clear image $I$ simultaneously from $B$. Basically, it is a highly ill-posed problem as there exist many different solution pairs of $(k,I)$ giving rise to the same $B$. Note that a typical undesired solution is that $I=B$ and $k$ being a delta blur kernel.

To make the blind deblurring problem well-posed, image prior and blur kernel model exploitation is the key of most effective methods. Well developed image priors include the gradient sparsity prior [5, 2, 3, 4, 1], normalized sparsity prior [6], patch prior [7], group sparsity prior [8], intensity prior [9], dark channel prior [10, 11], extreme channel prior [12], latent structure prior [13], local maximum gradient prior [14], class-specific prior [15], and learned image prior using a deep network [16], to name just a few. Meanwhile, blur kernel models include the non-uniform model with blur from multiple homographies [17, 18, 19, 20], depth variation model [21, 22], in-plane rotation model [23], and forward motion model [24]. Most of these methods exploit image prior and blur kernel model under the maximum a posterior (MAP) framework. Generally, since the related deblurring problems are non-convex [25], incorporating regularization to exploit image prior and/or kernel model helps to effectively increase the probability of achieving a good local solution. In addition, heuristic edge selection is also an effect way to help the MAP algorithms to avoid undesired trivial solutions. For a more detailed discussion, see [26, 27].

While image gradient sparsity is a popular and commonly used prior, image intensity or gradient based priors have shown good complementary effectiveness when jointly used with the image gradient prior [10, 9, 12, 14]. As priors and models designed for natural images usually do not generalize well to specific images such as text images [28], face images, and low-light images [29], simultaneously exploiting multiple priors has the potential to achieve satisfactory performance on both natural and specific images [10, 9, 12, 14].

Though the sophisticated priors [10, 12, 14] have shown promising effectiveness, there exist two limitations: i) Jointly using multiple priors makes the corresponding algorithms computationally expensive. ii) Since these priors involved subproblems in the corresponding algorithms cannot be solved explicitly, non-rigorous approximate solution is commonly used, which limits the best exploitation of the capability of such priors. These limitations motivate us to develop a more effective and efficient method in this work. The main contributions are as follows.

I-A Contribution

Firstly, we propose a novel local intensity based prior, namely the patch-wise minimal pixels (PMP) prior. The PMP is a collection of local minimal pixels. Intuitively, since the blur process has a smoothing effect on the image pixels, the intensity of a local minimal pixel would increase after the blur process. As a result, the PMP of clear images are much sparser than those of blurred ones. The PMP metric is as simple as the direct intensity prior, but is very effective in discriminating between clear and blurred images. It can be viewed as a simplification of the dark channel prior in [10], which facilitates efficient computation while being effective. A more detailed comparison with existing intensity priors [10, 9] is provided in Section II-B.

Secondly, a novel algorithm is proposed to exploit the sparsity of PMP under the MAP framework. The new algorithm flexibly imposes sparsity inducing on the PMP of the latent image in the MAP deblurring process. Compared with existing algorithms solving augmented MAP formulations directly based on half quadratic splitting, e.g., [10, 12, 14], the proposed algorithm greatly improves the computational efficiency in substance. More importantly, while the algorithms [10, 12, 14] use approximate solution in handling non-explicit subproblems, the new algorithm avoids such approximation. As a consequence, in comparison with state-of-the-art methods, it is practically more robust and can achieve competitive performance on both natural and specific images.

Finally, extensive experimental results on blind image deblurring have been provided to evaluate the performance of the proposed method. The results show that the proposed method can achieve state-of-the-art performance on both natural and specific images. In terms of the deblurring quality, robustness and computational efficiency, the proposed method is superior to the compared algorithms.

I-B Related Work

In recent years, single image blind deblurring has made much progress with the aid of various effective priors on images and blur kernels [30]. Most works are based on the variational Bayes and MAP frameworks [5, 26, 27, 31, 2, 6, 3, 32, 10, 9, 4, 11, 1, 12, 16, 33, 34, 35, 36]. Typically, such a blind deblurring method generally has two steps. First, blur kernel is estimated from the blurred observation under the MAP framework. Then, based on the estimated blur kernel, the latent sharp image is estimated via non-blind deconvolution methods, e.g., [37, 38, 39].

As the naive MAP method could fail on natural images, exploiting the statistical priors of natural images and selection of salient edges are the key of the success of state-of-the-art methods. The gradient sparsity prior of natural images is the most widely used prior. But it has been shown in [2] that, in practice, the methods using the gradient sparsity prior in the MAP framework favor blurry images rather than clear ones. This limitation can be mitigated by techniques as explicit sharp edge pursuit [5, 27, 33, 40] or heuristic edge selection [26]. However, the assumption of such techniques that there exist strong edges in the latent images may not always be satisfied.

There also exist various other image priors designed to reinvigorate MAP, such as normalized sparsity prior [6], internal patch recurrence [32], and direct intensity prior [9]. Though effective for either natural or specific images, such priors usually cannot yield satisfactory performance on both natural and specific images. The recently proposed dark-channel prior [10] and data driven learned prior [16] can achieve satisfactory performance on both natural and specific images, but the involved optimization algorithms are computationally expensive.

Particularly, the intensity based priors considered in [10, 9] are close to our proposed PMP prior. As PMP is computed based on local minimal intensities, it is fundamentally different from the intensity priors in [10, 9]. A detailed explanation on the difference is provided in Section II-B. Furthermore, our algorithm exploits the PMP sparsity prior in a different way from that in [10, 9] and, as a consequence, it can reduce the computational complexity significantly while has more robust (stable) performance. A detailed comparison on the algorithms is provided in Section IV.

Recently, data driven methods have also made much success with the aid of powerful deep learning techniques, e.g., [23, 41, 42, 43, 44, 45]. For example, Sun et al. [23] endeavored to employ a convolutional neural network (CNN) to estimate and remove non-uniform motion blur. Nah et al. [42] proposed a multi-scale CNN to recover the latent image in an end-to-end manner without any assumption on the blur kernel. Meanwhile, spatially variant recurrent neural network and scale-recurrent network have been designed for deblurring in [41, 45]. Moreover, Kupyn et al. [44] proposed an end-to-end learned method for deblurring based on conditional generative adversarial networks (GAN). In addition, end-to-end CNN model for video deblurring has been considered in [43]. Though these data driven methods can yield favorable performance in various scenarios, their success depends heavily on the consistency between the training data and the test data. This leads to the limitation of their generalization capability.

Outline: The rest of this paper is organized as follows. Section II introduces the sparsity property of PMP and discusses the connection between PMP and existing intensity priors. The new algorithm is developed in detail in Section III. Section IV presents comparison between the proposed algorithm and existing related algorithms. Section V provides experimental results. Finally, this paper concludes with a brief summery in Section VI.

Notations: $\otimes$ and $\nabla$ denote the convolution and gradient operation, respectively, $\left\lceil\cdot\right\rceil$ denotes the ceil operator, $\circ$ stands for Hadamard (element-wise) product, $\bar{x}$ denotes the conjugate of a complex quantity $x$. $X(i,j)$ denotes the $(i,j)$-th element of a matrix $X$, ${\cal F}(X)$ denotes the 2-D FFT of $X$, and ${\cal F}^{-1}(X)$ denotes the 2-D inverse FFT of $X$.

II Patch-Wise Minimal Pixels

This section first introduces the proposed PMP prior and presents analysis on its statistic property. Then, comparison with existing intensity priors is provided.

II-A Patch-Wise Minimal Pixels

PMP is a collection of local minimal pixels over non-overlapping patches. Let an image $I\in{\mathbb{R}^{m\times n\times c}}$ be divided into $P$ non-overlapped patches with a patch size of $r\times r$, for which $P=\left\lceil{\frac{m}{r}}\right\rceil\!\cdot\!\left\lceil{\frac{n}{r}}\right\rceil$. The PMP is defined as

[TABLE]

for $i=1,2,\cdots P$, where ${\Omega_{i}}$ denotes the index set of the pixel locations of the $i$-th patch. It is easy to see that ${\cal P}(I)\in{\mathbb{R}^{P}}$ which contains patch-wise (local) minimal pixels of $I$.

In what follows, we show that the PMP of blurred images are much less sparse than those of natural clear images. Fig. 1 compares the histogram statistic of PMP between clear images and their blurred counterparts over more than 5000 natural images from the VGG111http://www.robots.ox.ac.uk/$\sim$vgg/data/ dataset. The blurred images are synthesized from the clear ones using the blur kernels of the dataset [2]. It can be seen from Fig. 1 that the PMP of clear natural images have significantly more zero elements than those of blurred images. The PMP of clear images under a threshold (e.g., 0.9) follow a hyper Laplacian distribution and manifest a sparsity property. This sparsity property of PMP provides a natural metric to discriminate clear images from blurred ones. Based on this property, the proposed algorithm imposes sparsity inducing on PMP in the deblurring process to achieve more accurate kernel estimation.

Besides the empirical results shown in Fig. 1, the following result theoretically shows that blurred images have less sparse PMP than their clear counterparts.

Property 1: Let ${\cal P}(B)$ and ${\cal P}(I)$ denote the PMP of the blurred and clear images, respectively. Then

[TABLE]

This property can be directly derived via extending Property 1 in [10]. It indicates that the blur process increases the values of PMP, which gives rise to that the PMP of blurred images are less sparse than their clear counterparts.

In the following, without loss of generality, consider $c=1$ for simplicity. For the PMP ${\cal P}(I):{\mathbb{R}^{m\times n}}\to{\mathbb{R}^{P}}$ defined in (2), we further define its inverse operation for later use. Its inverse operation is defined by its transpose ${{\cal P}^{T}}(z):{\mathbb{R}^{P}}\to{\mathbb{R}^{m\times n}}$ for any $z\in{\mathbb{R}^{P}}$. Accordingly, we have

[TABLE]

where $M\in\mathbb{R}^{m\times n}$ is the binary mask corresponding to the PMP subset of $I$, with $M(i,j)=1$ if $I(i,j)$ is a minimal pixel of a patch and $M(i,j)=0$ otherwise.

II-B Comparison with Existing Intensity Sparsity Priors

Intensity sparsity has also been considered in [10, 9] for blind deblurring. However, the proposed PMP is fundamentally different from the intensity metrics considered in [10, 9], which is explained as follows.

The closest work exploiting intensity sparsity is the dark channel metric considered in [10]. For an image $I\in{\mathbb{R}^{m\times n\times c}}$, the dark channel is defined as

[TABLE]

for $i=1,2,\cdots m$ and $j=1,2,\cdots n$, where ${\Omega_{i,j}}$ denotes the index set of the pixel locations of the patch centered at the $(i,j)$-th pixel. It is easy to see that ${\cal D}(I)\in{\mathbb{R}^{m\times n}}$. While the dark-channel is computed in a convolution like manner with an output of size $m\times n$, the proposed PMP is computed on non-overlapped patches with a vector output ${\cal P}(I)\in{\mathbb{R}^{P}}$ of size $P=\left\lceil{\frac{m}{r}}\right\rceil\!\cdot\!\left\lceil{\frac{n}{r}}\right\rceil$ for a patch size of $r\times r$. As a result, the proposed PMP is much simpler than the dark-channel prior, thereby facilitating the design of more efficient algorithm.

Moreover, the work [9] uses sparsity inducing directly on the intensity of the image for text image deblurring. Since the intensity distribution of text images is close to two-tone, using $\ell_{0}$-regularization to promote the intensity sparsity has demonstrated outstanding effectiveness in text image deblurring. However, the distribution of the intensity values of natural images are more complex than that of text images, and the direct intensity sparsity is not applicable to natural images.

In comparison, the proposed PMP metric is as simple as the direct intensity metric [9], but is very effective in discriminating between clear and blurred natural images as shown in Fig. 1.

III Proposed Deblurring Algorithm Using PMP Sparsity Regularization

This section presents an efficient algorithm via flexibly incorporating the sparsity regularization of PMP into the conventional MAP framework. The new algorithm is a variant of the half quadratic splitting algorithm, but is different to the direct half quadratic splitting algorithm used in [10, 9, 12, 14], which is explained in detail in Section IV.

Recall that the well-known MAP formulation is given by

[TABLE]

where $\gamma$ and $\mu$ are positive weight parameters, and ${\cal L}$ is a data fidelity term, which restricts $k\otimes I$ to be consistent with the blurred image $B$. To make the problem well-posed, $G$ and $R$ are penalty functions to exploit the structures in the blur kernel and the latent image, respectively.

With the nature of that the gradient of natural images is sparse, $R(I)$ is usually selected as the $\ell_{0}$-norm penalty of $\nabla I$ (the gradient of $I$). Meanwhile, selecting both the loss function ${\cal L}$ and the penalty for the kernel as the ${\ell_{2}}$-norm yields

[TABLE]

The ${\ell_{2}}$-norm is not only optimal for Gaussian noise, but also enables the development of efficient algorithms because it facilitates fast computation of the involved subproblems via fast Fourier transform (FFT).

To further exploit the sparsity of the PMP of the latent image, e.g., ${\cal P}(I)\in{\mathbb{R}^{P}}$ for a patch size of $r\times r$, now consider a constrained extension of (7) as

[TABLE]

As introduced in Section II, $p(x)$ is a probability density function of a hyper Laplacian distribution for $x$ below a threshold such as 0.9. The constrained problem (8) is nonsmooth and nonconvex. Similar to most existing deblurring algorithms in solving MAP-like objective functions, we propose an efficient algorithm to solve (8) via alternatingly update the blur kernel and the latent image. In the proposed algorithm, the constraint in (8) is approximately imposed via iteratively sparsity inducing on ${\cal P}(I)$ in the latent image subproblem.

Note that a natural alternative of (8) to promote sparsity of the PMP in the MAP framework is the formulation as follows:

[TABLE]

where $\alpha$ is a positive weight parameter and the last term uses $\ell_{0}$-norm penalty to achieve sparsity inducing on the PMP of the latent image. This formulation can be solved directly by the half quadratic splitting algorithm in an alternating manner similar to the algorithms in [10, 9, 12, 14]. However, we show that the proposed algorithm for solving (8) is superior to the direct half quadratic splitting algorithm for solving (9), which will be explained in Section IV in detail.

III-A Estimating the Latent Image

Given an interim estimation of the blur kernel, denoted by $k^{i}$, the latent image is updated via optimizing the following problem:

[TABLE]

Using an auxiliary variable $G$ with respect to the image gradient $\nabla I$, the problem (10) can be approximated by

[TABLE]

where $\beta>0$ is a sufficient large penalty parameter such that it enforces $\left\|{\nabla I-G}\right\|_{2}^{2}\approx{\rm{0}}$, and hence $G\approx\nabla I$.

Without the constraint, such as in the case of the traditional MAP algorithm, the problem (11) can be typically solved using the block coordinate descent algorithm, which alternatingly updates the two variables $I$ and $G$. The proposed algorithm also solves (11) via alternating between the variables $I$ and $G$ in which the constraint is approximately imposed.

Specifically, since the constraint in fact imposes a sparsity regularization on the PMP of $I$, we use a simple thresholding/shrinkage step in the iteration procedure to impose sparsity inducing on the PMP of $I$. Given $I^{t}$ and at the $(t+1)$-th iteration of the latent image subproblem, denote the PMP subset of $I^{t}$ by $I_{s}^{t}:={\cal P}({I^{t}})$, we iteratively impose thresholding on $I_{s}^{t}$ and update $I$ and $G$ via the following steps.

First, let $\lambda>0$ be a threshold parameter. The PMP is thresholded as

[TABLE]

Subsequently, let ${\Omega^{t+1,j}}$ denote the index set of the PMP in ${I^{t+1,j}}$ and define the mask corresponding to the PMP subset to be

[TABLE]

Then, we update ${I^{t+1,j}}$ as

[TABLE]

where ${{\cal P}^{T}}$ is the inverse operation of ${{\cal P}}$ defined in Section II-A. With ${\tilde{I}^{t+1,j}}$ given in (14), the gradient-subproblem solves the following formulation

[TABLE]

which is a proximal minimization and from [46] the solution is explicitly given by

[TABLE]

Finally, the latent image is updated via solving the following problem

[TABLE]

which can be efficiently computed by means of FFT as (18) given in the next page, where $\nabla=({\nabla_{h}},{\nabla_{v}})$ and $G=({G_{h}},{G_{v}})$ are used, such that they correspond to image gradients in the horizontal and vertical directions, respectively.

This algorithm is summarized in Algorithm 1, which contains two loops. In Algorithm 1, $a$ is a positive increasing factor, which is set to $a=2$ in the experiments. Extensive numerical studies show that the inner loop usually converges within a few iterations. For example, we use $J=3$ in the experiments in Section V.

III-B Estimating the Blur Kernel

Similar to other existing state-of-the-art algorithms, the kernel estimation is performed in the gradient space. As it has been shown that gradient space based methods are more accurate than intensity space based ones [27, 3, 4]. Specifically, given an interim estimation of the latent image, denoted by $I^{i}$, the blur kernel is updated via solving

[TABLE]

Due to its quadratic form, the solution can be efficiently computed by means of FFT, give by

[TABLE]

Moreover, the estimated kernel is further refined via setting the negative elements to zero and normalization. In practical implementation, the multi-scale deconvolution scheme [27] is adopted to estimate the kernel in a coarse-to-fine manner. The main steps for kernel estimation at a single scale level are shown in Algorithm 2.

III-C Implementation Tricks

To make the augmented objective function in (11) accurately approaching that in (10), a sufficiently large value of $\beta$ is desired, ideally $\beta\to\infty$. However, with a very large value of $\beta$, an alternating algorithm directly minimizing (11) would be very slow and impractical. To address this problem, a standard trick is to use a continuation process for $\beta$. In other words, one starts with a properly small value of $\beta$ and gradually increase it by iteration until a large target value is reached. This trick is used in Algorithm 1 with a $a>1$.

The thresholding step of the PMP in (12) corresponds to a nonconvex $\ell_{0}$-regularization. In addition, the objective function in (11) is also nonconvex. Hence, with different initialization and/or parameter setting, a nonconvex algorithm would end up with one of its many local minimizers. In view of this, in implementing Algorithm 1, we use soft-thresholding instead of the hard-thresholding in the first few scales in the multi-scale procedure and then turn back to the hard-thresholding. As the soft-thresholding corresponds to the convex $\ell_{1}$-regularization, this “first loose and then tight” strategy makes the proposed algorithm more stable and performing satisfactorily.

IV Comparison with the Half Quadratic Splitting Algorithm Solving the Regularized MAP Formulation (9)

As mentioned in Section III, a natural alternative of (8), which can incorporate sparsity inducing of the PMP into the MAP framework, is given by (9). Compared with the formulation (8), the formulation (9) is even more explicit and can be solved by means of the half quadratic splitting algorithm. However, in this section we show that, compared with the direct half quadratic splitting algorithm solving (9), the proposed algorithm is not only superior in computational complexity, but also can avoid non-rigorous approximate solution in solving the regularized MAP problem.

To solve (9) with a given interim kernel estimation $k^{i}$, the latent image problem becomes

[TABLE]

Similar to [13], [14], using two auxiliary variables $G$ and $Z$ with respect to $\nabla I$ and ${\cal P}(I)$, respectively, the problem (21) is approximated by

[TABLE]

where $\beta$ and $\rho$ are positive penalty parameters. Given ${I^{t}}$ at the $(t+1)$-th iteration, the $G$- and $Z$-subproblems are proximity operators, which can be efficiently solved in an element-wise manner as

[TABLE]

and

[TABLE]

Then, the $I$-subproblem becomes

[TABLE]

In view of that there exist two augmentation terms in the nonconvex problem (22), to make the algorithm practically working well, a standard trick is to use a continuation process for each of $\beta$ and $\rho$ similar to the algorithms in [10, 9]. In such a manner, the main steps of the half quadratic splitting algorithm are sketched in Algorithm 3.

Although both Algorithms 1 and 3 contain two main loops, the former is much more efficient than the latter in practice. That is because the penalty parameters ${\rho_{max}}$ and ${\beta_{max}}$ in Algorithm 3 should be chosen sufficiently large to make (22) accurately approximating for (21), while a small value of $J$ (e.g., $J=3$) is enough for Algorithm 1 to give satisfactory performance.

Moreover, although the $I$-step in Algorithm 3 solves a quadratic problem (25), it cannot be efficiently solved via FFT similar to (18). A strategy to explicitly solve the least-square problem (25) in closed-form is to vectorize the variables and convert the convolution operation into linear multiplication. However, this least-square problem involves computing the inverse of high-dimensional matrices of size $(mn)\times(mn)$ with $m\times n$ be the size of $I$. Thus, it is computationally expensive to handle practical-sized inputs. Meanwhile, since ${\cal P}(I)$ is a subsampling function of $I$ and only contains a subset of the pixels of $I$, clearly, the problem (25) cannot be solved in closed-form via FFT similar to (18) and the algorithms as given in [6, 47, 32].

With the definition in (4), problem (25) can be equivalently rewritten as

[TABLE]

Now, let ${{\breve{I}}_{p}}=I\circ(1-M)$ be the complementary set of $I_{p}$ such that it satisfies ${I_{p}}+\breve{I}_{p}=I$. It is easy to see from (25) that $I_{p}$ and $\breve{I}_{p}$ are coupled through the kernel convolution operation. With this in mind, to enable FFT based efficient solution, we can iteratively solve (26) via alternating between $I_{p}$ and $\breve{I}_{p}$. For example, we first fix $\breve{I}_{p}$ to solve $I_{p}$ by

[TABLE]

and then fix $I_{p}$ to solve $\breve{I}_{p}$ by

[TABLE]

Thanks to the quadratic form of (26), iteratively solving (27) and (28) is guaranteed to converge to the global minimizer of (26) with any bounded initialization. Even though both (27) and (28) can be efficiently solved by means of FFT, iteratively solving them makes Algorithm 3 having three iteration loops, and hence results in an increase of computational complexity in terms of runtime.

Note that in the dark-channel based method [10], the $I$-subproblem has a similar formulation as (25), where ${\cal P}(I)$ is replaced by the dark-channel extraction operation, which is solved via FFT in close-form. The close-form solution is derived via implicitly using an approximation. That is

[TABLE]

where ${{\cal D}_{{I^{t}}}}$ denotes the dark-channel extraction operator based on ${I^{t}}$, ${\cal D}_{{I^{t}}}^{T}$ is the inverse operator of ${{\cal D}_{{I^{t}}}}$, and $\breve{I}_{d}^{t}$ is the complementary set of $I_{d}^{t}$ with $I_{d}^{t}$ being the subset pixels of $I^{t}$ which forms the dark-channel. In fact, the right term in (29) can be viewed as an approximation of

[TABLE]

where ${\cal D}_{{I^{t}}}^{T}({{\cal D}_{{I^{t}}}}(I))={I_{d}}$ and $I={I_{d}}+\breve{I}_{d}$ are used.

Similar approximation has also been used in [12, 14]. In comparison, the proposed Algorithm 1 completely avoids such approximation and would be more stable in practical applications, as demonstrated by the results in Fig. 2, 3 and Table II in the experiments.

V Experimental Results

We firstly investigate the robustness of the new algorithm in terms of sensitivity against the kernel size parameter, in comparison with the most close method [10]. Then, we evaluate the proposed algorithm on three benchmark datasets in comparison with state-of-the-art blind image deblurring methods. Furthermore, we conduct evaluation on face, natural, text, and low-light images. Matlab code for reproducing the results of the new algorithm is available at https://github.com/FWen/deblur-pmp.git.

For the new algorithm, $\mu=4\times{10^{-3}}$, $a=2$, $J=3$, ${\beta_{0}}=2\mu$, and ${\beta_{max}}={10^{5}}$ are used. For each scale, we use $max\_iter=5$ as a trade-off between accuracy and speed. The threshold parameter for PMP is initially set to $\lambda=0.1$ and gradually reduced to the mean of the PMP values. Note that, in general the optimal values of these parameters depend on the distribution of the clear images, the kernel models, the optimization algorithm used for the non-convex deblurring problem. Hence, it is difficult to select the optimal values of these parameters. Like most blind deblurring methods, such as the ones compared in the sequel, these parameters are selected empirically in practice.

The “first loose and then tight” strategy introduced in Section III-C is employed to make the algorithm performing practically well. The patch size is set dependent on the image size as $r=0.025\cdot mean(m,n)$. Similar to [5, 26, 2, 10], we first estimate the blur kernel by the proposed algorithm, and then obtain the final latent image based on the estimated kernel by a non-blind deblurring method. The non-blind deblurring algorithm [9] is employed for the final latent image estimation. The performance of the compared algorithms is evaluated in terms of peak-signal-to-noise ratio (PSNR), structural similarity (SSIM) [48] and cumulative error ratio of the deblurred images and kernel estimation.

Table I presents a quantitative evaluation of the proposed method versus the patch size $r$ on the dataset [50]. It can be seen that the new method is robust to the patch size.

V-A Robustness: Sensitivity to the Kernel Size Parameter

As discussed in Section IV, the new algorithm can avoid the non-rigorous approximation in solving non-explicit priors (e.g., dark-channel) involved subproblems in [10], [12], [14]. This brings the new algorithm a potential advantage of being more stable in practical applications. To illustrate this point, the first experiment compares the performance of the new algorithm with Pan et al. [10] through investigating their sensitivity against the kernel size parameter. The selection of the kernel size parameter has a substantial influence on the performance of most deblurring algorithms.

Table II and Fig. 2 present the PSNR and SSIM results of the two compared methods versus kernel size on deblurring three challenging kernels in the dataset [50] (the kernels 8, 9 and 10). These selected kernels are the most challenging kernels in the dataset [50], which have sizes of more than 100 pixels. There are four samples for each of the three kernels. Moreover, Fig. 3 shows deblurred results of the two algorithms when using different values of the kernel size parameter on three realistic images. It can be seen from Fig. 2, 3, and Table II that our method is less sensitive to the kernel size variation, which demonstrates its better robustness in practice.

V-B Evaluation on Benchmark Datasets

The second experiment uses the dataset by Kohler et al. [50], which contains 48 blurred samples corresponding to 4 clear images and 12 blur kernels. The compared algorithms include Cho and Lee [27], Xu and Jia [26], Shan et al. [31], Fergus et al. [5], Krishnan et al. [6], Whyte et al. [47], Hirsch et al. [49], and Pan et al. [10]. Fig. 4 presents a statistical analysis of the PSNR and SSIM results of the compared algorithms on deblurring the 48 blurred images. Table III shows the average PSNR and average SSIM of the algorithms. The PSNR and SSIM of each deblurred image are computed via comparing it with 199 clear images captured within the camera motion trajectory. The results of the methods [5, 26, 27, 31, 6, 47, 49] are those reported in [50], while the result of the method [10] is computed from the deblurred results provided by the authors at their website222http://vllab1.ucmerced.edu/$\sim$jinshan/projects/dark-channel-deblur/.

It can be seen that the new algorithm can achieve state-of-the-art performance in terms of the PSNR and SSIM results. Fig. 5 presents visual comparison on four challenging images with heavy blurs from the dataset [50], including the ‘Blurry1_8’, ‘Blurry2_9’, ‘Blurry3_10’, and ‘Blurry4_11’ images. It can be inferred from Fig. 5 that the proposed algorithm can achieve comparable or even better visual results compared with existing state-of-the-art methods [3, 10].

Fig. 6 further investigates the effectiveness of the proposed PMP regularization to show the results of the new algorithm with and without PMP regularization on the benchmark dataset [50]. The results demonstrate that the PMP regularization gives rise to distinct PSNR and SSIM improvement.

The third experiment uses the dataset by Levin et al. [2], which contains 32 blurred samples corresponding to 4 clear images and 8 blur kernels. The parameter $\mu$ is set to $5\times{10^{-3}}$ for all examples. Fig. 7 shows the cumulative error ratios of the compared algorithms, which are computed based on the sum of square difference (SSD) error as follows. Firstly, for a restored image, the SSD error is computed as the SSD between it and its clear counterpart using the best shift between them. Then, this SSD error is normalized with respect to the SSD error of the de-convolution result using the ground-truth kernel, which results in an error ratio. Finally, from the error ratios of all the deblurred images, the cumulative error ratio is computed to get the success rate. It is empirically noticed that deblurred results with error ratios above 2 are visually implausible. Table IV compares the average PSNR and average SSIM of the methods, whilst Fig. 8 presents statistical analysis of the PSNR and SSIM results

It can be observed from Table IV, Fig. 7 and Fig. 8 that, the proposed algorithm achieves the best PSNR and SSIM results while attaining 100% success at an error ratio of 2 on the dataset [2]. Fig. 9 illustrates the estimated kernels by the proposed algorithm on this dataset, whilst Fig. 10 presents the deblurred results of the proposed algorithm on four challenging samples in this dataset.

The fourth experiment further considers a much larger dataset of Sun et al. [7], which contains 640 blurred samples corresponding to 80 clear images and 8 blur kernels. Fig. 11 shows the cumulative error ratios of the algorithms on this dataset. The compared algorithms include Cho and Lee [27], Cho et al. [33], Krishnan et al. [6], Levin et al. [2], Xu and Jia [26], and Pan et al. [10]. Table V compares the average PSNR and average SSIM of the algorithms, whilst Fig. 12 presents statistical analysis of the PSNR and SSIM results. It can be seen that the proposed algorithm can achieve the state-of-the-art performance in terms of cumulative error ratio, PSNR and SSIM of revcovery.

V-C Evaluation on Natural and Specific Images

In what follows, the proposed method is further evaluated on face, natural, text, and low-light images. Here we only provide some typical results for each class due to limited space. More samples are provided online at https://github.com/FWen/deblur-pmp.git.

Face image: Face image deblurring is challenging for methods developed for natural images, since the lack of edges and textures in face images makes accurate kernel estimation challenging. Fig. 13 compares the proposed method with the methods [3, 10] on two realistic blurred face images. The results demonstrate that our method compares favorably or even better against the methods give in [3, 10].

Natural image: The results of the compared algorithms on two real natural images are shown in Fig. 14. Again, our algorithm compares competitively against the methods [3, 10]. It can be see that the proposed PMP regularization helps to significantly reduce the ringing artifacts in the deblurred image, which makes the proposed algorithm yielding state-of-the-art performance.

Text image: The results of the compared algorithms on two text images are illustrated in Fig. 15. Our algorithm performs comparably with the method [10]. When without using the PMP regularization, our algorithm may fail to reconstruct the correct blur kernel and yields a result with heavy ringing artifacts.

Low-light image: Low-light images usually cannot be well handled by most deblurring methods. A main reason is that low-light images often have saturated pixels which interfere with the kernel estimation process [29, 51]. Fig. 16 presents the results on a low-light image. The state-of-the-art low-light image deblurring method [29] is used in the comparison. Compared with the method [29] specifically designed for low-light images, our method gives a comparable result.

V-D Computational Complexity

Finally, we compare the computational complexity of our algorithm ($r=0.025\cdot mean(m,n)$) with those of the algorithms [3, 10, 4]. Our algorithm without PMP and our algorithm with different patch size $r\in\{4,8,16\}$ are also compared. As explained in Section IV, the proposed Algorithm 1 would be more efficient than the traditional half quadratic splitting (HQS) algorithm, e.g. Algorithm 3. To substantiate this, the runtime of the HQS algorithm using the PMP (Algorithm 3) is also compared (with $r\!=\!0.025\!\cdot\!mean(m,n)$). The simulations are conducted under Windows 10 on a desktop PC with an Intel Core i7-4790 CPU at 3.6 GHz with 16 GB RAM. For our method and the methods [10, 4], the runtime of the non-blind deblurring step is included in the results. Among these algorithms, one can observe from Table VI that the algorithm developed by Xu et al. [3] with C++ implementation is the fastest. However, in some cases, its restoration quality is inferior to our algorithm as illustrated earlier in the above figures. Our algorithm is much faster than the algorithms [10, 4]. Compared with the algorithm of Pan et al. [10], our algorithm is more than an order of magnitude faster. Note that, the algorithm [10] can be accelerated in the dark-channel computation step as mentioned by the authors. The result of this algorithm presented here is obtained by the code provided by the authors at their website, which is implemented without such acceleration.

VI Conclusion

This work proposed a local minimal intensity based prior, namely PMP, and an improved algorithm for blind image deblurring. The prior is simple yet effective in discriminating between clear and blurred images. Rather than directly using the half quadratic splitting algorithm, the new algorithm flexibly imposes sparsity inducing on the PMP in the deblurring procedure under the MAP framework. Particularly, it avoids non-rigorous approximate solution in existing algorithms in jointly handling multiple non-explicit priors, while being much more efficient. Extensive experiments on both natural and specific images demonstrated that it not only can achieve state-of-the-art deblurring quality, but also can improve the practical stability and computational efficiency substantially. In brief, in terms of both the practical robustness and computational efficiency, the proposed algorithm is superior to the compared algorithms in this work.