J Math Imaging Vis DOI 10.1007/s10851-014-0537-9

Motion Deblurring Using Non-stationary Image Modeling Wen-Ze Shao · Qi Ge · Hai-Song Deng · Zhi-Hui Wei · Hai-Bo Li

Received: 23 May 2013 / Accepted: 8 September 2014 © Springer Science+Business Media New York 2014

Abstract It is well-known that shaken cameras or mobile phones during exposure usually lead to motion blurry photographs. Therefore, camera shake deblurring or motion deblurring is required and requested in many practical scenarios. The contribution of this paper is the proposal of a simple yet effective approach for motion blur kernel estimation, i.e., blind motion deblurring. Though there have been proposed several methods for motion blur kernel estimation in the literature, we impose a type of non-stationary Gaussian prior on the gradient fields of sharp images, in order to automatically detect and purse the salient edges of images as the important clues to blur kernel estimation. On one hand, the prior is able to promote sparsity inherited in the non-stationarity of the precision parameters (inverse of variances). On the other hand, since the prior is in a Gaussian form, there exists a great possibility of deducing a conceptually simple and computationally tractable inference scheme. Specifically, the

well-known expectation–maximization algorithm is used to alternatingly estimate the motion blur kernels, the salient edges of images as well as the precision parameters in the image prior. In difference from many existing methods, no hyperpriors are imposed on any parameters in this paper; there are not any pre-processing steps involved in the proposed method, either, such as explicit suppression of random noise or prediction of salient edge structures. With estimated motion blur kernels, the deblurred images are finally generated using an off-the-shelf non-blind deconvolution method proposed by Krishnan and Fergus (Adv Neural Inf Process Syst 22:1033–1041, 2009). The rationality and effectiveness of our proposed method have been well demonstrated by the experimental results on both synthetic and realistic motion blurry images, showing state-of-the-art blind motion deblurring performance of the proposed approach in the term of quantitative metric as well as visual perception.

W.-Z. Shao (B) · Q. Ge College of Telecommunications and Information Engineering, Nanjing University of Posts and Telecommunications, Nanjing, China e-mail: [email protected]

Keywords Blind motion deblurring · Camera shake · Blur kernel estimation · Expectation–maximization · Split Bregman iteration

W.-Z. Shao Department of Computer Science, Technion-Israel Institute of Technology, Haifa, Israel H.-S. Deng School of Mathematics and Statistics, Nanjing Audit University, Nanjing, China Z.-H. Wei School of Computer Science and Technology, Nanjing University of Science and Technology, Nanjing, China H.-B. Li School of Computer Science and Communication, KTH Royal Institute of Technology, Stockholm, Sweden

1 Introduction As a fundamental problem in the fields of image processing and computational photography, blind motion deblurring has been attracted much attention and widely investigated in the past several years [1–9], aiming at restoring a √ √ M× M from its blurry and noisy version clear image f ∈ R √ √ √ √ y ∈ R M× M when the motion blur kernel k ∈ R D× D is unknown. Without loss of generality, the standard motion blurring degradation model can be formulated in the convolution form by

123

J Math Imaging Vis

y = k ⊗ f + n,

(1)

where ⊗ denotes the convolution operation, and n is additive white noise, or in the matrix-vector form by y = Kf + n,

(2)

where y, f, and n are M × 1 lexicographically ordered vectors of the intensities of the observed blurry image y, original image f , and additive white noise n, and the M × M matrix K is block-circulant with circulant blocks (BCCB) corresponding to k, meaning that the periodic boundary condition is used in (1). Besides, this paper assumes Gaussian statistics for the noise n, i.e., n∼N (0, α −1 I M×M ), where α −1 is the noise variance and α is called the precision parameter. Motion deblurring (1) or (2) is a notoriously ill-posed inverse problem, even though in the case that the blur kernel k is assumed to be known. Therefore, proper regularization should be imposed in order to obtain reasonable estimates for the blur kernels and the original images. Inspired by that, there have been proposed a number of blind motion deblurring methods in the literature most of which essentially originate from the Bayesian framework. For instance, Fergus et al. [1] use a mixture-of-Gaussians (MoG) prior for modeling natural image statistics and a mixture-of-exponentials (MoE) prior for modeling motion blur kernels with parameters involved in both priors learned in advance, and use the variational Bayesian (VB) approach of Miskin and Mackay [10] to estimate the posteriors of both image gradient fields and blur kernels. In the final, the original images are estimated by the Richardson–Lucy algorithm [23,24] and histogram equalization. Based on Fergus et al.’s original work, Levin et al. [5] provide a more thorough and profound analysis on blind motion deblurring, pointing out that the variational Bayesian methods are significantly more robust than the maximum a posteriori (MAP) ones. Additionally, a slightly simpler numerical scheme is provided for Fergus et al. [1] by Levin et al. [6]. As for the final image estimation, [6] exploits a sparse prior-based non-blind deblurring principle, which is solved by the iteratively reweighted least squares (IRLS) method [22]. And empirical studies in [6] conform with the theoretical analysis in [5]; that is, [6] outperforms the MAP methods of Shan et al. [2] and Cho and Lee [3]. Nevertheless, it has to be noted that the variational Bayesian inference for blind deblurring with sparse or superGaussian priors remains a challenging obstacle [6,9]. Along this line, a new blind motion deblurring approach is recently proposed by Babacan et al. [9] imposing a kind of general sparse priors on the image gradient fields. Interestingly, the empirical studies in [9] show that the non-informative prior p(ξ ) ∝ ξ −1 is the most effective one, where ξ is a random variable. Besides, Wipf and Zhang revisit in-depth the VB blind deconvolution methods more recently [27], providing a counter-intuitive perspective that the optimal image priors

123

