PSF Estimation via Gradient Cepstrum Analysis for Image Deblurring in Hybrid Sensor Network

Abstract

In hybrid sensor networks, information fusion from heterogeneous sensors is important, but quite often information such as image is blurred. Single image deblurring is a highly ill-posed problem and usually regularized by alternating estimating point spread function (PSF) and recovering blur image, which leads to high complexity and low efficiency. In this paper, we first propose an efficient PSF estimation algorithm based on gradient cepstrum analysis (GCA). Then, to verify the accuracy of the strategy, estimated PSFs are used for image deconvolution step, which exploits a novel total variation model coupling with a gradient fidelity term. We also adopt an alternating direction method (ADM) numerical algorithm with rapid convergence and high robustness to optimize the energy function. Both synthetic and real blur experiments show that our scheme can estimate PSF rapidly and produce comparable results without involving long time consuming.

1. Introduction

In real world, information fusion in hybrid sensor network (HSN) is necessary in different applications [1–3]. For example, in an emergency natural disaster scenario, HSN-based multimodal information integration for first responders is critical for search and rescue. Danger may appear anywhere at any time; therefore, first responders must monitor a large area continuously in order to identify potential danger and take actions. Due to the dynamic and complex nature of natural disaster, some victims may not be found with a single type of sensor modality; for example, image sensors can be used to spot the victims based on optical images; UWB radar sensors can be applied to penetrate the ground or sense-through-wall [4, 5], and acoustic sensors are needed to collect the voice from victims [6]. Overall, image sensors play important roles in HSN. Unfortunately, single has long been a fundamental research problem in many science and engineering areas such as computer vision, aerial imaging, and remote sensing, especially in many photographs capturing ephemeral moments that cannot be recaptured under some certain conditions. Assuming that the imaging system is linear space-invariant, the image degraded process can be modeled as convolution of a latent image with a blur kernel. The progress is described as

$$\begin{array}{cc}& g\left(x,y\right)=f\left(x,y\right)\otimes h\left(x,y\right)+n\left(x,y\right),\end{array}$$

(1)

where $(x,y)$ denotes the pixel coordinate in spatial domain. $g(x,y)$ , $f(x,y)$ , and $n(x,y)$ represent the observed blur image, the latent image, and the additive noise, respectively. ⊗ denotes the convolution operator, and $h(x,y)$ is the blur kernel, which is often described as a point spread function (PSF). Because both latent image and the blur kernel are unknown, single image deblurring reduces to blind image deconvolution, which is significantly ill-posed. Therefore, PSF estimation places an important part in resolving this challenging problem.

Early approaches typically assume simple parametric forms for PSF, such as a Gaussian filter in the frequency domain [7], linear or harmonic motion [8], or a sum of normal motion [9]. In practical cases, blur kernels are often more complex and should be described by more sophisticated parametric models. Cai et al. [10] present an approach that formulates the blind blurring as a new joint optimization problem, which simultaneously maximizes the sparsity of the blur kernel and the sparsity of the clear image under certain suitable redundant tight frame systems. Krishnan et al. [11] use a normalized sparsity measure scheme that compensates for the attenuation of high frequencies, which greatly stabilizes the PSF estimation process. Shi et al. [12] propose a novel local prior that is called gradient fidelity constraint on the global prior of heavy-tailed distribution for single image deblurring.

Recent state-of-the-art PSF estimation algorithms tend to deploy effective priors about the statistics of natural images for single image deblurring. Fergus et al. [13] use a zero-mean mixture-of-Gaussian model for natural image gradients distribution together with variational Bayes algorithm to estimate PSF, and then Richardson-Lucy is employed for deconvolution. The proposed method provides good results just for the case of small PSF. Krishan and Fergus [14] model the heavy-tailed distribution by hyper-Laplacian priors and adopt an alternating minimization scheme, which is fast optimized by a lookup table (LUT) algorithm. However, there are significant staircase effects in the recovered image. Levin et al. [15] also propose an efficient marginal likelihood optimization method based on the maximum a posteriori (MAP) that is $\mathrm{M}\mathrm{A}{\mathrm{P}}_h$ principle and proved more robust than the common $\mathrm{M}\mathrm{A}{\mathrm{P}}_{f,h}$ . Yet it involves a challenging marginalization over latent images and computational complexity.

