Split Bregman Iteration Algorithm for Image Deblurring Using Fourth-Order Total Bounded Variation Regularization Model

Abstract

We propose a fourth-order total bounded variation regularization model which could reduce undesirable effects effectively. Based on this model, we introduce an improved split Bregman iteration algorithm to obtain the optimum solution. The convergence property of our algorithm is provided. Numerical experiments show the more excellent visual quality of the proposed model compared with the second-order total bounded variation model which is proposed by Liu and Huang (2010).

1. Introduction

Image restoration problem is one of the earliest and most classic linear inverse problems [1–3]. In this class of problems, a noisy indirect observation 𝑓 of an original image 𝑢 is modeled as 𝑓=𝐴𝑢+𝑛, (1) where 𝐴 is a bounded linear operator representing the convolution and 𝑛 denotes the additive noise.

Equation (1) is a typically ill-posed inverse problem; that is, a small change in 𝑓 will lead to huge deviation in the solution 𝑢. Hence, to keep the numerical stability, the regularization method known as adding a regularization term to the energy function has been developed. The original scheme introduced by Tikhonov and Arsenin [4] is given by min_𝑢Ω|𝐿𝑢|²+𝜆2𝐴𝑢−𝑓²₂, (2) where the nonnegative function ∫_Ω|𝐿𝑢|² which regularizes the solution by enforcing certain prior constrains on original image is known as regularization/penalty and the ‖𝐴𝑢−𝑓‖²₂ which measures the violation of the relation between 𝑢 and its observation 𝑓 is known as fidelity. The scale 𝜆 is called the regularization parameter; it compromises fidelity with penalty. By this Tikhonov regularization method, we can compute stable approximations to the original solution. However, this smoothness penalty model does not preserve the edge, sparsity pattern, and texture well because of its isotropic smoothing properties. To overcome this shortcoming, Rudin et al. [5] proposed the total variation (TV) based regularization scheme (the ROF model) min_𝑢Ω||+𝜆𝐿𝑢2𝐴𝑢−𝑓²₂, (3) where Ω⊆ℝ² denotes a bounded subset with Lipschitz boundary, 𝑢∈𝕃¹(Ω), and 𝐿𝑢 represents the distributional derivative of 𝑢. Many computational methods [5–11] for solving (3) sprang up in recent years. For the moment we may think of the nonsmooth penalty term ∫_Ω|𝐿𝑢| as the 𝕎^1,1(Ω) seminorm. More precisely, _Ω|𝐿𝑢|=sup_Ω𝑢div𝜑d𝑥∣𝜑∈𝐶¹_𝐶Ω,ℝ^𝑑,𝜑≤1. (4) As a result, the bounded variation (BV) seminorm is endowed with ‖𝑢‖_BV=‖𝑢‖_𝕃1+∫_Ω|𝐿𝑢|. Then the Banach space 𝔹𝕍(Ω) is essentially an extension of 𝕎^1,1(Ω). This model has been extremely successful in a wide variety of image restoration problems, such as image denoising [5, 12], signal processing [13, 14], image deblurring [15], image decomposition, and texture extraction [16]. However, this model also usually produces staircase effects and new edges that do not exist in the true image.

In [17, 18], the authors concentrate especially on the full norm ‖𝑢‖_BV+𝜌‖𝑢‖²₂ as the regularization term, compared with the TV regularization, which is preferable due to its ability to preserve edges in the original image during the reconstruction process. Then specify the original problem to the following variation model [18]: min_𝑢Ω||+𝛼𝐿𝑢2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂,over𝑢∈𝕂(Ω)∩𝕏(Ω), (5) where 𝛼≥0,𝛽>0,𝕂 is a closed, convex subset of 𝕃²(Ω),and𝕏=𝕃²(Ω)∩𝔹𝕍(Ω). The space 𝕏 endowed with the norm ‖𝑢‖_𝕏=‖𝑢‖_𝕃²+‖𝑢‖_𝔹𝕍 is a Banach space. A quadratic regularization term is utilized in model (5) comparing with model (3), it serves two advantages. One advantage is that, for 𝛼>0, it provides a coercive term for the subspace of constant function which is in the kernel of the 𝐿-operator (in this model 𝐿 represents the gradient operator). The other advantage of the quadratic regularization term is that it provides the probability to discriminate the structure of stability results from that of the nonquadratic 𝔹𝕍 term.

