Efficient Algorithm for Isotropic and Anisotropic Total Variation Deblurring and Denoising

Abstract

A new deblurring and denoising algorithm is proposed, for isotropic total variation-based image restoration. The algorithm consists of an efficient solver for the nonlinear system and an acceleration strategy for the outer iteration. For the nonlinear system, the split Bregman method is used to convert it into linear system, and an algebraic multigrid method is applied to solve the linearized system. For the outer iteration, we have conducted formal convergence analysis to determine an auxiliary linear term that significantly stabilizes and accelerates the outer iteration. Numerical experiments demonstrate that our algorithm for deblurring and denoising problems is efficient.

1. Introduction

The purpose of image restoration is to recover original image 𝑢 from an observed data 𝑧 (noisy and blurred image) from the relation 𝑧=𝐾𝑢+𝑛, (1) where 𝐾 is a known linear blurring operator and 𝑢 is a Gaussian white noise. The well-known minimization problem for image denoising and deblurring based on isotropic total variation (TV) is proposed by Rudin et al. [1]: min_𝑢𝛼_Ω|∇𝑢|²1+𝛽𝑑𝑥𝑑𝑦+2‖𝐾𝑢−𝑧‖²_𝐿2. (2) Here 𝛼>0 is the penalty parameter, and 𝛽>0 is the diffusion regularization parameter and is typically small. The functional in (2) is strictly convex with a unique global minimizer. The paper [2] showed the well-posedness of problem (2) as 𝛽→0⁺ and the existence of a unique global minimizer. The corresponding Euler-Lagrange equation for (2) is given by −𝛼∇⋅∇𝑢|∇𝑢|²+𝛽+𝐾^∗(𝐾𝑢−𝑧)=0, (3) where 𝐾^∗ is the adjoint operator of 𝐾 with respect to standard 𝐿₂ inner product.

In addition to the isotropic models, the anisotropic models for a qualitative improvement at corners are shown in [3] with that a transfer of the discretization from the anisotropic model to the isotropic setting results in an improvement of rotational invariance. The anisotropic TV defined in [3] is TV₂(𝑢)=_Ω∇_𝑥𝑢+∇_𝑦𝑢𝑑𝑥𝑑𝑦. (4) The existence and uniqueness of solutions of the anisotropic total variation flow are shown in [4]. The explicit time-marching schemes are applied in [1, 5]. Zuo and Lin [6] proposed a generalized accelerated proximal gradient (GAPG) approach for solving TV-based image restoration problems by replacing the Lipschitz constant with an appropriate positive-definite matrix, resulting in faster convergence, and the TV regularization can be either isotropic or anisotropic. The convergence rate of 𝑂(𝑘²) is maintained by GAPG, where 𝑘 is the number of iterations. The introduction of an anisotropy to TV regularization [7] indeed leads to improved denoising: the staircasing effect is reduced while at the same time the creation of artifacts is suppressed.

There are different algorithms to solve (3). Earlier works include the time-marching scheme [1], the affine scaling algorithm [8], the Newton’s method with a continuation procedure on 𝛽 [9], a fixed point iteration with multigrid method [10–12], the combination of fixed point iteration and preconditioned conjugate gradient which method is proposed in [13], a fast algorithm based on the discrete cosine transform [14, 15], a new modular solver [16], a proximity algorithm [17], and a nonlinear primal-dual method [18]. Though the corresponding linear system of the nonlinear primal-dual method is twice as large as the primal system, it is shown to be better conditioned. In [19], the combination of the algebraic multigrid (AMG) methods [20, 21], Krylov subspace acceleration, and extrapolation of initial data are successfully combined to give an efficient algorithm for the purely denoising problem. In [22], the authors proposed the split Bregman method to solve image denoising and compressed sensing problem rapidly. Chang et al. [23] added a linear matrix to speed up the computational process.

In this paper, we propose an efficient and rapid algorithm for minimization problem (2) rather than solving the Euler-Lagrange equation (3) directly which combines the split Bregman method, the algebraic multigrid method, and Krylov subspace acceleration. Our algorithm consists of two steps. One uses the split Bregman method to convert (2) into linear system in the outer iteration, and the other adopts an AMG method previously developed by one of the authors for other applications [23, 24] to solve a linear system in the inner iteration. Moreover, the outer iteration is accelerated by Krylov subspace extrapolation. One 𝑉-cycle per inner iteration is sufficient in all our simulations. For the outer iteration, we have developed a stabilizing technique adapted to the blur operator 𝐾. This is done by adding a linear term on both sides of the equation [23]. Motivated by that, a wider linear term different from that in [23] is considered here. This linear stabilizing term is obtained via formal convergence analysis and is expressed explicitly in terms of the mask of 𝐾. The inclusion of the linear stabilizing term plays a crucial role in our scheme. The outer iteration indeed may diverge without a proper linear stabilizing term (Table 6, Section 5).

The rest of the paper is organized as follows. In Section 2, we discuss the purely image denoising problem by describing briefly the AMG method and the split Bregman method for the isotropic TV model and the anisotropic TV model. In Section 3, we discuss the image denoising and deblurring problem and present the formal convergence analysis and the framework of our scheme. In Section 4, we give out explicit formulae of the linear stabilizing term, together with other implementation details. Numerical results are also compared among several algorithms in Section 5. We end this paper with a brief summary in Section 6.

2. The Algebraic Multigrid