for blind deblurring need not reflect the statistics of natural images, which seems to be consistent with the empirical studies in [9]. Another commonly exploited blind deconvolution principle is the maximum a posteriori or penalized least squares [2,3,8,11–14]. Similar to variational Bayesian approaches [1,6,9], many methods of this type also exploit sparse priors [2,8,14,26] imposed on both the image gradient fields and motion blur kernels. For example, Cai et al. [14] introduce a sparse representation based motion deblurring method, the sparse priors of which are based on the framelet transform. Krishnan et al. [8] propose a normalized sparsity-based motion deblurring method with motion blur kernels regularized by the Laplace priors, and claim that their algorithm achieves similar performance to that of Fergus et al. [1], but significantly outperforms that of Shan et al. [2]. More recently, Xu et al. [26] propose a gradually approximate l 0 sparse representation method for natural image deblurring, and the motion blur kernels are solved by the ridge regression. Despite of that, there are several methods [3,11–13] not explicitly exploiting sparse priors. As a matter of fact, the success of these deblurring methods usually builds on two important properties of their estimation processes, namely, explicit suppression of random noise in homogenous regions and explicit prediction of salient edges in the images. For example, prior to the blur kernel estimation, Cho and Lee [3], and Xu and Jia [11] have utilized bilateral filtering or Gaussian filtering together with shock filtering to predict sharp edges. Our main contribution in this paper is the proposal of a simple yet effective approach for motion blur kernel estimation. The new approach is regularized by a type of nonstationary Gaussian (NSG) prior, which is essentially placed on the sparse gradient fields of natural images. The benefits of NSG prior are at least twofold. On one hand, the prior is capable of promoting sparsity which is inherited in the non-stationarity of the precision parameters in NSG, and therefore helpful to automatically purse the salient edges in images as important clues to blur kernel estimation. Due to the non-stationarity of the precision parameters, the proposed approach has essentially exploited a sparse prior in an implicit way, although the NSG prior is expressed in a non-sparse Gaussian form. On the other hand, since NSG is in a Gaussian form, there is a great possibility of deducing a conceptually intuitive and computationally tractable inference scheme. In specific, the well-known expectation– maximization (EM) algorithm [28] is exploited to alternatingly estimate the motion blur kernels, the image gradients and the parameters in the NSG prior. In difference from many existing methods, no hyperprior is imposed on any parameter in this paper; there are not any pre-processing steps involved in our proposed approach, either, such as explicit suppression of random noise in homogenous areas or prediction of

J Math Imaging Vis

salient edges in the images. With estimated motion blur kernels, deblurred images are finally generated using an off-theshelf non-blind deconvolution method [20]. The rationality and effectiveness of the proposed approach have been well demonstrated by the experimental results on both synthetic and real-world motion blurry images, showing the state-ofthe-art performance in the term of quantitative metric as well as visual perception. The rest of the paper is organized in the following. In Sect. 2, an EM-based numerical algorithm is proposed for motion blur kernel estimation by exploiting the NSG image prior, with finally deblurred images produced using the fast non-blind deconvolution method in [20]. Section 3 provides an in-depth discussion on the relationships between the proposed approach and state-of-the-art methods, MAP and VB mechanisms both included. In Sect. 4, the experimental results on both synthetic and real-world motion blurry images are provided and analyzed, with comparison against classical and state-of-the-art blind deblurring methods [1,3,5,6,8,25]. Finally, Sect. 5 concludes the paper.

likelihood function with respect to {fd }d∈|| , β, k in the following p({yd }d∈ |{fd }d∈ ; β, k) − ||M 2

= (2π )

β

||M 2



 β 2 yd − Kf d 2 , exp − 2

(5)

d∈

which essentially corresponds to the classical Richardson– Lucy blind deblurring algorithm [25]. As noted earlier, (4) is an ill-posed inverse problem due to the ill-conditionedness of K T K. Besides, the blind nature of (4) makes it more difficult and complex to implement the inversion process. Actually, this is why a large amount of image restoration literature has been built on the MAP and VB principles, which provide flexible frameworks to incorporate appropriate prior information into the estimation process. Taking into account of both the calculation convenience and modeling flexibility, this paper imposes a type of nonstationary Gaussian (NSG) prior (6) on the gradient fields of natural images  N (fd |0, ( d )−1 ) p({fd }d∈ ; {γ d }d∈ ) = d∈

2 Blind Motion Deblurring with Non-stationary Gaussian Image Modeling

=

M  

N (fd (m)|0, (γ d (m))−1 ),

d∈ m=1

(6)

2.1 Modeling for Motion Blur Kernel Estimation As several previous methods [1,2,5,6,8], this subsection focuses on motion blur kernel estimation in the derivative space. The reason lies in that, on one hand, the derivative representation fits better the NSG prior which assume independent relationships among the derivative components; on the other hand, a derivative deconvolution system is to be less illconditioned as the regularization is imposed on the unknowns themselves rather than their derivatives [6]. Specifically, we formulate (1) in the derivative space using the horizontal and vertical derivative filters ∇h = [1, −1], ∇v = [1, −1]T as follows yd = k ⊗ f d + n d ,

(3)

where yd = ∇d ⊗ y, f d = ∇d ⊗ f, n d = ∇d ⊗ n, d ∈ ,  = {h, v}, and the cardinality of  is 2, i.e., || = 2. For the convenience of presentation in the following text, the vectorized version of (3) is formulated as yd = Kf d + nd .

(4)

Here, we still assume that nd is of Gaussian statistics satisfying nd ∼N (0, β −1 I M×M ), where β is the precision hyperparameter. By now, our goal becomes estimating k and {fd }d∈ as well as β using the blurry and noisy {yd }d∈ . The simplest principle used for estimating {fd }d∈ , β, k is the maximum likelihood (ML), that is, by maximizing the

where γ d = (γ d (1), . . . , γ d (M))T are the precision para

meters, and  d = diag{γ d }. As a matter of fact, Vera et al. [15] propose a similar modeling scheme for nonblind image deconvolution, in which, however, the precision parameters γ d are estimated using a computationally intensive method, called type-II ML or evidence analysis (EA). In this paper, (6) is exploited to the problem of blind motion deblurring with {fd }d∈ , {γ d }d∈ , k, β estimated by the well-known EM algorithm. The advantage of exploiting (6) is at least twofold. On one hand, compared against the highly sparse priors in [1,6,9] as well as the normalized sparsity measure in [8], the theoretical analysis in Sect. 3 plus the empirical studies in Sect. 4 shows that the non-stationary hyper-parameters γ d make the NSG prior capable of adaptively achieving the sparsity promotion; that is, salient edges in the images can be automatically estimated as important clues to motion blurkernel estimation. On the other hand, the convexity and smoothness of the logarithm of the NSG prior lead to a computationally convenient scheme for blur kernel estimation. 2.2 Optimization Using Expectation–Maximization By now, we come to the issue of estimating {fd }d∈ , {γ d }d∈, β, k with the likelihood function (5) and the NSG prior (6),

123

J Math Imaging Vis

That is, the final estimates of {γ d }d∈ , β, k can be formulated as {γ dEA }d∈ , k EA , β EA







log β −1 I + K( d )−1 K T

= arg min {γ d }d∈ ,k,β



d∈

+ydT β −1 I + K( d )−1 K T

Fig. 1 A graphical model that describes the dependencies among random variables of the proposed model. Circular nodes denote hidden random variables, square nodes correspond to parameters of the model, and the observed random variables are represented by triple circled nodes