As mentioned in [19–21], this technique preserves edges well, but the obtained images for this model are often piecewise constant. To prevent the staircase effect, we penalize jumps more; this can be achieved by taking the second derivation into account. So in this paper we present an improved model; that is, the second-order diffusive term is replaced by fourth-order diffusive term in the model (5); new model substantially reduces the staircase effect, while preserving sharp jump discontinuities (edges). Here we rewrite the proposed model: min_𝑢Ω∇²𝑢+𝛼2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂,over𝑢∈𝕃²(Ω)∩𝔹𝕍²(Ω), (6) where the Frobenius norm of the Hessian matrix ∇²𝑢 is (see [22]) ∇²𝑢=𝑢²_𝑥𝑥+𝑢²_𝑦𝑦+𝑢²_𝑥𝑦+𝑢²_𝑦𝑥^1/2, (7) and the 𝔹𝕍² is defined as follows [23].

Definition 1. Let Ω⊂ℝ² be an open subset with Lipschitz boundary. 𝔹𝕍²(Ω) is a subspace of functions 𝑢∈𝕃¹(Ω) such that the following equation is satisfied: _Ω∇²𝑢=sup_Ω𝑢div²𝜓d𝑥∣𝜓∈𝐶²_𝑐Ω,ℝ^2×2,𝜓≤1<∞, (8) where div²𝜓=²_ℎ,𝑘=1𝜕_𝑘𝜕_ℎ𝜓^ℎ𝑘,𝜓(𝑥)=²_ℎ,𝑘=1𝜓^ℎ𝑘2. (9)𝐶²_𝑐(Ω) stands for the set of functions in 𝐶²(Ω) with compact support in Ω.

The proof of the existence, uniqueness, convergence, and stability of our proposed model (6) can be founded in [18].

The organization of the rest of paper is as follows. In Section 2, we give a detailed description of the Bregman iteration method. In Section 3, we elaborate on the analysis of the extended split Bregman iteration method for the proposed model. In Section 4, the convergence analysis is displayed. Numerical experiments intended for demonstrating the effectiveness of our model are provided in Section 5. Finally, concluding remarks are given in Section 6.

2. Bregman-Related Algorithms

Bregman iteration is a concept that originated in function analysis for finding extrema of convex function [24], which was initially introduced and studied by Osher et al. for image processing [25]. Now we will show the general formulation of Bregman iteration technique.

2.1. Bregman Iteration

In [26], Goldstein and Osher considered the generalized constrained minimizations of the form min_𝑢𝐽(𝑢)subjectto𝐻(𝑢)=0, (10) where 𝐽 and 𝐻 defined in ℝ^𝑛 are convex functions. The associated unconstrained problem is min_𝑢𝐽(𝑢)+𝜆𝐻(𝑢), (11) where 𝜆 is the positive parameter which should be chosen extremely large.