In this section, we will show how to directly use AMG in solving (3). We first show some denotations to discretize (3) (cf. [19]). We partition the domain Ω=(0,1)×(0,1) into 𝐿×𝐿 uniform cells with mesh size ℎ=1/𝐿. The cell centers are 𝑥_𝑙,𝑦_𝑘=1𝑙−21ℎ,𝑘−2ℎ,𝑙,𝑘=1,…,𝐿. (5) Following the ideas of [12, 19], we discretize (3) by standard five-point finite difference scheme to get −𝛼ℎ𝑐_{𝑙+1/2,𝑘}𝑢_𝑙+1,𝑘−𝑢_𝑙,𝑘ℎ−𝑐_{𝑙−1/2,𝑘}𝑢_𝑙,𝑘−𝑢_{𝑙−1,𝑘}ℎ−𝛼ℎ𝑐_{𝑙,𝑘+1/2}𝑢_𝑙,𝑘+1−𝑢_𝑙,𝑘ℎ−𝑐_{𝑙,𝑘−1/2}𝑢_𝑙,𝑘−𝑢_{𝑙,𝑘−1}ℎ+𝐾^∗(𝐾𝑢−𝑧)_𝑙,𝑘=0,𝑙,𝑘=1,…,𝐿, (6) with homogeneous Neumann boundary condition: 𝑢_0,𝑘=𝑢_1,𝑘,𝑢_𝐿+1,𝑘=𝑢_𝐿,𝑘,𝑢_𝑙,0=𝑢_𝑙,1,𝑢_𝑙,𝐿+1=𝑢_𝑙,𝐿, (7) where 𝑐_{𝑙+1/2,𝑘}=1𝐷_𝑥𝑢_{𝑙+1/2,𝑘}²+𝐷_𝑦𝑢_{𝑙+1/2,𝑘}²,+𝛽 (8) with 𝐷_𝑥𝑢_{𝑙+1/2,𝑘}=𝑢_𝑙+1,𝑘−𝑢_𝑙,𝑘ℎ,𝐷_𝑦𝑢_{𝑙+1/2,𝑘}=12𝑢_{𝑙+1,𝑘+1}−𝑢_{𝑙+1,𝑘−1}+𝑢2ℎ_𝑙,𝑘+1−𝑢_{𝑙,𝑘−1},2ℎ (9) and so forth.

To simplify the notation, we abbreviate (6) as 𝛼𝐿(𝑢)𝑢=−𝐾^∗(𝐾𝑢−𝑧), (10) where 𝐿(𝑣)𝑤=−∇_ℎ⋅∇_ℎ𝑤∇_ℎ𝑣²+𝛽, (11) which is fully nonlinear with wildly varying coefficient.

Consider the general system (10) in an 𝐿²×𝐿²: 𝐴𝑈=𝐹, (12) with the matrix 𝐴=𝛼𝐿(𝑢) in (10). In general, 𝐴 varies wildly near areas of high contrast of the image and need not be diagonally dominant. Nevertheless, 𝐴 is symmetric and positive definite. Moreover, the matrix 𝐾^∗𝐾 is wide-banded, and the spectra of the matrices 𝐿(𝑢) and 𝐾^∗𝐾 are quite differently distributed. It is difficult to compute (10). In [13], Vogel and Oman adopted the combination of a fixed point iteration and a product PCG to handle the nonlinear term and the linear system, respectively. In [15], Chan et al. proposed another preconditioner based on the fast cosine transform.

Algebraic multigrid method is considered in [19, 24, 25] as a linear solver. In contrast to geometric multigrid method, the construction of the coarse grids Ω^𝑚 is solely based on the algebraic information provided by the matrix 𝐴 in an AMG setting. The details and effectiveness can be found, for example, in [23, 24, 26]. For the restriction and coarse grid operators, a simple and popular approach is the Galerkin type construction. Namely, 𝐼_𝑚^𝑚+1=(𝐼^𝑚_𝑚+1)^∗ and 𝐴^𝑚+1=𝐼_𝑚^𝑚+1𝐴^𝑚𝐼^𝑚_𝑚+1. How to construct the interpolation operators 𝐼^𝑚_𝑚+1 is the most important part in the method. The convergence proof for this improved AMG method was given in [25] when 𝐴^𝑚 is symmetric positive definite. The robustness of these interpolation formulae is supported by plenty of numerical evidence [24, 25]. The detailed algorithm using the AMG method to solve (10) is in Algorithm 1.

Initialize: 𝑢0=𝑧 While 𝑢𝑘−𝑢𝑘−12𝑢𝑘2>𝜖,    Compute 𝑢 by 𝛼𝐿(𝑢𝑘)𝑢𝑘+1=−𝐾∗(𝐾𝑢𝑘−𝑧). end

Algorithm 1 

2.1. The Split Bregman Method for Isotropic TV Model