the graphical model of which is shown in Fig. 1. Notice that the image gradient fields {fd }d∈ are the hidden random variables, {yd }d∈ are the observed random variables, and the motion blur kernel k, the precision parameters {γ d }d∈ , β are all the deterministic model parameters. Based on the Bayes law, the posterior distribution of the hidden variables {fd }d∈ can be formulated as p({fd }d∈ |{yd }d∈ ; {γ d }d∈ , β, k) p({yd }d∈ |{fd }d∈ ; β, k) p({fd }d∈ ; {γ d }d∈ ) . = p({yd }d∈ ; {γ d }d∈ , β, k)

K T yd , μd = β d  −1  =  d + βK T K . d

= ln p({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k) p({f } |{y } ;{γ (t) } ,β (t) ,k (t) ) d d∈ d d∈ d d∈

    β ||M yd − Kf d 22 ln β − = 2 2 d∈  M M 1 1 2 + ln γ d (m) − γ d (m)(fd (m)) +const 2 2 m=1 m=1  β  ||M yd − Kf d 22 ln β − = 2 2



N (yd |0, β −1 I + K( d )−1 K).

d∈

123

+ −

ln γ d (m)

M   1 γ d (m) (fd (m))2 + const. 2

(12)

d∈ m=1

In (12), the expected values · are taken with respect to (t) p({fd }d∈ |{yd }d∈ ; {γ d }d∈ , β (t) , k (t) ), and 2     yd − Kf d 22 = yd − K(t) μ(t) d 2    T (t) K(t) K(t) , + trace



 (10)

1 2

M 

d∈ m=1

(9)

d∈

=

d∈

(8)

p({yd }d∈ ; {γ d }d∈ , β, k)  = p(yd |fd ; β, k)p(fd ; γ d )dfd

(11)

Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k)

(7)

If the parameters {γ d }d∈ , β, k have been provided, then the hidden variables {fd }d∈ can be estimated as (8). In practice, {γ d }d∈ , β, k have to be estimated along with {fd }d∈ , and one possible estimation strategy is the evidence analysis by maximizing the marginal likelihood function

 yd .

Nonetheless, direct optimization of (11) presents several computational difficulties, since its derivatives with respect to the model parameters {γ d }d∈ , β, k are difficult to compute. Besides, the minimization problem requires a constrained EA , k EA have to optimization algorithm since {γ EA d }d∈ , β be positive. Instead of (11), the EM algorithm is used in this section to overcome the computational difficulties. After initializing the model parameters {γ d }d∈ , β, k to some val(1) (1) ues {γ (1) d }d∈ , β , k , the algorithm proceeds by alternatingly performing the E-step and M-step as follows. E-step Based on (5) and (6), compute the expected value of the logarithm of the complete likelihood:

According to Tipping’s relevance vector machine [17], it , β, k) is a multivariis not hard to obtain that p(fd |yd ; γ d ate Gaussian distribution N (fd |μd , d ) with parameters defined as 

−1

d

 (t) 2 (fd (m))2 = (μ(t) (m,m), d (m)) + d

(13) (14)

J Math Imaging Vis

T (t)  (t) (t) K μ(t) = β yd , d d −1   T (t) (t) (t) (t) (t) K = d + β K , d

(15) (16)

(t) where μ(t) d and d are respectively the expectation and covariance matrix of the random variable fd . In practice, (15) can be efficiently computed using the conjugate gradient (CG) method just similar to[1,5,6,9]. Besides, for computational efficiency of (14), (t) d is approximated by a diagonal matrix via inverting only the diagonals of  (t) d +  T β (t) K(t) K(t) . M-step With Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k), the (t+1) model parameters γ d , β (t+1) , k (t+1) can now be updated via solving the following problem (t+1)

{γ d

}d∈ , β (t+1) , k (t+1)

= arg

max

{γ d }d∈ ,β,k

Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k). (17)

Based on (12) and (17), it is easy to deduce ∂ Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k) ∂β  1  ||M yd − Kf d 22 , − = 2β 2

(18)

(19)

Setting (18), (19) to zero and using (13), (14), we obtain the (t+1) updated {γ d }d∈ , β (t+1) as follows (t+1)

γd

(m) =

β (t+1)

1

(20) , (m, m) ||M =    .    (t) 2 (t) (t) (t) T K(t) yd − K μd  + trace d K (t)

(μd (m))2 +

(t) d

2

d∈

(21) To deduce an updated estimate for k (t+1) , we first borrow the idea of Levin et al. [5,6] and rewrite (12) in a computationconvenient form as Q ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k)  β   T (t) ||M T ln β − = k Ad k − 2(b(t) d ) k 2 2 M 1 ln γ d (m) 2 d∈ m=1

M   1 γ d (m) (fd (m))2 + const, 2 d∈ m=1

d

(t)

bd (i) =

M 

(t)

μd (m + i) yd (m).

(23) (24)

m=1 (t)

where Ad ∈ R D×D denotes the covariance matrix of all √ √ (t) the D × D windows in fd , bd denotes the correlation with yd , and i, j represents the kernel indexes. To be noted that, (23) is efficiently (t) computed in a similar way to (14) in practice; that is, d is approximated by a diagonal matrix. Then, the derivative of Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k) with respect to the motion blur kernel k can be obtained as Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k) ∂k    (t) (t) = Ad k − 2(bd )T . (25) Thus, estimating k (t+1) becomes a simple quadratic program which can be solved efficiently, following which√k (t+1)√ is   D D projected onto the constraint set C = k|k ≥ 0, i=1 j=1  k(i, j) = 1 . 2.3 Implementation Details Basically, motion blur kernel estimation can be implemented by alternatingly computing formulas (16), (15), (21), (20), (25). However, it should be noted here that motion blur-kernel estimation relies heavily on the estimation accuracy of β. Nonetheless, it is empirically found that (21) is not much effective in practice, direct use of which may sometimes lead to visually unpleasant deblurred images. Through a multitude of trials, a workable formula for adaptively updating β is finally found which is a slightly modified version of (21), given as follows β (t+1) =

||M    ,    (t) 2 (t) (t) (t) T K(t) K η+ yd − K μd  + trace d d∈

d∈

d∈

−

(t) + (m + i, m + j) ,

2

(26)      T  (t) 2 (t) where η  c+ yd −K(t) μd  + trace d K(t) K(t) ,

(t)

+

m=1

d∈

d∈

∂ Q (t) ({yd }d∈ , {fd }d∈ ; {γ d }d∈ , β, k) ∂γ d (m)   1 = − (fd (m))2 . γ d (m)

where k ∈ R√D is√the vectorized version of the motion blur kernel k ∈ R D× D , and M  (t) (t) Ad (i, j) = μ(t) d (m + i) μd (m + j)

(22)

2