Although above algorithms obtain good results, they are time consuming. Because optimizing the energy function usually requires quite sophisticated iterative numerical algorithms, especially the optimizing scheme that alternates between PSF estimation and image restoration. Moreover, if the initialized kernel is not well set in the matched style or proper size, it often fails to converge to the true global minimum, such as the well-known case that the kernel converges to a delta function but the recover image is still blurred. In this paper, we propose a novel efficient gradient cepstrum analysis (GCA) strategy into the single blurred image for PSF estimation. For a given blurred image, we first estimate its gradient cepstrum, which is proved to represent the property of the blur kernel and provide the magnitude information. The phase retrieval (PR) technique is adopted to recover the phase given just the magnitude of a two-dimension PSF Fourier transform. Then, the total variation regularized model coupling with an image gradient fidelity term is established to evaluate the accuracy of our PSF estimation strategy and an alternating direction method (ADM) with rapid and stable convergence is used to optimize the energy function in order not to cover up the efficiency of the proposed GCA strategy. Finally, both synthetic and real blurred images are tested to verify the performance of our method.

The rest of this paper is organized as follows. In Section 2, we make a description of the proposed GCA strategy. Section 3 describes our novel TV based image restoration model and the ADM numerical algorithm. Experimental results are shown in Section 4. Finally, a conclusion is made in Section 5.

2. PSF Estimation by Gradient Cepstrum Analysis Strategy

A block diagram of the degrade model and our proposed single image deblurring scheme including the PSF estimation and the fast image deconvolution is shown in Figure 1.

Figure 1

Block diagram of the degrade model and the proposed single image deblurring process.

In this section, we focus on the PSF estimation process and the details are described in the following subsections.

2.1. Definition and Property of the Image Cepstrum

The linear space-invariant degrade process is described as (1). In the case of ignoring the additive noise, it is formulated in frequency domain as

$$\begin{array}{cc}& G\left(u,v\right)=F\left(u,v\right)H\left(u,v\right),\end{array}$$

(2)

where $(u,v)$ represents the index in frequency domain, $G(u,v)$ , $F(u,v)$ , and $H(u,v)$ denote the discrete Fourier transform of the blurred image, PSF, and latent image, respectively. In image restoration technology, cepstrum of the image $g(x,y)$ is defined as

$$\begin{array}{cc}& C_g\left(p,q\right)=\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}\left[\left|\log \left|G\left(u,v\right)\right|\right|\right],\end{array}$$

(3)

where $(p,q)$ and $C_g(p,q)$ indicate the index in cepstrum domain and the cepstrum of image $g(x,y)$ , respectively. $\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}$ denotes the inverse Fourier transform. In practical applications, to avoid the ill-posed problem caused by $\left|G(u,v)\right|=\mathrm{0}$ , $C_g(p,q)$ is often defined as

$$\begin{array}{cc}& C_g\left(p,q\right)=\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}\left\{\log \left[\mathrm{1}+\left|G\left(u,v\right)\right|\right]\right\}.\end{array}$$

(4)

The cepstrum of degrade progress is described as

$$\begin{array}{cc}& C_g\left(p,q\right)=C_f\left(p,q\right)+C_h\left(p,q\right),\end{array}$$

(5)

where $C_f(p,q)$ and $C_h(p,q)$ represent the cepstrum of the latent image and the PSF. Obviously, the convolution in spatial domain is translated into the additive form, which is easy to separate the information of the PSF. Here the autocorrelation function of the image is associated and described as

$$\begin{array}{cc}& R_g\left(p,q\right)=\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}\left[\left|G\left(u,v\right)\right|\right].\end{array}$$

(6)

The difference between the cepstrum and the autocorrelation function is the weighted logarithmic, which is adopted to concentrate the energy of the transformed signal and expand the range of dynamic analysis and improve the accuracy of the transform. Essentially, image cepstrum has same properties as its autocorrelation, as shown in Figure 2. Figure 2(a) is the original image of a natural scene and its isotropic autocorrelation energy spectrum is shown in Figure 2(b). Figure 2(c) is the autocorrelation energy spectrum of motion blurred image; motion length and angle are 20 and 30, respectively. It is obvious that motion blurred image has better correlation in its blur direction.