In [22], the split Bregman method is applied to the ROF model for image denoising problem. For completeness, we show the main progress about the split Bregman method applied to the ROF model. Firstly, consider discretizing ROF model into _Ω||𝜇∇𝑢𝑑𝑥𝑑𝑦+2‖𝐾𝑢−𝑧‖²_𝐿2≈_𝑖,𝑗∇_𝑥𝑢²_𝑖,𝑗+∇_𝑦𝑢²_𝑖,𝑗+𝜇2‖𝐾𝑢−𝑧‖²_𝐿2. (13) Setting 𝑑_𝑥≈∇_𝑥𝑢 and 𝑑_𝑦≈∇_𝑦𝑢, the split Bregman formulation of the isotropic problem then becomes min_{𝑢,𝑑𝑥,𝑑𝑦}𝑑_𝑥,𝑑_𝑦2+𝜇2‖𝐾𝑢−𝑧‖²₂+𝜆2𝑑_𝑥−∇_𝑥𝑢−𝑏_𝑥²₂+𝜆2𝑑_𝑦−∇_𝑦𝑢−𝑏_𝑦²₂, (14) where 𝑑_𝑥,𝑑_𝑦2=_𝑖,𝑗𝑑²_{𝑥,𝑖,𝑗}+𝑑²_{𝑦,𝑖,𝑗}. (15) Consider it into the following three subproblems: 𝑢^𝑘+1=min_𝑢𝜇2‖𝐾𝑢−𝑧‖²₂+𝜆2𝑑^𝑘_𝑥−∇_𝑥𝑢−𝑏^𝑘_𝑥²₂+𝜆2𝑑^𝑘_𝑦−∇_𝑦𝑢−𝑏^𝑘_𝑦²₂,𝑑_𝑥^𝑘+1,𝑑_𝑦^𝑘+1=min_{𝑑𝑥,𝑑𝑦}𝑑_𝑥,𝑑_𝑦2+𝜆2𝑑^𝑘_𝑥−∇_𝑥𝑢^𝑘+1−𝑏^𝑘_𝑥²₂+𝜆2𝑑^𝑘_𝑦−∇_𝑦𝑢^𝑘+1−𝑏^𝑘_𝑦²₂,𝑏_𝑥^𝑘+1=𝑏^𝑘_𝑥+∇_𝑥𝑢^𝑘+1−𝑑_𝑥^𝑘+1,𝑏_𝑦^𝑘+1=𝑏^𝑘_𝑦+∇_𝑦𝑢^𝑘+1−𝑑_𝑦^𝑘+1. (16)(17)(18)

The formula (16) is 𝑙₂-norm; it can be solved using the variation approach. The corresponding Euler-Lagrange equation for (16) is 𝜇𝐾^∗𝐾−𝜆Δ𝑢^𝑘+1=𝜇𝐾^∗𝑍+𝜆∇^𝑇_𝑥𝑑^𝑘_𝑥−𝑏^𝑘_𝑥+𝜆∇^𝑇_𝑦𝑑^𝑘_𝑦−𝑏^𝑘_𝑦. (19) The difference scheme is 𝜇𝐾^∗𝐾𝑢^𝑘+1_𝑖𝑗+𝜆4𝑢^𝑘+1_𝑖𝑗−𝑢^𝑘+1_𝑖+1,𝑗−𝑢^𝑘+1_{𝑖−1,𝑗}−𝑢^𝑘+1_𝑖,𝑗+1−𝑢^𝑘+1_{𝑖,𝑗−1}=𝐹^𝑘∶=𝜇𝐾^∗𝑧_𝑖𝑗+𝜆𝑑^𝑘_𝑥−𝑏^𝑘_{𝑥𝑖−1,𝑗}−𝑑^𝑘_𝑥−𝑏^𝑘_{𝑥𝑖,𝑗}+𝑑^𝑘_𝑦−𝑏^𝑘_{𝑦𝑖,𝑗−1}−𝑑^𝑘_𝑦−𝑏^𝑘_{𝑦𝑖,𝑗}. (20) When 𝐾=𝐼, the difference scheme of the previous formula is simplified as 𝜇+4𝜆𝑢^𝑘+1_𝑖𝑗−𝜆𝑢^𝑘+1_𝑖+1,𝑗+𝑢^𝑘+1_{𝑖−1,𝑗}+𝑢^𝑘+1_𝑖,𝑗+1+𝑢^𝑘+1_{𝑖,𝑗−1}=𝜇𝑧_𝑖𝑗+𝜆𝑑^𝑘_{𝑥,𝑖−1,𝑗}−𝑑^𝑘_{𝑥,𝑖,𝑗}+𝑑^𝑘_{𝑦,𝑖,𝑗−1}−𝑑^𝑘_{𝑦,𝑖,𝑗}−𝑏^𝑘_{𝑥,𝑖−1,𝑗}+𝑏^𝑘_{𝑥,𝑖,𝑗}−𝑏^𝑘_{𝑦,𝑖,𝑗−1}+𝑏^𝑘_{𝑦,𝑖,𝑗}. (21) Using Gauss-Seidel method for (21), we have 𝑢^𝑘+1_𝑖𝑗=𝐺^𝑘_𝑖,𝑗𝜆∶=𝑢𝜇+4𝜆^{𝑘𝑖+1,𝑗}+𝑢^𝑘+1_{𝑖−1,𝑗}+𝑢^𝑘_𝑖,𝑗+1+𝑢^𝑘+1_{𝑖,𝑗−1}+𝑑^𝑘_{𝑥,𝑖−1,𝑗}−𝑑^𝑘_{𝑥,𝑖,𝑗}+𝑑^𝑘_{𝑦,𝑖,𝑗−1}−𝑑^𝑘_{𝑦,𝑖,𝑗}−𝑏^𝑘_{𝑥,𝑖−1,𝑗}+𝑏^𝑘_{𝑥,𝑖,𝑗}−𝑏^𝑘_{𝑦,𝑖,𝑗−1}+𝑏^𝑘_{𝑦,𝑖,𝑗}+𝜇𝑧𝜇+4𝜆^𝑖𝑗. (22)