and it has been found that c = 2 × 10−2 works well for a great number of motion blurry images (Remark: The choice of η is relevant to ||, and in this paper || = 2). Actually, our empirical choice of (26) originates from the idea of variational Bayesian inference [16,17] as the hyperparameter β is imposed with a hyperprior. In one of our on-going

123

J Math Imaging Vis

work, blind motion deblurring is formulated completely in the framework of variational Bayes, which has achieved more advanced modeling and more effective blind deblurring performance. Now, motion blur kernel estimation is summarized as Algorithm 1, corresponding to the alternating iterations among formulas (16), (15), (26), (20), (25).

l p -norm based deconvolution method proposed by Krishnan and Fergus [20], because the method is very fast and has been found robust to small kernel errors [8]. Specifically, the method uses the split Bregman iteration scheme and the fast Fourier transform (FFT) to solve the following cost function   fˆ = arg min λ k ⊗ f − y22 + ∇h f  p + ∇v f  p , f

Algorithm 1 Motion Blur Kernel Estimation with Nonstationary Image Modeling Inputs. observed motion blurry image y , the blur kernel √ The√ size D × D , the iterations T = 10. (0) (0) Initializations. With β (0) = 1e2,  d = 1e4 ·I M×M,μd ,k (0), (0) (1) (1) using (16), (26), (20) and calculate d ,  (1) d ,β , k (25). while t < T do (t) (t) Step 1. Given the current estimates {γ (t) d }d∈ , β , k , (t) update d using (16); (t) (t) Step 2. Given the current estimates {γ (t) d }d∈ , β , k , (t) update μd using (15); (t) (t) (t) Step 3. Given the current estimates d , μd , k , update β (t+1) using (26); (t) (t) Step 4. Given the current estimates d , μd , update (t+1) using (20); γd (t) (t) Step 5. Given the current estimates d , μd , update k (t+1) by solving (25) using quadratic programming with projection onto the constraint set C ;

end while Another point to be noted is that, as for large-scale motion blur kernels, a large number of iterations are often required in order to converge to a reasonable solution [8]. Therefore, in the same way as many previous blind deblurring mechanisms, the idea of multi-scale implementation is borrowed in this paper to estimate motion blur kernels, i.e., using a coarse-to-fine pyramid of image resolutions.√The number of √ D × D and scales S is determined by the kernel size √ the size ratio of 2 between consecutive resolutions in each dimension, such that the kernel size at the coarsest scale is (0) 3 × 3. It is to be noted that, except the finest scale, μd and k (0) in Algorithm 1 are obtained by up-sampling from their estimates in the previous scale, and as for the finest scale, (0) they are simply set as μd = yd , k (0) = I√ D×√ D . 2.4 Non-blind Deconvolution With the estimated motion blur kernels, motion deblurring reduces to a non-blind deconvolution issue, and there have been proposed many non-blind deblurring algorithms in the literature, e.g., [18,19]. In this paper, we choose a simple

123

(27) where the regularization parameter λ controls the tradeoff between the data fidelity and constrained terms, k is the estimated motion blur kernel, 0 < p ≤ 1, and ∇h , ∇v are the same derivative filters as Sect. 2.1. 3 Reformulating and Understanding the Proposed Approach In this section, discussions on the relationships between the proposed method and several state-of-the-art MAP and VB approaches are provided. First of all, by reformulating the proposed approach into an energy minimization functional, we are to demonstrate its theoretical superiority over the current MAP approaches. On the second, the differences and similarities between the proposed approach and the VB methodologies are also discussed, demonstrating that our approach is simpler in both terms of problem modeling and algorithm deduction. According to the proposed numerical scheme in Sect. 2.2, (1) (1) upon given {γ (1) d }d∈ , β , k , the estimates of {μd }d∈ , {γ d }d∈ , β, k can be reformulated in the following form   {μd }d∈ , {γ d }d∈ , β, k = arg = arg

×

min

{μd }d∈ ,{γ d }d∈ ,β,k



G {μd }d∈ , {γ d }d∈ , β, k



min

{μd }d∈ ,{γ d }d∈ ,β,k

⎧ ||M ⎪ ⎪ ⎨ − 2 ln β +

β 2

  T    yd − Kμd 2 + trace dK K 2

d∈

M ⎪ 1   ⎪ ln γ d (m) + ⎩−2 d∈ m=1

1 2

M   d∈ m=1

γ d (m)



⎫ ⎪ ⎪ ⎬

 . 2  ⎪ μd (m) + d (m, m) ⎪ ⎭

(28) More specifically, the cost functions that are used to update (t) μd and k (t+1) in (15) and (20) are given by $ 2   (t) μd = arg min yd − K(t) μd  μd 2 % M 1  (t) 2 + (t) γ d (m)(μd (m)) , (29) β m=1 $   (t) 2 (t+1) = arg min yd − Kμd  k k 2 ⎛ ⎞⎫ (t) ⎬  +trace ⎝ K T K⎠ , (30) ⎭ d

J Math Imaging Vis (t) (t) where γ d , d , β (t) , k (t) are the previous estimates of  γ d , d , β, k. With (29) and (30), it comes to our mind that the MAP motion deblurring methods in [3,11] possess simpler updated equations as follows $  2   (31) fd(t) = arg min yd − K(t) fd  + ν fd 22 , 2 fd $    (t) 2 (32) k(t+1) = arg min yd − Kf d  +  k22 . k

2

where ν,  are positive tuning parameters. Nonetheless, the success of [3,11] is achievable only as (31) is assisted by additional steps of bilateral/Gaussian filtering and shock filtering in order to predict salient edges as key clues to motion blur kernel estimation. Note that, in our approach there are no any pre-processing operations involved to suppress random noise or predict salient edges. The reason is that the regularization term in (29) assures that the proposed approach is capable of automatically detecting and pursing the salient edges in images as important clues to blur kernel estimation. This advantage is essentially rooted in the non-stationarity of precision parameters {γ d }d∈ . Meanwhile, compare (30) against (32) and find that, the proposed approach has incorporated the uncertainty of fd into the blur kernel estimation  via the covariance matrix d . While, [3,11] have not, and in fact, it is also the case with the MAP methods [2,8,14,26] which explicitly employ the sparse priors. In this view, therefore, our proposed method is theoretically more robust than the current MAP approaches. As for differences between the proposed method and the VB principles, the most remarkable point is that no hyperpriors are imposed on any parameters in this paper, and there is no need to estimate the posterior distributions of the parameters involved as most VB principles [1,5–7,9,10] do. Therefore, our proposed approach is conceptually simpler compared against the VB blind motion deblurring methods in [1,6,9]. In spite of the differences in problem modeling and algorithm deduction, it is interesting to note that the final kernel estimation algorithms of the proposed and VB approaches [1,6,9] are similar in form as (29) and (30). Hence, it can be conjectured that our method has intense theoretical connections with VB principles. However, compared with most of the current VB approaches, e.g., [1,5–7,10], estimates of the precision parameters {γ d }d∈ , β are computationally distinct. Here, we take the influential MoG-based VB kernel estimation method proposed by Levin et al. [6] for example. In [6], {γ d }d∈ is estimated as