Definition 2. The Bregman distance of function 𝐽 between 𝑢 and 𝑣 is 𝐷^𝑝_𝐽(𝑢,𝑣)=𝐽(𝑢)−𝐽(𝑣)−⟨𝑝,𝑢−𝑣⟩,𝑝∈𝜕𝐽(𝑣), (12) where ⟨⋅,⋅⟩ stands for duality product and 𝑝 is in the subdifferential of 𝐽 at 𝑣 with 𝜕𝐽(𝑣)∶=𝑝∈BV(Ω)^∗∣𝐽(𝑢)≥𝐽(𝑣)+𝑝,𝑢−𝑣. (13) We assume that 𝐻 is differentiable then problem (11) can be iteratively solved by [24] 𝑢^𝑘+1=argmin_𝑢𝐷^𝑝𝑘_𝐽𝑢,𝑢^𝑘(𝑢),𝑝+𝜆𝐻_𝑘+1=𝑝^𝑘−𝜆∇𝐻𝑢^𝑘+1. (14) The convergence analysis of this Bregman iterative scheme was provide in [25]. The computational performance of Bregman iteration relays on how fast the subproblem argmin_𝑢𝐷^𝑝𝑘_𝐽𝑢,𝑢^𝑘+𝜆𝐻(𝑢) (15) can be solved. Let 𝐻(𝑢,𝑓)=(1/2)‖𝐴𝑢−𝑓‖²₂, where 𝐴 is a linear operator. As show in [25, 27], iteration (14) can be reformulated as a simplified form 𝑢^𝑘+1=argmin_𝑢𝐽(𝑢)+𝜆𝐻𝑢,𝑓^𝑘,𝑓^𝑘+1=𝑓^𝑘+𝑓−𝐴𝑢^𝑘+1,𝑓₀=𝑓. (16) This Bregman iteration which was proposed by Osher et al. for TV based image denoising [25] has two advantages [26]; the first one is that this method converges very quickly, especially for problem where 𝐽(𝑢) contains 𝕃1-regularization term. The second advantage is that the parameter 𝜆 in (14) remains constant; so for the purpose of fast convergence, we can choose 𝜆 which minimizes the condition number of the subproblem. Due to the highefficiency and robustness of Bregman iteration method, it has been widely used for image reconstruction [27–29].

2.2. Split Bregman

Goldstein and Osher [26] proposed the split Bregman iteration based on the split formulation provided in [30] to attack the general 𝕃1-regularized optimization problem (11). In a recent paper [31], this method is used to solve general variational models for image restoration. The crucial point of the split Bregman method is that it could separate the 𝕃1 and 𝕃2 portions in (11). They converted (11) into the constrained optimization problem min_𝑢‖𝑑‖₁+𝜆𝐻(𝑢)suchthat𝑑=𝜙(𝑢), (17) where 𝐻(𝑢) and 𝜙(𝑢) stand for the convex functions. Then transform it into an unconstrained problem min_𝑢‖𝑑‖₁(𝑢)+𝛾+𝜆𝐻2𝑑−𝜙(𝑢)²₂. (18) This problem is similar to (11); so they enforced the simplified Bregman iterative algorithm to problem (18) 𝑢^𝑘+1,𝑑^𝑘+1=min_𝑢,𝑑‖𝑑‖₁(𝑢)+𝛾+𝜆𝐻2𝑑−𝜙(𝑢)−𝑏^𝑘2₂,𝑏^𝑘+1=𝑏^𝑘+𝜙𝑢^𝑘+1−𝑑^𝑘+1. (19) This is called two-phase split Bregman iterative algorithm.

Then pay attention to the subproblem 𝑢^𝑘+1,𝑑^𝑘+1=min_𝑢,𝑑‖𝑑‖₁(𝑢)+𝛾+𝜆𝐻2𝑑−𝜙(𝑢)−𝑏^𝑘2₂. (20) Due to the “split” of the 𝕃1 and 𝕃2 components of this function, the minimization problem was performed by iteratively minimizing with respect to 𝑢 and 𝑑 separately Step  1:    𝑢^𝑘+1=min_𝑢𝛾𝜆𝐻(𝑢)+2‖𝑑^𝑘−𝜙(𝑢)−𝑏^𝑘‖²₂, Step  2:    𝑑^𝑘+1=min_𝑢‖𝑑‖₁+𝛾2‖𝑑−𝜙(𝑢^𝑘+1)−𝑏^𝑘‖²₂.The speed of the split Bregman method is relayed on how quickly the two steps can be solved. A wide variety optimization techniques can be used to solve Step 1, for instance, the Fourier transform method and conjugate gradient method. Step 2 could be attacked with the extremely fast shrinkage formula, namely, 𝑑_𝑗^𝑘+1=shrink𝜙(𝑢)_𝑗+𝑏^𝑘_𝑗,1𝛾, (21) where ()=𝑥shrink𝑥,𝜏|𝑥|∗max(|𝑥|−𝜏,0). (22)

3. The Proposed Algorithm