Problem (17) can be solved using the fast shrinkage operators: 𝑑_𝑥^𝑘+1=shrink∇_𝑥𝑢^𝑘+𝑏^𝑘_𝑥,1𝜆=max𝑠^𝑘−1𝜆⋅∇,0^𝑥𝑢^𝑘+𝑏^𝑘_𝑥𝑠^𝑘,𝑑_𝑦^𝑘+1=shrink∇_𝑦𝑢^𝑘+𝑏^𝑘_𝑦,1𝜆=max𝑠^𝑘−1𝜆⋅∇,0^𝑦𝑢^𝑘+𝑏^𝑘_𝑦𝑠^𝑘,𝑠^𝑘=∇_𝑥𝑢^𝑘+𝑏^𝑘_𝑥²+∇_𝑦𝑢^𝑘+𝑏^𝑘_𝑦². (23)

The optimal value of 𝑑 can be explicitly gotten using the above shrinkage operators. The formula (18) can be discretized directly. Thus the split Bregman algorithm of the isotropic TV denoising is as in Algorithm 2.

Initialize: 𝑢0=𝑧,𝑑0𝑥=0,𝑑0𝑦=0,𝑏0𝑥=0,𝑏0𝑦=0. While 𝑢𝑘−𝑢𝑘−12𝑢𝑘2>𝜖,    𝑢𝑘+1=𝐺𝑘 by (22),    𝑑𝑥𝑘+1=max(𝑠𝑘−1𝜆∇,0)⋅𝑥𝑢𝑘+1+𝑏𝑘𝑥𝑠𝑘,    𝑑𝑦𝑘+1=max(𝑠𝑘−1𝜆∇,0)⋅𝑦𝑢𝑘+1+𝑏𝑘𝑦𝑠𝑘,    𝑏𝑥𝑘+1=𝑏𝑘𝑥+(∇𝑥𝑢𝑘+1−𝑑𝑥𝑘+1),    𝑏𝑦𝑘+1=𝑏𝑘𝑦+(∇𝑦𝑢𝑘+1−𝑑𝑦𝑘+1). End

Algorithm 2 

2.2. The Split Bregman Method Based on AMG for Isotropic TV Model

It is easy to find that the accuracy of the split Bregman method is not high because the Gauss-Seidel iteration of 𝑢 is executed only once per loop. At the same time, the AMG algorithm has fast convergence and high precision. It can be applied to large signal-to-noise ratio of image denoising, which has strong stability. So, we consider the combination of the two methods, and get a new denoising algorithm which is the split Bregman method based on AMG (here we call it as Algorithm 3). This algorithm contains four steps.

Step 1: Decomposing the denoising model of ROF to three subproblems (16)–(18) using the split    Bregman method. Step 2: To get 𝑢 by solving (22) using the AMG method. Step 3: To get 𝑑𝑥,𝑑𝑦, using the shrinkage operator (23). Step 4: To directly discretize 𝑏𝑥,𝑏𝑦, using difference scheme.

Algorithm 3 

Moreover, we use Krylov technique to accelerate subspace extrapolation.

2.3. The Split Bregman Method for Anisotropic TV Model

The anisotropic total variation of 𝑢 [27] is represented by ∇_𝑥𝑢₁+∇_𝑦𝑢₁. (24) The anisotropic TV model for denoising can be formulated as the following minimization problem: min_{𝑢,𝑑𝑥,𝑑𝑦}∇_𝑥𝑢₁+∇_𝑦𝑢₁+𝜇2‖𝑢−𝑧‖²₂, (25) where 𝜇 is an appropriately chosen positive parameter. Define an operator [27] 1cut𝑐,𝜆|𝑐|−1=𝑐−max𝜆⋅𝑐,0|𝑐|=1𝜆1,𝑐>𝜆,1𝑐,−𝜆1≤𝑐≤𝜆,−1𝜆1,𝑐<−𝜆. (26) With the previous operator, we have Algorithm 4 to solve the anisotropic TV denoising problem.

Initialize: 𝑢1=𝑧,𝑏0𝑥=0,𝑏0𝑦=0. While 𝑢𝑘−𝑢𝑘−12𝑢𝑘2<𝜖,    𝑏𝑘𝑥=cut(∇𝑥𝑢𝑘+𝑏𝑥𝑘−1,1𝜆),    𝑏𝑘𝑦=cut(∇𝑦𝑢𝑘+𝑏𝑦𝑘−1,1𝜆),    𝑢𝑘+1𝜆=𝑧−𝜇(∇𝑇𝑥𝑏𝑘+∇𝑇𝑦𝑏𝑘). End

Algorithm 4 

It should be noticed that there is no need to solve any equations like the isotropic model. This algorithm need not solve any elliptic equation of 𝑢 in [27], so the procedure is rapid, simple, and not costly.

3. Convergence Analysis