(t+1)

γd

(m) =

L  1 2 σ l=1 l

πl σl e L  l=1

−

πl σl e

(t) (t) (μd (m))2 + d (m,m) 2σl2

−

(t) (t) (μd (m))2 + d (m,m) 2 2σl

.

(33)

where πl , σl are the prior parameters in the mixture of L Gaussians. Clearly, (33) is more complex than (20) in both terms of deduction and computation. Furthermore, the addiL , have to be tional parameters involved in (33), i.e., {πl , σl }l=1 provided or estimated in advance, which is, however, avoided in the proposed method. As for β, in order to speed up convergence, the authors [6] start with a low precision parameter β0 in each scale and gradually increase it during iteration until a desired value of β is reached which is provided manually at the beginning of each scale. As stated in Sect. 2.3, in the proposed approach the precision parameter β has been updated automatically in each scale rather than provided in a hand-tuned way. To briefly sum up, the estimation accuracy of the motion blur kernel k in the proposed approach is highly relevant to that of the precision parameters {γ d }d∈ and β. Another key point of the proposed method is the incorporation of the uncertainty of fd into the process of blur kernel estimation.

4 Experimental Results In this section, experiment results on both synthetic and realworld motion blurry images are provided and analyzed. All the experiments are performed with MATLAB v7.0 on a portable computer equipped with an Intel i7-4600M CPU (2.90 GHz) and 8 GB memory, running Windows 7 (Service Pack 1). 4.1 Experiments on Levin et al.’s Benchmark Dataset [5] In this subsection, the proposed method is tested on the benchmark blurry image dataset proposed by Levin et al. [5], which can be downloaded from the link: www.wisdom.weizmann. ac.il/~levina/papers/LevinEtalCVPR2011Code.zip. The dataset contains 32 motion blurry images generated utilizing 4 base images of size 255 × 255 and 8 motion blur kernels of sizes ranging from 13 × 13 to 27 × 27, estimated by recording the trace of focal reference points on the boundaries of the original images [5,27]. Except for the motion blurry images, deblurred results corresponding to four different blind motion deblurring methods are also included in the dataset. These methods include: the stationary Gaussian image prior-based VB blind deblurring method in [5], the influential MoG image prior-based VB blind deblurring methods in [1,6], and the implicit sparse image priorbased MAP blind deblurring method in [3]. The metric of SSD (Sum of Squared Difference) as introduced in [5] is used to make evaluations which quantifies errors between the estimated and the original images. As advocated in [5] as well as many other state-of-the-art methods [6–9], the SSD ratio between the images deconvolved respectively with the estimated blur kernel and the ground truth blur-kernel is used

123

J Math Imaging Vis Fig. 2 The cumulative histogram of the SSD error ratios achieved by Fergus et al. [1], Cho and Lee [3], Levin et al. [5,6], and our proposed approach with final image deconvolution using [20]. The success percentages, i.e., SSD error ratio below 3, of different methods are: 69 % [1], 75 % [3], 19 % [5], 84 % [6], 94 % (Ours)

as the final evaluation measure. The underlying reason is to take into account of the fact that harder blur-kernels give a larger deblurring error even as the ground truth kernel is known since the corresponding non-blind deblurring problem is also harder. Another point to be noted is that, the finally deblurred images of different methods [1,3,5,6] included in the benchmark dataset are obtained using the non-blind deblurring method in [22] rather than [20] as in our proposed method. In fact, [22] utilizes a similar non-blind deconvolution model as [20], i.e., the formula (27), except that [22] also considers second-order derivative filters and solves (27) based on the iteratively reweighted least squares technique, which is, however, computationally more demanding than [20]. According to the above observation as well as for the completely fair assessment of different methods of motion blur kernel estimation, we use the estimated blur kernels as provided in the dataset and the non-blind deblurring method in [20] (with the same parameter settings) to regenerate the finally deblurred images for each of the four methods to be compared. Figure 2 shows the cumulative histogram of the SSD error ratios across 32 test images achieved by our proposed method and other four methods [1,3,5,6]. The r ’th bin in the figure counts the percentage of the motion blurry images in the dataset achieving error ratio below r [5]. For example, the bar corresponding to bin 3 indicates the percentage of test images with the SSD error ratio below 3. For each bin, the higher the bar the better the deblurring performance. As pointed out by Levin et al. [6], the finally deblurred images are usually visually plausible when their error ratios are below 3, and in this case the blind motion deblurring is considered to be successful. Obviously, our method has achieved the highest percentage of success among the five methods, i.e., 94 %. As for

123

other four methods, the percentages of success are respectively given as follows: 69 % [1], 75 % [3], 19 % [5], 84 % [6]. Besides, our method also achieves the lowest average SSD deconvolution error ratio, i.e., 1.5392, and the average values of other methods are: 3.7443 [1], 2.4177 [3], 7.0478 [5], 2.0201 [6]. It is particularly observed that, our method has performed far better than the Gaussian image prior-based blind deblurring method [5], demonstrating the necessity and superiority of achieving non-stationarity of the precision parameters in the prior (6). As mentioned earlier, although the prior (6) is expressed in a non-sparse Gaussian form, it is capable of promoting sparsity and in this sense, our proposed method has utilized sparse priors essentially in an implicit way. To demonstrate this statement, the estimated {fd }d∈ , i.e., salient edge images, and the finally deblurred images corresponding to the proposed and other four blur kernel estimation methods are shown in Fig. 3 (We just take Image01-kernel05 in the dataset for example; the blur kernels are magnified just for better visual perception). Besides, the peak signal-to-noise ratio (PSNR) is calculated based on the SSD to measure the quality of deblurred images. From the images (h) and (i) in Fig. 3, it can be seen that our proposed approach has succeeded in pursing the salient edges in the original image, which are utilized as important clues for motion blur kernel estimation. It is also observed that, as for this test case, the PSNR of the deblurred image (Fig. 3g) obtained by our approach is slightly higher than that of the image (Fig. 3b) generated using the ground truth kernel, meaning that a lower SSD error between the estimated blur kernel and the ground truth one does not necessarily imply better blind motion deblurring performance. As advised by one of the reviewers, the average SSD errors of motion blur kernels for different blur kernel estimation methods are provided as follows: 0.0086 [1], 0.0094 [3], 0.0132