Due to [5], 𝐿₁ estimation procedures are more appropriate for image restoration, and TV norm is essentially 𝐿₁ norm of derivatives. However, the fourth-order total variation term ∫_Ω|∇²𝑢| in model (6) is continuous expression; by utilizing the similar method of discretizing total variation in [32], we can deduce the discrete fourth-order total variation (let 𝐷 represent the operator ∇²): _Ω|𝐷𝑢|=‖𝐷𝑢‖₁. (23)

So the proposed model (6) can be rewritten as follows: min_𝑢‖𝐷𝑢‖₁+𝛼2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂. (24)

Then we perform the split Bregman iteration to solve problem (24); this yields a constrained problem min_𝑢‖𝑑‖₁+𝛼2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂,s.t.𝑑=𝐷𝑢. (25) Obviously, problem (25) is equivalent to the following unconstraint problem: min_𝑢‖𝑑‖₁+𝛼2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂+𝜆2‖𝑑−𝐷𝑢‖²₂. (26) Concretely, the extended split Bregman iterative for solving (26) is depicted as 𝑢^𝑘+1=argmin_𝑢𝛼2‖𝑢‖²₂+𝛽2𝐴𝑢−𝑓²₂+𝜆2𝐷𝑢−𝑑^𝑘+𝑏^𝑘2₂,𝑑^𝑘+1=argmin_𝑢‖𝑑‖₁+𝜆2𝑑−𝐷𝑢^𝑘+1−𝑏^𝑘2₂, (27) with the update formula for 𝑏^𝑘𝑏^𝑘+1=𝑏^𝑘+𝐷𝑢^𝑘+1−𝑑^𝑘+1. (28)

For more precisely, we derive our algorithm as follows.

Algorithm 3. We have the following steps.(1) Set input value 𝑢⁰=0,𝑑⁰=𝑏⁰=0. Set 𝑛∶=0.(2) Update 𝑢^𝑘+1 from 𝑢^𝑘+1=𝛼𝐼+𝛽𝐴^𝑇𝐴+𝜆𝐷^𝑇𝐷⁻¹×𝛽𝐴^𝑇𝑓+𝜆𝐴^𝑇𝑑^𝑘−𝑏^𝑘. (29)(3) Update 𝑑^𝑘+1 from 𝑑^𝑘+1=shrink𝐷𝑢^𝑘+1+𝑏^𝑘,1𝜆. (30)(4) Update 𝑏^𝑘+1 from 𝑏^𝑘+1=𝑏^𝑘+𝐷𝑢^𝑘+1−𝑑^𝑘+1. (31)(5) If stopping criterion holds, output the 𝑢^𝑘+1. Otherwise, set 𝑛∶=𝑛+1 and go to Step 2.

4. Convergence Analysis

In this section, we concentrate on the rigorous convergence proof of our iterative algorithm. Our analysis below is similar to that in [33, 34], where the authors presented the analysis of the unconstrained split Bregman iteration in detail.

We note that the subproblems (27) are convex and differentiable; so the first-order optimality conditions for 𝑢^𝑘+1 and 𝑑^𝑘+1 are obtained as follows: 0=𝛼𝑢^𝑘+1+𝛽𝐴^𝑇𝐴𝑢^𝑘+1−𝑓+𝜆𝐷^𝑇𝐷𝑢^𝑘+1−𝑑^𝑘+𝑏^𝑘,0=𝑝^𝑘+1+𝜆𝑑^𝑘+1−𝐷𝑢^𝑘+1−𝑏^𝑘,𝑏^𝑘+1=𝑏^𝑘+𝐷𝑢^𝑘+1−𝑑^𝑘+1, (32) where 𝑝^𝑘∈𝜕‖𝑑^𝑘‖₁. The condition (32) will be used for analyzing the convergence property of our algorithm.

Theorem 4. Assume that there exists a unique solution 𝑢^∗ of problem (26) and 𝛼>0,𝛽>0,and𝜆>0. Then the sequence {𝑢^𝑘} generated by extended split Bregman iteration (27)–(40) holds: lim_𝑘→∞𝛼2𝑢^𝑘2₂+𝛽2𝐴𝑢^𝑘−𝑓²₂+𝐷𝑢^𝑘₁=𝛼2𝑢^∗2₂+𝛽2𝐴𝑢^∗−𝑓²₂+𝐷𝑢^∗₁,lim_𝑘→∞𝑢^𝑘−𝑢^∗₂=0. (33)