We consider the discretized system with general noise and blur: 𝜇𝐾^∗𝐾𝑢^𝑘+1_𝑖𝑗+𝜆4𝑢^𝑘+1_𝑖𝑗−𝑢^𝑘+1_𝑖+1,𝑗−𝑢^𝑘+1_{𝑖−1,𝑗}−𝑢^𝑘+1_𝑖,𝑗+1−𝑢^𝑘+1_{𝑖,𝑗−1}=𝜇𝐾^∗𝑧_𝑖𝑗+𝜆𝑑^𝑘_𝑥−𝑏^𝑘_{𝑥𝑖−1,𝑗}−𝑑^𝑘_𝑥−𝑏^𝑘_{𝑥𝑖,𝑗}+𝑑^𝑘_𝑦−𝑏^𝑘_{𝑦𝑖,𝑗−1}−𝑑^𝑘_𝑦−𝑏^𝑘_{𝑦𝑖,𝑗}=𝐹_𝑖𝑗. (27) To simplify the notation, we abbreviate (27) as −𝜆Δ𝑢=−𝜇𝐾^∗𝐾𝑢+𝐹. (28) In general, the matrix 𝐾^∗𝐾 would be dense. There is one way to completely avoid the matrix inversion. A natural approach would then be −𝜆Δ𝑢^(𝑠+1)=−𝜇𝐾^∗𝐾𝑢^(𝑠)+𝐹. (29) It turns out that the iteration (29) is not robust and may diverge even for weak blur operator corresponding to the mask (1/64)(1,1,4,1,1)^𝑇(1,1,4,1,1) with 𝜎=10. To overcome this instability, we propose to add a linear term 𝐵+𝜂𝐼 on both sides of (29) and get −𝜆Δ+𝐵+𝜂𝐼𝑢^(𝑠+1)=−𝜇𝐾^∗𝐾𝑢^(𝑠)−𝐵−𝜂𝐼𝑢^(𝑠)+𝐹, (30) where 𝜂 is arbitrary small positive number and 𝐵 is a symmetric and positive definite matrix to be determined through the following formal analysis. Now we will analyze the convergence property of (29).

Theorem 1. If the matrix 𝐵 satisfies 2𝐵−𝜇𝐾^∗𝐾≥0, (31) the equation (30) is locally convergent.

Let 𝑢^∞ be the solution of (29); that is, −𝜆Δ𝑢^(∞)=−𝜇𝐾^∗𝐾𝑢^(∞)+𝐹. (32) We then define 𝑒^(∞)=𝑢^(𝑠)−𝑢^(∞), and the error equation is −𝜆Δ+𝐵+𝜂𝐼𝑒^(𝑠+1)=−𝜇𝐾^∗𝐾𝑢^(𝑠)−𝐵−𝜂𝐼𝑒^(𝑠), (33) which is equivalent to 𝜂𝑒^(𝑠+1)−𝑒^(𝑠)=−𝑃𝑒^(𝑠+1)−𝑄𝑒^(𝑠) (34) with 𝑃=(−𝜆Δ+𝐵),𝑄=𝐵−𝜇𝐾^∗𝐾. (35) It is easy to show that the matrices 𝑃 and 𝑄 are symmetric. Taking the inner product on both sides of (34) with 𝑃𝑒^(𝑠+1)−𝑄𝑒^(𝑠) gives 𝜂𝑒^(𝑠+1)−𝑒^(𝑠),𝑃𝑒^(𝑠+1)−𝑄𝑒^(𝑠)=−𝑃𝑒^(𝑠+1)−𝑄𝑒^(𝑠)2₂≤0, (36) and we will now rearrange the formula (36) as 𝜂𝑒^(𝑠+1)−𝑒^(𝑠)𝑇×𝑃−𝑄2𝑒^(𝑠+1)+𝑒^(𝑠)+𝑃+𝑄2𝑒^(𝑠+1)−𝑒^(𝑠)≤0. (37) Rewrite the previous formula as 𝑒^(𝑠+1)𝑇(𝑃−𝑄)𝑒^(𝑠+1)−𝑒^(𝑠)𝑇(𝑃−𝑄)𝑒^(𝑠)+𝑒^(𝑠+1)−𝑒^(𝑠)𝑇(𝑃+𝑄)𝑒^(𝑠+1)−𝑒^(𝑠)≤0. (38) If (𝑃−𝑄) and (𝑃+𝑄) are symmetric and positive definite, it follows from (38) that 𝑒^(𝑠+1)2_𝑃−𝑄−𝑒^(𝑠)2_𝑃−𝑄+𝑒^(𝑠+1)−𝑒^(𝑠)2_𝑃+𝑄≤0, (39) where ‖𝑣‖²_𝐴=𝑣^𝑇𝐴𝑣 with positive definite matrix 𝐴. The behavior of (39) leads to 𝑒^(𝑠+1)2_𝑃−𝑄≤𝑒^(𝑠)2_𝑃−𝑄. (40) In other words, the iteration (30) is locally convergent if both 𝑃−𝑄=−𝜆Δ+𝜇𝐾^∗𝐾,𝑃+𝑄=−𝜆Δ+2𝐵−𝜇𝐾^∗𝐾 (41) are symmetric and positive definite which means 𝑃−𝑄=−𝜆Δ+𝜇𝐾^∗𝐾>0,𝑃+𝑄=−𝜆Δ+2𝐵−𝜇𝐾^∗𝐾>0. (42) Since −𝜆Δ and 𝜇𝐾^∗𝐾 are symmetric and positive definite, a sufficient condition for (42) is given by (31).

Remark. We can see that (31) indeed is a good guideline for devising the iteration of 𝑢. We will elaborate on the choice of 𝐵 better suited for numerical purposes in Section 4. Indeed, the iteration (30) may fail to converge when (31) is not satisfied. See Section 4.2 for more details. In addition, the preprocessing of initial data proposed in [19] is no longer required when (31) is satisfied.