J Math Imaging Vis

Fig. 3 Blind motion deblurring for Image01-kernel05 in benchmark dataset [5] (The blur kernels are magnified for better visual perception). a Blurry image, b–g Finally deblurred images using the ground truth motion blur kernel and the motion blur kernels estimated respec-

tively by Fergus et al. [1], Cho and Lee [3], Levin et al. [5], Levin et al. [6] and our proposed approach, h, i Salient edge images in the vertical and horizontal directions estimated by our proposed approach along with the motion blur kernel in (g)

[5], 0.0077 [6], 0.0069 (Ours). It is seen clearly that our proposed approach has achieved the highest accuracy of blur kernel estimation compared with other four methods [1,3,5,6]. In spite of that, the provided average SSD deconvolution error ratios show that the blind deblurring performance of the proposed method is not as good as indicated by the average SSD kernel errors. Since what we are really interested in practice are the finally deblurred images, it is more reasonable to compare different blur kernel estimation methods using the

finally deblurred images, just as Levin et al. [5,6] as well as many other methods do. To sum up, though the proposed method builds on the simple NSG image prior, it has achieved better blind deblurring performance than the influential MoG image prior-based VB blind deblurring methods in [1] and [6], the implicit sparse image prior-based MAP blind deblurring method in [3], as well as the non-sparse Gaussian image priorbased VB blind deblurring method in [5]. In the mean-

123

J Math Imaging Vis

Fig. 4 Ground truth images. a Cameraman, b House, c Lena, d Boat

while, it is seen that [6] performs better than the methods in [1,3,5].

Fig. 5 Synthetic blind motion deblurring with Cameraman. a Blurry Cameraman and ground truth motion blur kernel (size 13 × 13), b–d Finally deblurred images and estimated motion blur kernels by [6] (PSNR: 28.4306 dB, SSIM: 0.9098), [8] (PSNR: 28.5471 dB, SSIM: 0.9171) and our proposed method (PSNR: 29.1161 dB, SSIM: 0.9158)

4.2 Experiments on Standard Test Images in Image Processing In this subsection, the second group of synthetic experiments are performed to test the performance of the proposed method using four standard test images in the community of image processing, including Cameraman, House, Lena, and Boat as shown in Fig. 4, with size 256 × 256 and image values scaled to [0, 255]. The proposed method is compared against the classical Richardson–Lucy based ML blind deblurring method in [25], the MoG prior-based VB blind deblurring method in [6], and the normalized sparsity-based MAP blind deblurring method in [8]. Here, PSNR and SSIM (Structural SIMilarity) are used to quantitatively measure the deblurring performance of different methods. We refer readers to [21] for more details on SSIM. The achieved PSNR and SSIM by different methods are listed in Table 1, from which it is obviously observed that our proposed approach has performed the best among the four methods, while the ML approach [25] has completely failed in all the experiments with the worst performance. We take Fig. 4c for example. The values of PSNR and SSIM achieved by the four methods are respectively: 34.2115 dB, 0.9380 (the proposed approach), 31.3433 dB, 0.9126 (the MoG priorbased approach [6]), 25.9211 dB, 0.8118 (the normalized

123

Fig. 6 Synthetic blind motion deblurring with House. a Blurry House and ground truth motion blur kernel (size 21×21), b–d Finally deblurred images and estimated motion blur kernels by [6] (PSNR: 29.1059 dB, SSIM: 0.8733), [8] (PSNR: 29.9817 dB, SSIM: 0.8834) and our proposed method (PSNR: 36.0149 dB, SSIM: 0.9372)

J Math Imaging Vis

Fig. 7 Synthetic blind motion deblurring with Lena. a Blurry Lena and ground truth motion blur kernel (size 19 × 19), b–d Finally deblurred images and estimated motion blur kernels by [6] (PSNR: 31.3433 dB, SSIM: 0.9126), [8] (PSNR: 25.9211 dB, SSIM: 0.8118) and our proposed method (PSNR: 34.2115 dB, SSIM: 0.9380)

Fig. 8 Synthetic blind motion deblurring with Boat. a Blurry Boat and ground truth motion blur kernel (size 15 × 15), b–d Finally deblurred images and estimated motion blur kernels by [6] (PSNR: 29.3963 dB, SSIM: 0.8744), [8] (PSNR: 26.8675 dB, SSIM: 0.8161) and our proposed method (PSNR: 29.9120 dB, SSIM: 0.8626)

sparsity-based approach [8]), and 16.0381 dB, 0.3625 (the ML approach [25]). The advantage of the proposed approach as well as [6,8] over the ML approach [25] demonstrates well the necessity of imposing some constraints on the inversion of (5) to overcome the ill-posedness of the motion blur kernel estimation problem. For motion deblurring problems corresponding to the four images in Fig. 4, the deblurred images by our proposed method, the MoG prior-based method [6] and the normalized sparsity-based method [8], as well as the corresponding blurry images, the ground truth and estimated motion blur kernels are shown, respectively, in Figs. 5, 6, 7, and 8. It is clearly seen that, the deblurred images by our proposed method are of much better visual perception than [8]; particularly, there are many ringing artifacts along salient edges in the deblurred images by [8] for the images Lena and Boat, while the proposed method performs much better. In fact, that is just why our proposed approach has achieved higher PSNR and SSIM than [8]. Although deblurred images by the proposed method and [6] do not differ much visually, it is seen from Table 1 that [6] achieves lower PSNR and SSIM than our proposed method in nearly all the experiments. In addition, to analyze the proposed EM-based method empirically in more details, the evolution of the parameters  { d }d∈ , {μd }d∈ , {γ d }d∈ , β, k in Algorithm 1 should be studied, too. It is known from Sect. 3 that, the estimation

accuracy of the motion blur kernel k in our proposed approach is highly relevant to that of the parameters {γ d }d∈ and β as well as incorporating the uncertainty of fd into the process of blur kernel estimation. For this reason, the finally estimated precision parameter β in the four experiments is respectively provided, i.e., 0.0114, 0.0108, 0.0113 and 0.0108. As for the evolution of the parameters γ h , γ v and k, we just take Fig. 4d for example. Figures 9, 10 and 11 provide the evolution of the parameters γ h , γ v and k in the finest resolution, i.e., their iterative estimates in T = 10 iterations. It is observed from Figs. 9 and 10 that, our proposed method has succeeded in achieving the non-stationarity of the precision parameters γ h and γ v during the EM iteration, therefore making the NSG prior (6) capable of promoting sparsity. It is also seen that, the proposed approach converges well since the values of γ h and γ v change little in the 10 iterations of the finest scale. In Fig. 11, the iterative estimates of the motion blur kernel k demonstrate more directly that the proposed approach converges to reasonable estimates in practice, showing the effectiveness and robustness of the EM procedure in our NSG-based blind motion deblurring. 4.3 Experiments on Real-World Motion Blurry Images In this subsection, a group of experimental results on realworld motion blurry images are provided in Figs. 12, 13,