Proof. Reference [17] has shown that problem (26) exists a unique 𝑢^∗ so the first order optimality condition holds 0=𝐷^𝑇𝑝^∗+𝛼𝑢^∗+𝛽𝐴^𝑇𝐴𝑢^∗−𝑓, (34) where 𝑝^∗∈𝜕‖𝑑^∗‖₁ and 𝑑^∗=𝐷𝑢^∗; let 𝑏^∗=𝑝^∗/𝜆; then (32) is equivalent to 0=𝛼𝑢^∗+𝛽𝐴^𝑇𝐴𝑢^∗−𝑓+𝜆𝐷^𝑇𝐷𝑢^∗−𝑑^∗+𝑏^∗,0=𝑝^∗+𝜆𝑑^∗−𝐷𝑢^∗−𝑏^∗,𝑏^∗=𝑏^∗+𝐷𝑢^∗−𝑑^∗. (35) This illustrates that 𝑢^∗,𝑑^∗, and 𝑏^∗ are a fixed points of (27).
Let 𝑢^𝑘,𝑑^𝑘, and𝑏^𝑘 denote the sequence generated by algorithm (27), and 𝑢^𝑘_𝑒=𝑢^𝑘−𝑢^∗,𝑑^𝑘_𝑒=𝑑^𝑘−𝑑^∗, and 𝑏^𝑘_𝑒=𝑏^𝑘−𝑏^∗ represent the errors, respectively. Then subtracting every equation in (32) from the corresponding equations in (35), we give the result 0=𝛼𝑢_𝑒^𝑘+1+𝛽𝐴^𝑇𝐴𝑢_𝑒^𝑘+1+𝜆𝐷^𝑇𝐷𝑢_𝑒^𝑘+1−𝑑^𝑘_𝑒+𝑏^𝑘_𝑒,0=𝑝^𝑘+1−𝑝^∗+𝜆𝑑_𝑒^𝑘+1−𝐷𝑢_𝑒^𝑘+1−𝑏^𝑘_𝑒,𝑏_𝑒^𝑘+1=𝑏^𝑘_𝑒+𝐷𝑢_𝑒^𝑘+1−𝑑_𝑒^𝑘+1, (36) where 𝑝^𝑘+1∈𝜕‖𝑑^𝑘+1‖₁; we take the inner product of the first and second equations in (36) with respect to 𝑢^𝑘_𝑒 and 𝑑^𝑘_𝑒 separately; then followed by the same manipulations applied to the third equation of (36), we obtain 0=𝛼𝑢_𝑒^𝑘+12₂+𝛽𝐴𝑢_𝑒^𝑘+12₂+𝜆𝐷𝑢_𝑒^𝑘+12₂+𝜆𝐷^𝑇𝑏^𝑘_𝑒−𝐷^𝑇𝑑^𝑘_𝑒,𝑢_𝑒^𝑘+1,0=𝑝^𝑘+1−𝑝^∗,𝑑_𝑒^𝑘+1+𝜆𝑑_𝑒^𝑘+12₂−𝜆𝐷𝑢_𝑒^𝑘+1+𝑏^𝑘_𝑒,𝑑_𝑒^𝑘+1,𝑏_𝑒^𝑘+12₂−𝑏^𝑘_𝑒²₂=2𝑏^𝑘_𝑒,𝐷𝑢_𝑒^𝑘+1−𝑑_𝑒^𝑘+1+𝐷𝑢_𝑒^𝑘+1−𝑑_𝑒^𝑘+12₂. (37) Obviously, the third equation of (37) can be rewritten as 𝜆2𝑏^𝑘_𝑒²₂−𝑏_𝑒^𝑘+12₂𝜆=−2𝐷𝑢_𝑒^𝑘+1−𝑑_𝑒^𝑘+12₂−𝜆𝑏^𝑘_𝑒,𝐷𝑢_𝑒^𝑘+1−𝑑_𝑒^𝑘+1. (38) Combining the first two equations of (37) with (38), we have 𝜆2𝑏^𝑘_𝑒²₂−𝑏_𝑒^𝑘+12₂+𝜆2𝑑^𝑘_𝑒²₂−𝑑_𝑒^𝑘+12₂=𝜆2𝐷𝑢_𝑒^𝑘+1−𝑑^𝑘_𝑒²₂+𝑝^𝑘+1−𝑝^∗,𝑑^𝑘+1−𝑑^∗+𝛼𝑢_𝑒^𝑘+12₂+𝛽𝐴𝑢_𝑒^𝑘+12₂. (39) By summing the equation bilaterally from 𝑘=0 to 𝑁, we get 𝜆2𝑏⁰_𝑒²₂−𝑏_𝑒^𝑁+12₂+𝜆2𝑑⁰_𝑒²₂−𝑑_𝑒^𝑁+12₂=𝜆2^𝑁_𝑘=0𝐷𝑢_𝑒^𝑘+1−𝑑^𝑘_𝑒²₂+^𝑁_𝑘=0𝑝^𝑘+1−𝑝^∗,𝑑^𝑘+1−𝑑^∗+𝛼^𝑁_𝑘=0𝑢_𝑒^𝑘+12₂+𝛽^𝑁_𝑘=0𝐴𝑢_𝑒^𝑘+12₂. (40) Because 𝑝^𝑘+1∈𝜕‖𝑑^𝑘+1‖₁,𝑝^∗∈𝜕‖𝑑^∗‖₁, and ‖⋅‖₁ is convex, then we have 𝑝^𝑘+1−𝑝^∗,𝑑^𝑘+1−𝑑^∗≥0∀𝑘. (41) The fact that all terms involved in (40) are nonnegative leads to the following expression: 𝜆2𝑏⁰_𝑒²₂+𝑑⁰_𝑒²₂≥𝜆2^𝑁_𝑘=0𝐷𝑢_𝑒^𝑘+1−𝑑^𝑘_𝑒²₂+𝛼^𝑁_𝑘=0𝑢_𝑒^𝑘+12₂+𝛽^𝑁_𝑘=0𝐴𝑢_𝑒^𝑘+12₂+^𝑁_𝑘=0𝑝^𝑘+1−𝑝^∗,𝑑^𝑘+1−𝑑^∗. (42) With the assumption 𝛼>0,𝛽>0,and𝜆>0 and letting 𝑁→∞, inequality (42) suggests the following conclusions: ^+∞_𝑘=0𝑢_𝑒^𝑘+12₂<+∞,_+∞𝑘=0𝐴𝑢_𝑒^𝑘+12₂<+∞,_+∞𝑘=0𝐷𝑢_𝑒^𝑘+1−𝑑^𝑘_𝑒²₂<+∞,_+∞𝑘=0𝑝^𝑘+1−𝑝^∗,𝑑^𝑘+1−𝑑^∗<+∞. (43) And hence lim_𝑘→+∞𝑢^𝑘+1−𝑢^∗₂=0,lim_𝑘→+∞(𝐴𝑢^𝑘+1−𝑓)−(𝐴𝑢^∗−𝑓)₂=0,lim_𝑘→+∞𝐷𝑢^𝑘+1−𝑑^𝑘₂=0,lim_𝑘→+∞𝑝^𝑘−𝑝^∗,𝑑^𝑘−𝑑^∗=0. (44)(45)(46)(47) The Bregman distance satisfies 𝐷^𝑝_𝐽(𝑢,𝑣)+𝐷^𝑞_𝐽()=,(𝑣)(𝑢).𝑣,𝑢𝑞−𝑝,𝑢−𝑣∀𝑝∈𝜕𝐽,𝑞∈𝜕𝐽 (48) Associating this equation with (47), we get lim_𝑘→+∞𝑑^𝑘₁−𝑑^∗₁−𝑝^∗,𝑑^𝑘−𝑑^∗=0. (49) Again, together with (46) and the continuous property of ‖⋅‖₁, the following expression hold: lim_𝑘→+∞𝐷𝑢^𝑘₁−𝐷𝑢^∗₁−𝑝^∗,𝐷𝑢^𝑘−𝐷𝑢^∗=0,𝑝^∗∈𝜕𝐷𝑢^∗₁. (50) With the similar skills, we have lim_𝑘→+∞𝐴𝑢^𝑘−𝑓²₂−𝐴𝑢^∗−𝑓²₂−2𝐴^𝑇𝐴𝑢^∗−𝑓,𝑢^𝑘−𝑢^∗=0,lim_𝑘→+∞𝑢^𝑘2₂−𝑢^∗2₂=lim_𝑘→+∞2𝑢^∗,𝑢^𝑘−𝑢^∗. (51) Combining (50)-(51), we obtain lim_𝑘→+∞𝐷𝑢^𝑘₁+𝛼2𝑢^𝑘2₂+𝛽2𝐴𝑢^𝑘−𝑓²₂−𝐷𝑢^∗₁+𝛼2𝑢^∗2₂+𝛽2𝐴𝑢^∗−𝑓²₂−⟨𝐷^𝑇𝑝^∗+𝛼𝑢^∗+𝛽𝐴^𝑇𝐴𝑢^∗−𝑓,𝑢^𝑘−𝑢^∗⟩=0. (52) Finally, from (44) and (52), the main results are shown as follows: lim_𝑘→+∞𝐷𝑢^𝑘₁+𝛼2𝑢^𝑘2₂+𝛽2𝐴𝑢^𝑘−𝑓²₂=𝐷𝑢^∗₁+𝛼2𝑢^∗2₂+𝛽2𝐴𝑢^∗−𝑓²₂,lim_𝑘→+∞𝑢^𝑘−𝑢^∗₂=0. (53) This proves Theorem 4.