Figure 2

Natural image and its autocorrelation energy spectrum. (a) Original image. (b) Autocorrelation energy spectrum of (a). (c) Autocorrelation energy spectrum of blurred image.

The cepstrum has many important properties. We only discuss the two-dimensional cepstrum properties related to the field of image processing.

Property 1.

Convert the two-dimensional convolution into an addition operation in cepstrum domain.

Property 2 (isogonal rotation).

That is, in polar coordinates system, if

$$\begin{array}{cc}& g\left(r,\theta \right)=C_g\left(r,\gamma \right),\end{array}$$

(7)

one can get

$$\begin{array}{cc}& g\left(r,\theta -\alpha \right)=C_g\left(r,\gamma -\alpha \right).\end{array}$$

(8)

It indicates that if the original function rotates one angle α, its cepstrum rotates with the same direction and degree.

Property 3 (origin symmetry).

If $g(x,y)$ is a real function, then its cepstrum $C_g(p,q)$ is symmetric about the origin. That is,

$$\begin{array}{cc}& C_g\left(p,q\right)=C_g\left(-p,-q\right).\end{array}$$

(9)

Property 4.

Periodicity

$$\begin{array}{cc}& C_g\left(p+M,q+N\right)=C_g\left(p,q\right),\end{array}$$

(10)

where M, N represent the image size in p and q axis, respectively.

These two-dimensional cepstrum properties provide a reliable theoretical basis of our algorithm.

2.2. PSF Estimated from the Blurred Image Gradient Cepstrum

PSF, that is, the blur kernel of the convolution operation, actually describes the path of camera shake. As the camera shake path is random, the corresponding two-dimensional image of shake path is sparse and random. Fortunately, cepstrum is not sensitive to the effects of the random transfer path, but also without being affected by other random components in imaging progress. Take the random noise, for example; it will produce some frequency components, which are not periodic in frequency domain. There are not obvious peaks corresponding to the random noise in cepstrum domain. Therefore, PSF information can be extracted by cepstrum analysis of the blurred image.

In PSF estimation, image gradient provides more useful information than the image itself. Compared with the cepstrum of the blurred image itself, its gradient cepstrum has the advantage of isolating the information of the image itself effectively and estimating the PSF accurately. Considering that the Laplace operator is isotropic and rotation-invariant, we choose it as gradient operator. It can reflect second-order differential properties of the image, thus extracting points, lines, and boundaries for the image. It is also known as a boundary extraction operator and defined as

$$\begin{array}{cc}& {\nabla }^{\mathrm{2}}=\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}.\end{array}$$

(11)

There is

$$\begin{array}{cc}& {\nabla }^{\mathrm{2}}g\left(x,y\right)={\nabla }^{\mathrm{2}}f\left(x,y\right)\otimes h\left(x,y\right),\end{array}$$

(12)

and then

$$\begin{array}{cc}& C_{{\nabla }^{\mathrm{2}}g}=C_{{\nabla }^{\mathrm{2}}f}+C_h.\end{array}$$

(13)

In this paper, $C_{{\nabla }^{\mathrm{2}}f}$ is proved to be approximate to a delta function in the following experiment; then

$$\begin{array}{cc}& C_{{\nabla }^{\mathrm{2}}g}\approx C_h.\end{array}$$

(14)

It means that we can estimate the PSF by analyzing the cepstrum of the gradient of the blur image, just as the above mentioned strategy, gradient cepstrum analysis (GCA). Here, we choose ${\nabla }^{\mathrm{2}}$ to be a very compact gradient operator,

$$\begin{array}{cc}& {\nabla }^{\mathrm{2}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{1}& \mathrm{1}\\ \mathrm{1}& -\mathrm{8}& \mathrm{1}\\ \mathrm{1}& \mathrm{1}& \mathrm{1}\end{array}\right].\end{array}$$

(15)