123

J Math Imaging Vis Table 1 Achieved PSNR (dB) and SSIM by different methods for synthetic blind motion deblurring using four standard test images in Fig. 4 Image

PSNR

SSIM

[25]

[6]

[8]

Ours

[25]

[6]

[8]

Ours

16.1423

28.4306

28.5471

29.1161

0.4566

0.9098

0.9171

0.9158

House

17.1045

29.1059

29.9817

36.0149

0.3599

0.8733

0.8834

0.9372

Lena

16.0381

31.3433

25.9211

34.2115

0.3625

0.9126

0.8118

0.9380

Boat

17.1751

29.3963

26.8675

29.9120

0.3632

0.8744

0.8161

0.8626

Cameraman

The significance of bold values indicate the largest values of PSNR and SSIM achieved by the methods compared in each motion deblurring problem

Fig. 9 The iterative estimates (10 iterations) of the precision parameters γ h in the finest resolution as blind motion deblurring of Fig. 8a using Algorithm 1

Fig. 10 The iterative estimates (10 iterations) of the precision parameters γ v in the finest resolution as blind motion deblurring of Fig. 8a using Algorithm 1

Fig. 11 The iterative estimates (10 iterations) of the motion blur kernel k in the finest resolution as blind motion deblurring of Fig. 8a using Algorithm 1

123

J Math Imaging Vis

Fig. 12 Realistic blind motion deblurring with Board. a Blurry Board, b–d Finally deblurred images by [6,8] and our proposed method

Fig. 13 Realistic blind motion deblurring with Girl. a Blurry Girl, b–d Finally deblurred images by [6,8] and our proposed method

and 14. As for the motion blurry Board image, it is observed from Fig. 12 that our proposed approach and the MoG priorbased method [6] produce somewhat more faithful estimated

Fig. 14 Realistic blind motion deblurring with Flower. a Blurry Flower, b–d Deblurred images by [6,8] and our proposed method

motion kernels with size 19×19 than [8], and achieve slightly better deblurring performance than [8]. As for the blurry Girl image in Fig. 13, although finally deblurred images by the proposed approach, [6] and [8] are of visually similar perception, the results differ actually more than they appear, since the estimated motion blur kernel with size 21 × 21 by the proposed approach is closer to the truth compared against that of [6] and [8]. As for the blurry Flower image in Fig. 14, it is obvious that [8] performs much worse than our proposed approach and the MoG prior-based method [6], with a lot of ringing artifacts and hence visually unpleasant. The reason is, the estimated motion blur kernel of [8] is not as good as those of other two methods. Till now, intensive experimental results on both synthetic and real-world motion blurry images have demonstrated the rationality and effectiveness of the proposed method, as well as its advantage over the classical ML approach in [25], the MoG image prior-based VB methods in [1,6], the implicit sparse image prior-based MAP method in [3], the non-sparse Gaussian image prior-based Bayesian method in [5] and the normalized sparsity-based MAP method in [8]. In spite of that, the proposed approach is not as efficient as the MAP methods in [3] and [8]. For an image of size 256 × 256 and a motion blur kernel of size 15 × 15, our method generally takes about 60 s, while [3] (C++) takes about 1 s and [8] (MATLAB) takes about 30 s. Another point to be noted is that, there are many parameters to be manually tuned or provided in [3] and [8] for visually acceptable deblurring results,

123

J Math Imaging Vis

while the model parameters in our approach are estimated automatically by the EM algorithm. 5 Conclusions In this paper, a simple yet effective motion blur kernel estimation method is proposed for blind motion deblurring. The method is regularized by a type of non-stationary Gaussian prior placed on the sparse gradient fields of natural images, so as to automatically detect and purse the salient edges in images as important clues to blur kernel estimation. By exploiting the EM algorithm, motion blur kernels, sparse gradient fields and precision parameters involved in the image prior as well as the observation model are automatically estimated in an iterative way. Numerous experimental results have demonstrated the rationality and effectiveness of our proposed approach and its superiority to state-of-the-art methods. Acknowledgments Many thanks are given to the anonymous reviewers for their serious, pertinent and helpful comments significantly strengthening this paper. Wen-Ze Shao is grateful to Professor Michael Elad for the financial support allowing him to work at Department of Computer Science, Technion-Israel Institute of Technology, as well as Professor Yi-Zhong Ma and Dr. Min Wu for their kind supports in the past years. He would also like to show many thanks to Mr. Ya-Tao Zhang and other kind people for helping him through his lost and sad years. The work is supported in part by the National Natural Science Foundation (NSF) of China for Youth under Grant No. 61402239, the NSF of Jiangsu Province for Youth under Grant No. BK20130868, the NSF of Jiangsu Higher Education Institutions under Grant No. 13KJB510022 and 13KJB120005, and the NSF of Nanjing University of Posts and Telecommunications under Grant No. NY212014, NY213007.

References 1. Fergus, R., Singh, B., Hertzmann, A., Roweis, S.T., Freeman, W.T.: Removing camera shake from a single photograph. ACM Trans. Graph 25(3), 787–794 (2006) 2. Shan, Q., Jia, J., Agarwala, A.: High-Quality Motion Deblurring from a Single Image. SIGGRAPH, pp. 1–10. ACM, New York (2008) 3. Cho, S., Lee, S.: Fast motion deblurring. ACM Trans. Graph. (SIGGRAPH ASIA 2009), 28(5), article no. 145 (2009) 4. Almeida, M., Almeida, L.: Blind and semi-blind deblurring of natural images. IEEE Trans. Image Process. 19, 36–52 (2010) 5. Levin, A., Weiss, Y., Durand, F., Freeman, W.T.: Understanding and evaluating blind deconvolution algorithms. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1964–1971 (2009) 6. Levin, A., Weiss, Y., Durand, F., Freeman, W.T.: Efficient marginal likelihood optimization in blind deconvolution. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 2657–2664 (2011) 7. Amizic, B., Molina, R., Katsaggelos, A.K.: Sparse Bayesian blind image deconvolution with parameter estimation. EURASIP J. Image Video Process. 2012(20), 1–15 (2012)

123