4. Numerical Experiments and Discussion

4.1. Images and Blurring Operators

The numerical examples given in this paper are based on the three images shown in Figure 1. The first one is a satellite image (Image I), the second one is made up of simple shapes (Image II) used in the literature [15], and the last one is Lena image (Image III). We have experiments on restoring the three images blurred by the following three different blurring operators [23].

(a)

(b)

(c)

Figure 1 

Original images. Image I, Image II, and Image III.

Blur operator I: 1.5()7841,2,3,16,3,2,1^𝑇()=1,2,3,16,3,2,11.5.784123163212463264236948963163248256483216369489632463264212316321 (43)

Blur operator II: an out-of-focus blur with the kernel of convolution given by the kernel 𝑣_𝑘,𝑙=𝑐13,if𝑘²+𝑙²≤4,0,otherwise, (44) where 𝑐=1.2188.

Blur operator III: a truncated Gaussian blur given by 𝑣_𝑘,𝑙=𝑐𝑒^{−𝜏(𝑘2+𝑙2)}if|𝑘|,|𝑙|≤5,0,otherwise, (45) where 𝑐=0.176 and 𝜏=0.38.

4.2. The Linear Stabilizing Term

The condition (31) serves as a general guideline for the choice of the matrix. In this section, we will consider a quinta-diagonal matrix 𝐵 to be compared with a diagonal matrix 𝐵 [23].

4.2.1. Diagonal Matrix 𝐵

In order to solve (30) efficiently, an obvious candidate is diagonal matrices of the form 𝐵=𝜀𝐼. (46) The condition (31) then reads 𝜀≥𝜀^∗=12𝜆_max𝜇𝐾^∗𝐾=12𝜆_max𝜇𝐾². (47)

To proceed with the estimate of 𝜆_max(𝐾), we notice that our sample blurring operators I–III can be represented by a mask of the form 𝑣𝑉=^{−𝑟,−𝑟}𝑣_{−𝑟,−𝑟+1}⋯𝑣_−𝑟,𝑟𝑣_{−𝑟+1,−𝑟}𝑣_{−𝑟+1,−𝑟+1}⋯𝑣_{−𝑟+1,𝑟}𝑣⋮⋮⋮⋮_𝑟,−𝑟𝑣_{𝑟,−𝑟+1}⋯𝑣_{𝑟,𝑟(2𝑟+1)∗(2𝑟+1)}, (48) with 𝑣_𝑘,𝑙=𝑣_𝑙,𝑘=𝑣_−𝑘,𝑙=𝑣_𝑘,−𝑙. It is easy to see that, under the Neumann boundary conditions, the largest eigenvalue for such 𝐾 is given by 𝜆_max(𝐾)=^𝐿2_𝑖=1𝐾_𝑖,𝑗=^𝑟_{𝑘,𝑙=−𝑟}𝑣_𝑘,𝑙, (49) where the second term of (49) is independent and the corresponding eigenfunction is a constant. In summary, we have the following sufficient condition for (31): 𝐵=𝜀𝐼,𝜀≥𝜀^∗=12𝜇^𝑟_{𝑘,𝑙=−𝑟}𝑣_𝑘,𝑙². (50)

Similar conditions can be derived for variants of (50), such as 1𝐵=2diag𝐾^∗𝐾+12𝛿𝐼,𝛿≥𝛿^∗=max_{𝑖𝑗≠𝑖}𝜇𝐾^∗𝐾_𝑖,𝑗=𝜇^𝑟_{𝑘,𝑙=−𝑟}𝑣_𝑘,𝑙²−^𝑟_{𝑘,𝑙=−𝑟}𝑣²_𝑘,𝑙,𝛿≥𝛿^∗=max_{𝑖𝑗≠𝑖}𝜇𝐾^∗𝐾_𝑖,𝑗−𝐾^∗𝐾_𝑖,𝑖=𝜇^𝑟_{𝑘,𝑙=−𝑟}𝑣_𝑘,𝑙²−^𝑟_{𝑘,𝑙=−𝑟}𝑣²_𝑘,𝑙, (51) or 𝐵=diag𝐾^∗𝐾+𝛾𝐼,𝛾≥𝛾^∗=12max_{𝑖𝑗≠𝑖}𝜇𝐾^∗𝐾_𝑖,𝑗−2𝐾^∗𝐾_𝑖,𝑖=12𝜇^𝑟_{𝑘,𝑙=−𝑟}𝑣_𝑘,𝑙²−^𝑟_{𝑘,𝑙=−𝑟}𝑣²_𝑘,𝑙. (52)

In our examples, the diagonal entries of 𝐾^∗𝐾 are constant-valued except for the near-boundary pixels. We expect the performances of (50), (51), and (52) to be comparable to each other. In this paper, most of our numerical experiments are conducted using (52) with 𝛾=1. See Section 5 for more details.

4.2.2. Quinta-Diagonal Matrix 𝐵

In addition to diagonal matrices, we have also explored quinta-diagonal matrices 𝐵_𝑎,𝑏 given by the mask 0𝑏0𝑏𝑎𝑏0𝑏0. (53)