To illustrate the property more clearly, we carry out experiments on the synthetic blurred image; the result is shown in Figure 3. In Figure 3(a), the synthetic blurred image is degraded by a motion blur PSF with the blur length 20 pixels and the blur angle 30 degrees. Figures 3(b) and 3(c) show the cepstrum of the blur image and the original image, respectively. It can be seen that the cepstrum of the blurred image has rich information in the blur direction. This phenomenon matches the property of image autocorrelation that pixels of the blur image are more correlated in the direction of the blurring trajectory of the PSF than other directions [16]. The cepstrum of the original image gradient is shown in Figure 3(d). It is approximate to a delta function and matches the statistical property that gradients of the natural images are approximately independent to each other [13]. Figure 3(e) is the cepstrum of the gradient of the blur image. Obviously, it isolates the cepstral information of the original image effectively and preserves the information in the blurred direction. Figure 3(f) is the cepstrum of the PSF itself. It is centrosymmetric. The direction remains the same and the length is twice than the true PSF. The reproduced PSF information conforms to the physical meaning of the cepstrum; that is, the cepstrum of the signal can reproduce the original signal information. It is evident that there are great similarities between the cepstrum of the gradient of the blur image and the cepstrum of the PSF itself, as shown in Figures 3(e) and 3(f). Therefore, the cepstrum of the gradient of the blur image is approximate to the cepstrum of the PSF itself under the case that $C_{{\nabla }^{\mathrm{2}}f}$ is treated as a delta function and its influences on the PSF estimation will be ignored. The upper right or low left portion of the blur image gradient cepstrum is chosen as the magnitude of the estimated PSF. The phase information is recovered by the following phase retrieval technology.

Figure 3

(a) The synthetic blurred image. (b) The cepstrum of (a). (c) The cepstrum of the original image. (d) The cepstrum of the original image gradient. (e) The cepstrum of the gradient of blur image (a). (f) The cepstrum of the PSF itself.

2.3. Phase Retrieval

Phase spectrum of the image plays an important role in reflecting the image composition contents even more important than the power spectrum. The estimated PSF from the cepstrum of the blur image gradient only remain the amplitude information but the phase information is lost. The PR [17] technique is adopted here to derive the phase of the PSF. It recovers the phase just from the magnitude of a signal's Fourier transform, which is the iterative Fourier transform algorithm. The iterative process includes three constraints, that is, positivity, compact support, and module constraints. To be exact, positivity and compact support constraints are spatial domain constraints and the module belongs to the frequency domain constraints. The flow chart for the algorithm is shown in Figure 4.

Figure 4

Flow chart of the iterative Fourier transform algorithm.

The PR algorithm can be expressed as follows.

•  

Initialize ${\phi }_{\mathrm{0}}(u,v)=\mathrm{0}$ and $H_{\mathrm{0}}(u,v)=\left|H(u,v)\right|e^{i·\mathrm{0}}$ .

 

While iteration $k<\mathrm{100}$ , Do

• (1)

Apply the module constraint in Fourier domain $H_k^{\prime }(u,v)=\left|H_{\mathrm{0}}(u,v)\right|e^{i{\phi }_k(u,v)}$ .

(2)

Transform to spatial domain. $h_k^{\prime }(p,q)=\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}\left[H_k^{\prime }(u,v)\right]$

 