8. Krishnan, D., Tay, T., Fergus, R.: Blind deconvolution using a normalized sparsity measure. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 233–240 (2011) 9. Babacan, S.D., Molina, R., Do, M.N., Katsaggelos, A.K.: Bayesian blind deconvolution with general sparse image priors. In: Fitzgibbon A. et al. (eds.) ECCV 2012, Part VI, LNCS, vol. 7577, pp. 341–355 10. Miskin, J.W., MacKay, D.J.C.: Advances in Independent Component Analysis. Ensemble learning for blind image separation and deconvolution. Springer, Berlin (2000) 11. Xu, L., Jia, J.: Two-phase kernel estimation for robust motion deblurring. In: Daniilidis, K., Maragos, P., Paragios, N. (eds.): ECCV 2010, Part I, LNCS, vol. 6311, pp. 157–170 12. Joshi, N., Szeliski, R., Kriegman, D.J.: PSF estimation using sharp edge prediction. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1–8 (2008) 13. Joshi, N.: Enhancing photographs using content-specific image priors. PhD thesis, University of California, San Diego (2008) 14. Cai, J.F., Ji, H., Liu, C., Shen, Z.: Blind motion deblurring from a single image using sparse approximation. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 104–111 (2009) 15. Vera, E., Vega, M., Molina, R., Katsaggelos A.K.: A novel iterative image restoration algorithm using nonstationary image priors. In: The 18th IEEE International Conference on Image Processing (ICIP), pp. 3457–3460 (2011) 16. Tzikas, D., Likas, A., Galatsanos, N.P.: Variational Bayesian sparse kernel-based blind image deconvolution with Student’s-t priors. IEEE Trans. Image Process. 18(4), 753–764 (2009) 17. Tipping, M.E.: Sparse Bayesian learning and the relevance vector machine. J. Mach. Learn. Res. 1, 211–244 (2001) 18. Huang, Y., Ng, M.K., Wen, Y.: A fast total variation minimization method for image restoration. SIAM J. Multiscale Model. Simul. 7(2), 774–795 (2008) 19. Shao, W.-Z., Deng, H.-S., Wei, Z.-H.: Multi-Parseval frame-based nonconvex sparse image deconvolution. Opt. Eng. 51(6), 0670081–067008-12 (2012) 20. Krishnan, D., Fergus, R.: Fast image deconvolution using hyperLaplacian priors. Adv. Neural Inf. Process. Syst. 22, 1033–1041 (2009) 21. Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P.: Image quality assessment: from error measurement to structural similarity. IEEE Trans. Image Process. 13(4), 600–612 (2004) 22. Levin, A., Fergus, R., Durand, F., Freeman, W.T.: Image and depth from a conventional camera with a coded aperture. ACM Trans. Graph. 26(3), 1–9 (2007). Article 70 23. Lucy, L.: An iterative technique for the rectification of observed distributions. Astron. J. 79, 745–754 (1974) 24. Richardson, W.: Bayesian-based iterative method of image restoration. J. Opt. Soc. Am. 62(1), 55–59 (1972) 25. Fish, D.A., Brinicombe, A.M., Pike, E.R., Walker, J.G.: Blind deconvolution by means of the Richardson–Lucy algorithm. J. Opt. Soc. Am. A 12, 58–65 (1995) 26. Xu, L., Zheng, S., Jia, J. : Unnatural L 0 sparse representation for natural image deblurring. In: IEEE Conference Computer Vision and Pattern Recognition (CVPR), pp. 1107–1114 (2013) 27. Wipf, D., Zhang, H.: Revisiting Bayesian blind deconvolution. arXiv: 1305.2362 (2013) 28. Dempster, A., Laird, N., Rubin, D.: Maximum likelihood from incomplete data via the EM algorithm. J. R Stat. Soc. A 39(1), 1–38 (1977)

J Math Imaging Vis Wen-Ze Shao was born in Ganyu, Jiangsu Province, China. He received the B.S. degree in Science of Information and Computation in June 2003, and the Ph.D. degree in Pattern Recognition and Intelligent System in July 2008, both from Nanjing University of Science and Technology (NUST), Nanjing, China. From June 2003 to December 2011, he served as an officer in the People’s Liberation Army (PLA) of China. In January 2012, he joined Nanjing University of Posts and Telecommunications (NUPT) as an assistant professor. Since May 2014, he has been also working as an academic visitor (Postdoc) at the Department of Computer Science in Israel Institute of Technology (Technion). He is particularly interested in the fields of natural image modeling, synthesis and analysis sparse modeling, statistical machine learning, variational PDE and Bayesian methods with applications to innovative and interactive imaging and vision problems. Qi Ge was born in Nantong, Jiangsu Province, China. She received the B.S. degree in Science of Information and Computation in 2006 and the M.S. degree in Applied Mathematics in 2009, both from Nanjing University of Information and Engineering, Nanjing, China, and the Ph.D. degree in Pattern Recognition and Intelligent System in 2013 from NUST, Nanjing, China. Now, she serves as an assistant professor at College of Telecommunications and Information Engineering at NUPT. She is particularly interested in the fields of variational PDE approaches, probabilistic graphical models, sparse representation and their applications to image restoration, segmentation, and so on.

Zhi-Hui Wei was born in Huai’an, Jiangsu Province, China. He received the B.S. degree in 1983 and the M.S. degree in 1986, both in Mathematics, and the Ph.D. degree in Electrical Engineering in 2003, all from Southeast University, Nanjing, China. Now, he is a Professor of Image Processing and Pattern Recognition, and Mathematics in NUST, Nanjing, China. He has been working on the fields of statistical image modeling, geometrical multi-scale analysis, variational PDE, sparse representation, compressed sensing, superresolution, image compression, and so on. Hai-Bo Li was born in Xuzhou, Jiangsu Province, China. He received the B.E. degree in Wireless Engineering in 1985 and the M.S. degree in Communication and Electronic Systems in 1988, both from NUPT, Nanjing, China, and the Ph.D. degree in Information Theory in 1993 from Lin ping University, Sweden. In 1997, he became docent (UK equivalent senior lecturer, US equivalent associate professor) in image coding, and in 1999 he became the youngest lifetime professor of signal processing in Umea University, Sweden. Now, he is a Professor of innovative media technology in KTH Royal Institute of Technology. His research interest is mainly media signal processing, including image coding, video compression, motion estimation, facial and hand gesture recognition for the next generation of mobile phones, invisible interaction technology, and so on.

Hai-Song Deng was born in Xuzhou, Jiangsu Province, China. She received the B.S. degree in Mathematics in 2003 from Xuzhou Normal University, Xuzhou, China, the M.S. degree in Statistics in 2006, and the Ph.D. degree in Management Science and Engineering in 2011 both from NUST, Nanjing, China. Since 2011, she served as an assistant professor at School of Science in Nanjing Audit University. Her research fields include Bayesian variable selection, sparse optimization, surrogate modeling, design and analysis of computer experiments, and so on.

123