𝐵_𝑎,𝑏 has the same support as laplace operator Δ; therefore the computational cost in the AMG step is comparable to diagonal 𝐵. In view of (33) and (31), it is clear that the optimal 𝐵_(𝑎,𝑏)^∗ is the one that solves min_𝑎,𝑏𝜌−𝜆Δ+𝐵_(𝑎,𝑏)⁻¹𝐵_(𝑎,𝑏)−𝜇𝐾^∗𝐾subjectto2𝐵_(𝑎,𝑏)−𝜇𝐾^∗𝐾≥0, (54) where 𝜌(𝐴) is the spectral radius of 𝐴. The actual minimizer of (54) depends on the true image 𝑢^(∞) and there is no simple way of finding it. An accessible approximation is given by the following modified minimization problem: min_𝑎,𝑏𝜌𝐵_(𝑎,𝑏)−𝜇𝐾^∗𝐾subjectto2𝐵_(𝑎,𝑏)−𝜇𝐾^∗𝐾≥0. (55) In our examples, the common eigenbasis of 𝐵_𝑎,𝑏 and 𝐾 is given by 𝑒_𝑚,𝑛𝑥_𝑘,𝑦_𝑙=cos𝑚𝜋𝑥_𝑘cos𝑛𝜋𝑦_𝑙,𝑚,𝑛=1,…,𝐿. (56) It is straightforward to compute the corresponding eigenvalues for (𝐵_(𝑎,𝑏)−𝜇𝐾^∗𝐾) and the numerical solution of (55) by varying 𝑎 and 𝑏 over a suitable range. The optimal (𝑎,𝑏)^∗, together with the critical values of 𝜀, 𝛿, and 𝛾, is given in Table 1. It is not clear a priori why (55) should result in better performance. Nevertheless, we find from our experiment that 𝐵_(𝑎,𝑏)^∗ indeed performs no worse than diagonal 𝐵 and can be significantly faster in some cases. See more details in Section 5.

Table 1 

Critical values of the parameters, 𝜀^∗,𝛿^∗,𝛾^∗ and the solution to (55).

4.3. The Stopping Criterion and Acceleration Technique

The outer iteration is stopped upon a relative decrease of the (normalized) residual by a factor of 10⁻⁴ for the blurring operators, namely, when the condition 𝑟^(𝑁)_𝑙2𝑟⁽¹⁾_𝑙2≤10⁻⁴ (57) is reached. Here the normalization of the residual is given by 𝑟^(𝑠+1)=𝐷^(𝑠+1)−1×(−𝜆Δ+𝐵)𝑣^(𝑠+1)+𝜇𝐾^∗𝐾−𝐵𝑣^(𝑠+1)−𝐹, (58) where 𝑣^(𝑠+1) is an approximate solution of (28) obtained with one 𝑉-cycle iteration (see also Section 5) and 𝐷^(𝑠+1)=diag(−𝜆Δ+𝐵). Note that the larger entries of the normalizing factor (𝐷^(𝑠+1))⁻¹ correspond to regions where 𝑣^(𝑠+1) is less smooth. It amplifies the (unnormalized) residual on regions where 𝑣^(𝑠+1) either has a jump or requires further denoising, therefore doing a better job measuring the image quality. Numerical experiments confirmed that the normalized residual is indeed a good indicator for the quality of denoising and deblurring.

To accelerate and stabilize the outer iteration, we have incorporated the Krylov subspace method [19, 28, 29] in our implementation. The approximate solution is optimized every 𝑝 steps on a subspace of dimension 𝑀≤𝑝. To be more precise, we take two fixed integers 0≤𝑀≤𝑝, and for any integer 𝑛>0, let 𝑢𝑐₁,…,𝑐_𝑀=𝑢^(4𝑛)+^𝑀_𝑚=1𝑐_𝑚𝑢^{(4𝑛+1−𝑚)}−𝑢^{(4𝑛−𝑚)}. (59) The residual of 𝑢(𝑐₁,…,𝑐_𝑀) can be approximated by 𝑟̃𝑐₁,…,𝑐=_𝑀^def𝑟^(4𝑛)+^𝑀_𝑚=1𝑐_𝑚𝑟^{(4𝑛+1−𝑚)}−𝑟^{(4𝑛−𝑚)}. (60)

One can minimize 𝑟(̃𝑐₁,…,𝑐_𝑀) with respect to (𝑐₁,…,𝑐_𝑀) to get min_{𝑐1,…,𝑐𝑀}𝑟̃𝑐₁,…,𝑐_𝑀𝑙²=𝑟̃𝑐^∗₁,…,𝑐^∗_𝑀𝑙² (61) and reset 𝑢^(𝑝𝑛) to 𝑢(𝑐^∗₁,…,𝑐^∗_𝑀) before going to (𝑝𝑛+1) the outer iteration. The numerical results demonstrate that the Krylov acceleration method is an efficient way to accelerate the outer iteration. See Section 5.

5. Numerical Results

To generate an observed image 𝑧, a Gaussian white noise with mean 0 and variance 𝜎 is added to the blurred image. The level of noise can be measured by the signal-to-noise ratio (SNR): ‖‖SNR=𝐾𝑢−𝑧_𝑙²‖𝑢‖_𝑙². (62) The signal-to-blur ratio (SBR) ‖‖SBR=𝐾𝑢−𝑢_𝑙²‖𝑢‖_𝑙² (63) measures the strength of the blurring operator 𝐾. The two kinds of error 𝑒2=‖𝑢^𝑘+1−𝑢^𝑘‖_𝐿² and 𝑒0=‖𝑢^𝑘+1−𝑢^𝑘‖_𝐿^∞ are used to judge the iterative error as the iteration ends. The cost of the iterative process is denoted by “time." These indexes are captioned in each example for readers’ convenience.