Apply the compact constraint and get $h_{k+\mathrm{1}}(p,q)=\left\{\begin{array}{cc}{h}_k^{\prime }(p,q)& (p,q)\in \eta \\ 0& (p,q)\notin \eta \end{array}\right.$

 

Apply the positivity constraint $h_{k+\mathrm{1}}(p,q)>\mathrm{0}$ .

(3)

Transform to frequency domain $H_{k+\mathrm{1}}(u,v)=\mathrm{F}\mathrm{F}\mathrm{T}\left[h_{k+\mathrm{1}}(p,q)\right]$ .

 

End Do

 

Finally, $H_{\mathrm{e}\mathrm{s}\mathrm{t}}(u,v)=\left|H_{\mathrm{0}}(u,v)\right|e^{j{\phi }_{k+\mathrm{1}}(u,v)}$ and $h_{\mathrm{e}\mathrm{s}\mathrm{t}}=\left|\mathrm{F}\mathrm{F}{\mathrm{T}}^{-\mathrm{1}}\left[H_{\mathrm{e}\mathrm{s}\mathrm{t}}(u,v)\right]\right|$ .

Where $H_k^{\prime }(u,v)$ and $h_k^{\prime }(p,q)$ are temporary variables and do not represent usual matrix operations. η denotes the set of points that meet the spatial domain constraint. $H_{\mathrm{e}\mathrm{s}\mathrm{t}}(u,v)$ is the estimated frequency result at the end of the iteration and $h_{\mathrm{e}\mathrm{s}\mathrm{t}}$ represents the estimated PSF. In general case, iterative Fourier transforms for the image PR technique will take a long computing time, but in this paper, the PSF size is much smaller than that of the image and constrained by the module and the compact support. Therefore, the weakness of long time consuming is not exposed in our PSF estimation.

3. The Optimization Approach

In this section, we propose a total variation regularized model coupling with an image gradient fidelity term to evaluate the accuracy of our PSF estimation strategy and an alternating direction method (ADM) with rapid and stable convergence is used to optimize the energy function.

3.1. The Proposed Model

It is well known that the error between the estimated PSF and the true degrade PSF is inevitable no matter how excellent the PSF estimation method is, which will lead to undesirable ringing effects. Total variation image restoration algorithm has become one of the standard techniques with the advantage of suppressing ringing effects and preserving sharp edges and object boundaries [18]. However, it still has a potential shortcoming that does not meet the morphological principle of image processing, which is often referred to the staircase effect in the recovered image [19]. Although the staircase effect is reluctantly acceptable in the image features extraction or target detection, it does not meet the human visual requirement. In order to use the unique advantage of TV regularization and avoid its disadvantage, we propose a novel TV based regularized model coupling with a gradient fidelity term to recover the latent image from the single blurred image.

The model is defined as

$$\begin{array}{cc}& \mathrm{m}\mathrm{i}\mathrm{n}{\Phi }_{\mathrm{r}\mathrm{e}\mathrm{g}}\left(f\right)+\beta {\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}\left(f,g\right)+\gamma {\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}\left(\nabla f,\nabla g\right),\end{array}$$

(16)

where ${\Phi }_{\mathrm{r}\mathrm{e}\mathrm{g}}(f)$ regularizes the solution by TV norm, and

$$\begin{array}{cc}& \mathrm{T}\mathrm{V}\left(f\right)={{\displaystyle \sum }}_i‖D_{i}f‖,\end{array}$$

(17)

i indexes pixels in the image. D is the first-order finite difference operator $D_x=\left[\mathrm{1},-\mathrm{1}\right],D_y={\left[\mathrm{1},-\mathrm{1}\right]}^T$ . ${\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}(f,g)$ called the energy fidelity term is to measure the violation of the relation between f and the observation g. It is constrained by the noise and based on the assumption that noise random variables for all pixels are independent and identically distributed and each of them follows a Gaussian distribution. We define

$$\begin{array}{cc}& {\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}\left(f,g\right)=\frac{{‖f\otimes h-g‖}^{\mathrm{2}}}{\mathrm{2}}.\end{array}$$

(18)

${\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}(\nabla f,\nabla g)$ is the gradient fidelity term that constrains the gradient of the recovered image to be similar with the gradient of the latent image. It plays an important part in maintaining image smoothing and suppressing the staircase effect. Then

$$\begin{array}{cc}& {\Phi }_{\mathrm{f}\mathrm{i}\mathrm{d}}\left(\nabla f,\nabla g\right)=\frac{{‖\nabla f-\nabla g‖}^{\mathrm{2}}}{\mathrm{2}}.\end{array}$$

(19)

β and γ are constraint coefficients chosen as constants. Then we obtain the energy function

$$\begin{array}{cc}& E\left(f\right)={{\displaystyle \sum }}_i\left[‖D_{i}f‖+\frac{\beta }{\mathrm{2}}{‖f\otimes h-g‖}^{\mathrm{2}}+\frac{\gamma }{\mathrm{2}}{‖\nabla f-\nabla g‖}^{\mathrm{2}}\right].\end{array}$$

(20)

3.2. Numerical Algorithm

Optimization of the above energy function refers to the TV problem in $L^2$ norm. There are many numerical methods to solve it. A well-known linearized gradient method [20] is proposed to solve the Euler-Lagrangian equation via fix-point iteration. However, it is not efficient especially for large or severely ill-conditioned problem. Because each iteration is required to solve a linear system that becomes more and more difficult as the kernel becomes more ill-conditioned. Beck and Teboulle [21] present a fast iterative shrinkage-thresholding algorithm that simplifies the computation and has an optimal global rate of convergence. But it is not able to take advantage of problem structures. To further solve numerical difficulties, Goldstein and Osher [22] apply the classical augmented Lagrangian method (ALM), which updates the iterative scheme and guarantees the overall convergence but still requires solving each unconstrained subproblem almost exactly. Tao [23] adopts a variant of the classical augmented Lagrangian method that is the alternating direction method (ADM) for total variation structured optimization, which is proved to be more stable, efficient, and particularly faster than the ALM algorithm. The versatility of the ADM [24] has demonstrated its efficiency and robustness especially for large scale problems that are corresponding to the size solvable by interior-point methods. In this paper, we use the ADM to optimize the energy function and realize fast image deconvolution with the purpose of not covering the fast advantage of our PSF estimation strategy.

We introduce an auxiliary variable $w=(w_{\mathrm{1}},w_{\mathrm{2}}{)}^T$ , where $w_{\mathrm{1}}=D_{x}f$ , $w_{\mathrm{2}}=D_{y}f$ and the iterative parameter λ. Let $\nabla f\approx D_{i}f$ and the augmented Lagrangian function of (20) is expressed as

$$\begin{array}{cc}& L_A\left(f,w,\lambda \right)={{\displaystyle \sum }}_i\left[‖w_i‖-{\lambda }_i^T\left(w_i-D_{i}f\right)+\frac{\alpha }{\mathrm{2}}{‖w_i-D_{i}f‖}^{\mathrm{2}}+\frac{\gamma }{\mathrm{2}}{‖w_i-{\left(\nabla g\right)}_i‖}^{\mathrm{2}}\right]+\frac{\beta }{\mathrm{2}}{‖f\otimes h-g‖}^{\mathrm{2}}.\end{array}$$

(21)

Applying the ADM to optimize $L_A$ , its iterative scheme is

$$\begin{array}{cc}& w^{k+\mathrm{1}}⟵\underset{w}{\arg \hspace{0.17em}\min }{L}_A\left(f^k,w,{\lambda }^k\right)\\ f^{k+\mathrm{1}}⟵\underset{w}{\mathrm{a}\mathrm{r}\mathrm{g}\hspace{0.17em}\mathrm{m}\mathrm{i}\mathrm{n}}{L}_A\left(f,w^{k+\mathrm{1}},{\lambda }^k\right)\\ {\lambda }^{k+\mathrm{1}}⟵{\lambda }^k-\alpha \left(w^{k+\mathrm{1}}-D{f}^{k+\mathrm{1}}\right).\end{array}$$

(22)

We fix $f^k$ and ${\lambda }^k$ and optimize w. Equation (21) is transformed to minimizing

$$\begin{array}{cc}& L_A\left(f^k,w,{\lambda }^k\right)=‖w_i‖+\frac{\alpha }{\mathrm{2}}{‖w_i-\left(D_i{f}^k+\frac{\mathrm{1}}{\alpha }{\left({\lambda }^k\right)}_i\right)‖}^{\mathrm{2}}+\frac{\gamma }{\mathrm{2}}{‖w_i-\nabla g‖}^{\mathrm{2}}.\end{array}$$

(23)

By solving its Euler-Lagrange equation, we obtain the solution

$$\begin{array}{cc}& w_i^{k+\mathrm{1}}=‖D_i{f}^k+\frac{\mathrm{1}}{\alpha }{\left({\lambda }^k\right)}_i+\frac{\gamma }{\alpha }{\left(\nabla g\right)}_i‖-\frac{\mathrm{1}+\gamma }{\alpha }\frac{w_i^k}{‖w_i^k‖},\end{array}$$

(24)

which is given explicitly by the two-dimensional shrinkage

$$\begin{array}{cc}& w_i^{k+\mathrm{1}}=\mathrm{s}\mathrm{h}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{k}\left(D_i{f}^k+\frac{\mathrm{1}}{\alpha }{\left({\lambda }^k\right)}_i+\frac{\gamma }{\alpha }{\left(\nabla g\right)}_i,\frac{\mathrm{1}+\gamma }{\alpha }\right).\end{array}$$

(25)

Then, fix ${\lambda }^k$ and $w^{k+\mathrm{1}}$ and optimize f. Equation (21) is transformed to minimizing

$$\begin{array}{cc}& L_A\left(f,w^{k+\mathrm{1}},{\lambda }^k\right)={{\displaystyle \sum }}_i\left({\lambda }_i^T{D}_{i}f+\frac{\alpha }{\mathrm{2}}{‖w_i-D_{i}f‖}^{\mathrm{2}}\right)+\frac{\beta }{\mathrm{2}}{‖f\otimes h-g‖}^{\mathrm{2}}.\end{array}$$

(26)

The minimization of (26) with respect to f is a least squares problem and its corresponding normal equation is

$$\begin{array}{cc}& \left(D^{T}D+\frac{\beta }{\alpha }{h}^{T}h\right)f=D^T\left(w^{k+\mathrm{1}}-\frac{\mathrm{1}}{\alpha }{\lambda }^k\right)+\frac{\beta }{\alpha }{h}^{T}g.\end{array}$$

(27)

In order to ensure the nonsingularity of the coefficient matrix in (27), it shall obey the standard assumption; that is, $N(h)\cap N(D)=\mathrm{0}$ , where $N(·)$ represents the null space of a matrix. Both $D^{T}D$ and $h^{T}h$ are block circulant matrix when f is under the periodic boundary condition. Therefore, $D^{T}D+(\beta /\alpha )h^{T}h$ can be diagonalizable by the two-dimension discrete Fourier transform and thereby obtain its decomposed eigenvalues, which requires being computed only once. We compute the right-hand side of (27) and apply a forward Fourier transform to get $\mathrm{F}\mathrm{F}\mathrm{T}(f^{k+\mathrm{1}})$ . Then, the two dimensional inverse discrete Fourier transform is applied to get the new iterated $f^{k+\mathrm{1}}$ . Therefore, each iteration only needs two FFTs (including once inverse FFT), which greatly reduces the computing time. The ADM numerical algorithm can be described as follows.

•  

Initialization: $f=g$ , $\lambda ={\lambda }_{\mathrm{0}}$ , $\alpha >\mathrm{0}$ , $\beta >\mathrm{0}$ , $\gamma >\mathrm{0}$ ;

 

Iteration:

• (1)

Compute $w^{k+\mathrm{1}}$ according to (22) for given $(f^k,{\lambda }^k)$ .

(2)

Compute $w^{k+\mathrm{1}}$ via (27).

(3)

Update ${\lambda }^{k+\mathrm{1}}$ via (22).

 

Stop: Until $‖f^{k+\mathrm{1}}-f^k‖/\mathrm{m}\mathrm{a}\mathrm{x}\left\{‖f^k‖,\mathrm{1}\right\}\le \mathrm{t}\mathrm{o}\mathrm{l}$ . Here $\mathrm{t}\mathrm{o}\mathrm{l}$ is set to 10⁻³.

4. Experimental Results

In order to evaluate the performance of our scheme, we carry out a series of experiments on both synthetic and real blurred images and compare with some state-of-the-art algorithms.

Firstly, we carry out the experiment on the synthetic blurred image from [15] and their corresponding PSFs. The size of the tested kernel is 27 × 27 and results are shown in Figure 5. Figures 5(a) and 5(b) are the cepstrum of the PSF itself and the cepstrum of the gradient of blur image. Figure 5(c) is the real PSF and Figure 5(d) is the estimated PSF by our scheme. It can be seen that our estimated PSF is extremely similar to the true PSF.

Figure 5

The GAC of the blur image. (a) The cepstrum of the PSF itself. (b) The cepstrum of the gradient of blur image. (c) The true PSF. (d) Our estimated PSF.

It should be noted that original result might be a mirrored or shifted version of the true PSF, because all of them have the same cepstrum function. However, it is easy to validate and correct since wrong PSFs give very different deblurred result. Experimental results on other tested PSFs are shown in Figure 6. From left to right, their sizes are 13 × 13, 17 × 17, 21 × 21, 23 × 23, and 27 × 27, respectively. Figures 6(a) and 6(b) are a series of true PSFs and responding estimated results from [15]. Figure 6(c) shows estimated PSFs by Hu et al. in [16] and our results are displayed in Figure 6(d). Clearly, our results match well with true PSFs in visual effects.

Figure 6

Estimated PSFs by different algorithms. (a) True PSFs. (b) Estimated PSFs by Fergus et al. [13]. (c) Estimated PSFs by Krishan and Fergus [14]. (d) Our results.

We take the estimated PSFs to image deconvolution and adopt signal to noise (SNR) of restored images to evaluate the accuracy of our proposed PSF estimated method, as shown in the bottom-right of images in Figure 7. Figure 7(a) is a clear image and degraded by convolving with the mentioned 13 × 13 PSF shown in the bottom-right in Figure 7(b). Figures 7(c) and 7(d) are results by Fergus and Levin's methods, respectively. Figure 7(e) provides the result by Hu et al.'s method that uses gradient domain correlation based on a patch-based image degradation model. Figure 7(f) shows our result. Compared with other three state-of-art algorithms, it is evident that our scheme can also get the ideal restoration in visual although the SNR is slightly lower.

Figure 7

Estimated PSFs and corresponding restored images by different algorithms. (a) The clear image. (b) Blurred image and the degrade PSF. (c) The result of Fergus et al.'s method in [13]. (d) The result of Levin et al.'s method in [15]. (e) The result of Hu et al.'s method in [16]. (f) Our result.

We also test other PSFs in Figure 6. SNR between original image and recovered images using different estimated PSFs and computation time are listed in Tables 1 and 2, respectively. As can be seen from the result data, our scheme cannot only obtain a good recovery result but also shorten the compute time remarkably.

Table 1

SNR using different estimated PSFs in each size in Figure 5.

PSF size	SNR of degrade images (dB)	SNR using different estimated PSFs (dB)
		Fergus et al. [13]	Levin et al. [15]	Hu et al. [16]	Ours
17 × 17	8.32	12.83	13.06	12.20	12.78
21 × 21	8.90	9.32	9.47	9.17	9.35
23 × 23	5.34	6.20	5.79	5.86	6.23
27 × 27	4.16	5.06	4.90	4.87	5.11

Table 2

Processing time using different algorithms for PSFs in each size in Figure 5.

PSF size	Processing time of references (s)			Processing time of our scheme (s)
	Fergus et al. [13]	Levin et al. [15]	Hu et al. [16]	Kernel estimation	Image deconvolution	Total
13 × 13	30.52	127.80	6.32	1.45	2.57	4.02
17 × 17	33.93	138.57	7.69	1.57	3.76	5.33
21 × 21	44.03	146.95	8.73	1.69	4.41	6.10
23 × 23	53.47	154.20	8.91	1.78	4.45	6.23
27 × 27	56.20	155.79	9.48	1.89	5.38	7.27

Figure 8 shows a real blur image deblurring result using above algorithms. It is a large scale image obtained from [13] with size 1000 × 1256, as shown in Figure 8(a). Figure 8(b) shows Fergus's result, whose PSF is estimated by varying image resolution in a coarse-to-fine manner, and then the estimated PSF is used to recover the image by the standard Richardson-Lucy algorithm. Figure 8(c) shows Levin's result by alternating between estimating the PSF and solving for the image. Both of their methods obtain good results but are still time consuming, especially for the large scale image, which will cost several minutes. Our result is shown in Figure 8(d); extracted and zoomed color squares are shown in Figure 8(e). It can be seen that our method gets comparable result and only costs ten seconds.

Figure 8

Experiment on the real blur image. (a) The blurred image; (b) Fergus et al. [13] result, $t=129$  s; (c) Levin et al. [15] result, $t=407$  s; (d) our result, $t=12$  s; (e) color squares from results are zoomed.

5. Conclusion

Traditional single image deblurring algorithms require presetting or estimating the size of PSF according to the iteration stopping criterion, which is time consuming. In this paper, an efficient PSF estimated strategy has been proposed based on image gradient cepstrum analysis. PSF power information is obtained from cepstrum properties of the blur image gradient and its phase information is retrieved by PR technique. To verify the accuracy of the strategy, the estimated PSF is used for recovering image with a novel total variation image restoration model that is coupling with a gradient fidelity term. An alternating direction method (ADM) numerical algorithm with rapid convergence and high robustness is adopted to optimize the energy function. Both synthetic and real blur experiments show that our scheme can rapidly produce comparable results with some state-of-the-art algorithms and greatly shorten the computing time.