5. Numerical Experiments

In this section, a number of experiments are performed to demonstrate the effectiveness and efficiency of our proposed split Bregman iteration (PSBI) algorithm for the fourth-order diffusive model (6), which will be compared with the split Bregman iteration (SBI) for second-order diffusive model provided in [17]. All experiments are generated in MATLAB 7.10 environment on a desktop with Windows 7 operating system, 3.00 GHz Intel Pentium(R) D CPU, and 1.00 GB memory.

The performance of all algorithms is measured by the improved signal-to-noise ratio (ISNR) and mean squared error (MSE) defined as ISNR=10log₁₀𝑓−𝑢₀²𝑢−𝑢₀²,MSE=𝑢−𝑢₀²,𝑚𝑛 (54) where 𝑓,𝑢₀, and 𝑢 denote the degraded, original, and recovered images, and 𝑚 and 𝑛 are the sizes of image, respectively. The higher ISNR or the lower MSE, the better quality of the deblurred images. Moreover, the stopping criterion of all algorithms measured by the difference between the consecutive iterations of the deblurred image satisfies the following inequality: 𝑢^𝑘+1−𝑢^𝑘₂𝑢^𝑘+1₂≤5×10⁻⁴. (55)

Three classical grayscale images with resolution of 256×256 pixels in Figure 1 are used for synthetic degradations in our experiments. The blur kernels used for blurring are Gaussian blur (size of 7×7 pixels, variance of 3), linear motion blur (length of 10 pixels and direction of 30 degrees), out-of-focus blur (size of 10 × 10 pixels, defocus radius of 4) and uniform blur (size of 7×7 pixels). Periodic boundary conditions are considered in the following experiments.