In all our examples, the iteration begins with the observed image 𝑢⁽⁰⁾=𝑧 and all images are contaminated with small noise 𝜎=10 and large noise 𝜎=40. We choose 𝐵 according to (52) with 𝛾=1.0 except in Table 6 where we demonstrate the necessity of the linear stabilizing term and the overall efficiency of our scheme by varying 𝛾. It is enough to apply one single 𝑉-cycle in the AMG step with the Gauss-Seidel iteration as the smoother. We apply the Krylov acceleration procedure every 4 steps with the Krylov dimension 𝑀=2. The rate of convergence of the outer iteration can be significantly improved by the Krylov subspace technique.

5.1. Image Denoising

Table 2 shows the contrast results of Image II with the split Bregman method for isotropic TV model (Algorithm 2), the split Bregman method based on AMG (Algorithm 3), and the split Bregman method for anisotropic TV model (Algorithm 4). Figure 2 shows the denoising image for Image II with noise contamination of 𝜎=10 (SNR = 14.90%). Table 3 shows the contrast results of Image III containing 512∗512 pixels. Figure 3 shows the denoising image for Image III (512∗512) with noise contamination 𝜎=40 (SNR = 26.725%). It is clear that our algorithm (Algorithm 3) is efficient and suitable for elliptic equation.

Table 2 

Results of Figure 2.

Table 3 

Results of Figure 3.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

(a) Original image. (b) Noisy image contaminated by 𝜎=10 (SNR =14.90%). (c) Denoised image with 3 iterations of Algorithm 2. (d) Denoised image with 50 iterations of Algorithm 2. (e) Denoised image with 50 iterations of Algorithm 4. (f) Denoised image with 3 iterations of Algorithm 3.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

(a) Original image. (b) Noisy image contaminated by 𝜎=40 (SNR = 26.725%). (c) Denoised image with 3 iterations of Algorithm 2. (d) Denoised image with 50 iterations of Algorithm 2. (e) Denoised image with 50 iterations of Algorithm 4. (f) Denoised image with 3 iterations of Algorithm 3.

5.2. Image Denoising and Deblurring

The compared results are given in Tables 4-5. Figure 4 shows the denoising and deblurring images for Image I contaminated by noise with 𝜎=40 and blurring operator III (SNR=51.34%,SBR=45.31%). Figure 5 shows the denoising and deblurring images for Image II contaminated by noise with 𝜎=10 and blurring operator I (SNR=5.96%,SBR=48.79%).

Table 4 

Results of Figure 4.

Table 5 

Results of Figure 5.

Table 6 

Number of outer nonlinear iterations 𝑁 for various 𝛾. Here 𝛼=0.05, 𝜎=29.72 for Image I and 𝛼=2.5, 𝜎=89.16 for Image II.

(a)

(b)

(c)

(d)

Figure 4 

(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 20 iterations of Algorithm 1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.

(a)

(b)

(c)

(d)

Figure 5 

(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 10 iterations of Algorithm 1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.

We continue our numerical experiments by varying the parameter 𝛾 and compare their performance. The number of iterations of a typical example is shown in Table 6. The combination of 𝛾≈𝛾^∗ and 𝑀=2 gives the optimal result among possible choices of 𝛾 and 𝑀 in general. The number of iterations is insensitive to the variation of 𝛾 near 𝛾^∗ and gradually grows as 𝛾 increases. For simplicity of presentation, we take 𝛾=1 and 𝑀=2 in all other numerical examples. We have also included the results for 𝐵=𝐵_(𝑎,𝑏)^∗ in Table 6 for comparison. In general, 𝐵=𝐵_(𝑎,𝑏)^∗ performs no worse than (52) with 𝛾=𝛾^∗. In some cases, 𝐵=𝐵_(𝑎,𝑏)^∗ can be up to 25% faster. In addition, if the blur operator is mild, our method can be used to solve 𝐾𝑢=𝑧 directly.

According to the previous figures and tables, we get some conclusions here. First of all, in the second part of (42), we have neglected the −𝜆Δ term to obtain a sufficient condition (31), from which we further derived the critical value 𝛾^∗ in (52). We therefore expect 𝛾^∗ to be slightly overestimated. In other words, the outer iteration should converge for 𝛾≥𝛾^∗, but not necessarily for 𝛾<𝛾^∗. This is in accordance with our numerical result. See the row with 𝑀=0 in Table 6. Secondly, the Krylov subspace technique helps to stabilize the outer iteration for 𝛾<𝛾^∗, but eventually fails if 𝛾 becomes too small. This result demonstrates the necessity of the linear stabilizing term.

6. Conclusion

In this paper, we propose a new algorithm by adding a linear term for the total variation-based image denoising and deblurring which combine the split Bregman method, the algebraic multigrid method, and Krylov subspace acceleration. Through formal convergence analysis, we derived an explicit formula for a linear stabilizing term. Numerical experiments demonstrate that our algorithm is efficient and robust for deblurring and denoising problems.

Acknowledgments

The authors thank the referees for very useful remarks and suggestions. The research is supported by NSFC (no. 10801049 and no. 11271126) and Foundation of North China Electric Power University.