(a)

(b)

(c)

Figure 1 

Images used for synthetic degradations: (a) “Cameraman,” (b) “Boat,” (c) “Lena.”

In the first experiment, the original image “Cameraman” with complex background is blurred by Gaussian blur, out-of-focus blur, linear motion blur, and uniform blur, respectively. The blurred images are showed in Figures 2(a)–2(d), Figures 2(e)–2(h) show the deblurred results corresponding to the split Bregman iteration method for all the blur cases, and the selected parameter values corresponding to the SBI method are the same with those the authors given out in [17]. In order to give explicit comparison, we show the restored images obtained by our PSBI method in Figures 2(i)–2(l), and the selected parameter values are 𝛼=1𝑒+5,𝛽=1𝑒+12,and𝜆=2.1𝑒+6, respectively. Through comparing the “sky” in those deblurred images we can see that the PSBI method recovers more details (removes the stair effect) than SBI method. However, due to the fourth-order diffusive term, our proposed method costs a little more time than the SBI method. Table 1 gives the CPU times, ISNR values, and MSE values of the SBI method and PSBI method for all the blur kernels mentioned above, which shows that our algorithm works well for images with complex background.

Table 1 

Computational cost, ISNR value, and MSE value for different deblur cases.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 2 

Deblurring of “Cameraman." (a)–(d) Gaussian blurred image, linear motion blurred image, out-of-focus blurred image and uniform blurred image, respectively. (e)–(h) Images deblurred by SBI. (i)–(l) Images deblurred by PSBI.

In the next test, we use Gaussian blur and out-of-focus blur to degrade the “Boat” image and then run the two algorithms many times to obtain the best results. In PSBI method, the selected parameters and iterations are 𝛼=1.1𝑒+5,𝛽=1𝑒+12,𝜆=2𝑒+6,andite=46 for Gaussian blur and 𝛼=1𝑒+5,𝛽=1𝑒+12,𝜆=2.15𝑒+6,andite=25 for out-of-focus blur, respectively. Figure 3 shows the recovered results of the two methods. And for the purpose of better illustrations, we provide a close-up of image region in Figure 4. We can see from Figure 4 that the restored images estimated by PSBI are better when compared to the ones deblurred by SBI; for example, see the “letter” in the stern of “Boat.” Meanwhile, Table 2 illustrates the ISNR values obtained by PSBI are higher than those obtained by SBI while the MSE values are lower. So, it is clear that our method can effectively reduce the “staircase” which always appeals in the two-order diffusive model.

Table 2 

Image deblurring using SBI and PSBI for “Boat.”

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Comparison with SBI method. First column: out-of-focus blurred image and Gaussian blurred image. Second column: images deblurred by SBI. Third column: images deblurred by PSBI.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Close-ups of selected section of Figure 3. First column: out-of-focus blurred image and Gaussian blurred image. Second column: images deblurred by SBI. Third column: images deblurred by PSBI.

Figures 5 and 6 show the comparison between PSBI and SBI algorithm for image “Lena.” Firstly we blur the “Lena” by linear motion blur and uniform blur and then select the deblurred result that looks best after carefully tuning the parameters; for linear motion blur, the parameter values and iterations are 𝛼=1𝑒+05,𝛽=1.15𝑒+12,𝜆=2.2𝑒+6,andite=28, and for uniform blur, they are 𝛼=1𝑒+05,𝛽=1𝑒+12,𝜆=2.1𝑒+6,andite=25. It can be seen from Figures 5 and 6 that many staircase effects in the smooth regions can be detected obviously in the images restored by the SBI model while seldom appear in our results, and the upper-left corners in image “Lena” deblurred by PSBI are more visual comparing with that obtained by SBI model. Table 3 shows that the restored images obtained by our model have higher ISNR values and lower MSE values than the second-order model.

Table 3 

Image deblurring using SBI and PSBI for “Lena.”

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Comparison with SBI method. First column: motion blurred image and uniform blurred image. Second column: images deblurred by SBI. Third column: images deblurred by PSBI.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Close-ups of selected section of Figure 5. First column: motion blurred image and uniform blurred image. Second column: images deblurred by SBI. Third column: images deblurred by PSBI.

6. Conclusion

In this paper, we propose the fourth-order total bounded variation regularization based image deblurring model and exploit the split Bregman iteration method to solve this new model. The convergence analysis of our algorithm is provided. Numerical experiments show that our algorithm works well for images with complex background or simple background. In our synthetic experiments, the fourth-order diffusive model yields better results than the second-order diffusive model. It is believed that the proposed model can be extended to further applications in image processing and computer vision.

Acknowledgments

The authors would like to thank the referees very much for their valuable comments and suggestions.