# Efficient Algorithm for Isotropic and Anisotropic Total Variation Deblurring and Denoising

Authors: Yuying Shi; Qianshun Chang
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 1. Introduction

The purpose of image restoration is to recover original imageu from an observed data z (noisy and blurred image) from the relation
(1)z=Ku+n,
where K is a known linear blurring operator and u 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]:
(2)minu(α∫Ω‍|∇u|2+βdxdy+12∥Ku-z∥L22).
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
(3)-α∇·(∇u|∇u|2+β)+K*(Ku-z)=0,
where K* is the adjoint operator of K with respect to standard L2 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
(4)TV2(u)=∫Ω‍(|∇xu|+|∇yu|)dxdy.
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 O(k2) is maintained by GAPG, where k 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 V-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 K. 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 K. 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 Section2, 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 L×L uniform cells with mesh size h=1/L. The cell centers are
(5)(xl,yk)=((l-12)h,(k-12)h),l,k=1,…,L.
Following the ideas of [12, 19], we discretize (3) by standard five-point finite difference scheme to get
(6)-αh(cl+1/2,kul+1,k-ul,kh-cl-1/2,kul,k-ul-1,kh)-αh(cl,k+1/2ul,k+1-ul,kh-cl,k-1/2ul,k-ul,k-1h)+(K*(Ku-z))l,k=0,l,k=1,…,L,
with homogeneous Neumann boundary condition:
(7)u0,k=u1,k,uL+1,k=uL,k,ul,0=ul,1,ul,L+1=ul,L,
where
(8)cl+1/2,k=1((Dxu)l+1/2,k)2+((Dyu)l+1/2,k)2+β,
with
(9)(Dxu)l+1/2,k=ul+1,k-ul,kh,(Dyu)l+1/2,k=12(ul+1,k+1-ul+1,k-12h+ul,k+1-ul,k-12h),
and so forth.To simplify the notation, we abbreviate (6) as
(10)αL(u)u=-K*(Ku-z),
where
(11)L(v)w=-∇h·(∇hw|∇hv|2+β),
which is fully nonlinear with wildly varying coefficient.Consider the general system (10) in an L2×L2:
(12)AU=F,
with the matrix A=αL(u) in (10). In general, A varies wildly near areas of high contrast of the image and need not be diagonally dominant. Nevertheless, A is symmetric and positive definite. Moreover, the matrix K*K is wide-banded, and the spectra of the matrices L(u) and K*K 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 Ωm is solely based on the algebraic information provided by the matrix A 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, Imm+1=(Im+1m)* and Am+1=Imm+1AmIm+1m. How to construct the interpolation operators Im+1m is the most important part in the method. The convergence proof for this improved AMG method was given in [25] when Am 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.Algorithm 1
Initialize:u0=z
While∥uk-uk-1∥2∥uk∥2>ϵ,
Computeu by αL(uk)uk+1=-K*(Kuk-z).
end


### 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
(13)∫Ω‍|∇u|dxdy+μ2∥Ku-z∥L22≈∑i,j‍(∇xu)i,j2+(∇yu)i,j2+μ2∥Ku-z∥L22.
Setting dx≈∇xu and dy≈∇yu, the split Bregman formulation of the isotropic problem then becomes
(14)minu,dx,dy∥(dx,dy)∥2+μ2∥Ku-z∥22+λ2∥dx-∇xu-bx∥22+λ2∥dy-∇yu-by∥22,
where
(15)∥(dx,dy)∥2=∑i,j‍dx,i,j2+dy,i,j2.
Consider it into the following three subproblems:
(16)uk+1=minuμ2∥Ku-z∥22+λ2∥dxk-∇xu-bxk∥22+λ2∥dyk-∇yu-byk∥22,(17)(dxk+1,dyk+1)=mindx,dy∥(dx,dy)∥2+λ2∥dxk-∇xuk+1-bxk∥22+λ2∥dyk-∇yuk+1-byk∥22,(18)bxk+1=bxk+(∇xuk+1-dxk+1),byk+1=byk+(∇yuk+1-dyk+1).The formula (16) is l2-norm; it can be solved using the variation approach. The corresponding Euler-Lagrange equation for (16) is
(19)(μK*K-λΔ)uk+1=μK*Z+λ∇xT(dxk-bxk)+λ∇yT(dyk-byk).
The difference scheme is
(20)μ(K*Kuk+1)ij+λ(4uijk+1-ui+1,jk+1-ui-1,jk+1-ui,j+1k+1-ui,j-1k+1)=Fk:=μ(K*z)ij+λ[(dxk-bxk)i-1,j-(dxk-bxk)i,j+(dyk-byk)i,j-1-(dyk-byk)i,j].
When K=I, the difference scheme of the previous formula is simplified as
(21)(μ+4λ)uijk+1-λ(ui+1,jk+1+ui-1,jk+1+ui,j+1k+1+ui,j-1k+1)=μzij+λ(dx,i-1,jk-dx,i,jk+dy,i,j-1k-dy,i,jk-bx,i-1,jk+bx,i,jk-by,i,j-1k+by,i,jk).
Using Gauss-Seidel method for (21), we have
(22)uijk+1=Gi,jk:=λμ+4λ(ui+1,jk+ui-1,jk+1+ui,j+1k+ui,j-1k+1+dx,i-1,jk-dx,i,jk+dy,i,j-1k-dy,i,jk-bx,i-1,jk+bx,i,jk-by,i,j-1k+by,i,jk)+μμ+4λzij.Problem (17) can be solved using the fast shrinkage operators:
(23)dxk+1=shrink(∇xuk+bxk,1λ)=max(sk-1λ,0)·∇xuk+bxksk,dyk+1=shrink(∇yuk+byk,1λ)=max(sk-1λ,0)·∇yuk+byksk,sk=|∇xuk+bxk|2+|∇yuk+byk|2.The optimal value ofd 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.Algorithm 2
Initialize:u0=z,dx0=0,dy0=0,bx0=0,by0=0.
While∥uk-uk-1∥2∥uk∥2>ϵ,
uk+1=Gk by (22),
dxk+1=max(sk-1λ,0)·∇xuk+1+bxksk,
dyk+1=max(sk-1λ,0)·∇yuk+1+byksk,
bxk+1=bxk+(∇xuk+1-dxk+1),
byk+1=byk+(∇yuk+1-dyk+1).
End


### 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 ofu 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.Algorithm 3
Step1: Decomposing the denoising model of ROF to three subproblems (16)–(18) using the split
Bregman method.
Step2: To get u by solving (22) using the AMG method.
Step3: To get dx,dy, using the shrinkage operator (23).
Step4: To directly discretize bx,by, using difference scheme.Moreover, we use Krylov technique to accelerate subspace extrapolation.


### 2.3. The Split Bregman Method for Anisotropic TV Model

The anisotropic total variation ofu [27] is represented by
(24)∥∇xu∥1+∥∇yu∥1.
The anisotropic TV model for denoising can be formulated as the following minimization problem:
(25)minu,dx,dy∥∇xu∥1+∥∇yu∥1+μ2∥u-z∥22,
where μ is an appropriately chosen positive parameter. Define an operator [27]
(26)cut(c,1λ)=c-max(|c|-1λ,0)·c|c|={1λ,c>1λ,c,-1λ≤c≤1λ,-1λ,c<-1λ.
With the previous operator, we have Algorithm  4 to solve the anisotropic TV denoising problem.Algorithm 4
Initialize:u1=z,bx0=0,by0=0.
While∥uk-uk-1∥2∥uk∥2<ϵ,
bxk=cut(∇xuk+bxk-1,1λ),
byk=cut(∇yuk+byk-1,1λ),
uk+1=z-λμ(∇xTbk+∇yTbk).
EndIt should be noticed that there is no need to solve any equations like the isotropic model. This algorithm need not solve any elliptic equation ofu in [27], so the procedure is rapid, simple, and not costly.


## 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
(13)∫Ω‍|∇u|dxdy+μ2∥Ku-z∥L22≈∑i,j‍(∇xu)i,j2+(∇yu)i,j2+μ2∥Ku-z∥L22.
Setting dx≈∇xu and dy≈∇yu, the split Bregman formulation of the isotropic problem then becomes
(14)minu,dx,dy∥(dx,dy)∥2+μ2∥Ku-z∥22+λ2∥dx-∇xu-bx∥22+λ2∥dy-∇yu-by∥22,
where
(15)∥(dx,dy)∥2=∑i,j‍dx,i,j2+dy,i,j2.
Consider it into the following three subproblems:
(16)uk+1=minuμ2∥Ku-z∥22+λ2∥dxk-∇xu-bxk∥22+λ2∥dyk-∇yu-byk∥22,(17)(dxk+1,dyk+1)=mindx,dy∥(dx,dy)∥2+λ2∥dxk-∇xuk+1-bxk∥22+λ2∥dyk-∇yuk+1-byk∥22,(18)bxk+1=bxk+(∇xuk+1-dxk+1),byk+1=byk+(∇yuk+1-dyk+1).The formula (16) is l2-norm; it can be solved using the variation approach. The corresponding Euler-Lagrange equation for (16) is
(19)(μK*K-λΔ)uk+1=μK*Z+λ∇xT(dxk-bxk)+λ∇yT(dyk-byk).
The difference scheme is
(20)μ(K*Kuk+1)ij+λ(4uijk+1-ui+1,jk+1-ui-1,jk+1-ui,j+1k+1-ui,j-1k+1)=Fk:=μ(K*z)ij+λ[(dxk-bxk)i-1,j-(dxk-bxk)i,j+(dyk-byk)i,j-1-(dyk-byk)i,j].
When K=I, the difference scheme of the previous formula is simplified as
(21)(μ+4λ)uijk+1-λ(ui+1,jk+1+ui-1,jk+1+ui,j+1k+1+ui,j-1k+1)=μzij+λ(dx,i-1,jk-dx,i,jk+dy,i,j-1k-dy,i,jk-bx,i-1,jk+bx,i,jk-by,i,j-1k+by,i,jk).
Using Gauss-Seidel method for (21), we have
(22)uijk+1=Gi,jk:=λμ+4λ(ui+1,jk+ui-1,jk+1+ui,j+1k+ui,j-1k+1+dx,i-1,jk-dx,i,jk+dy,i,j-1k-dy,i,jk-bx,i-1,jk+bx,i,jk-by,i,j-1k+by,i,jk)+μμ+4λzij.Problem (17) can be solved using the fast shrinkage operators:
(23)dxk+1=shrink(∇xuk+bxk,1λ)=max(sk-1λ,0)·∇xuk+bxksk,dyk+1=shrink(∇yuk+byk,1λ)=max(sk-1λ,0)·∇yuk+byksk,sk=|∇xuk+bxk|2+|∇yuk+byk|2.The optimal value ofd 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.Algorithm 2
Initialize:u0=z,dx0=0,dy0=0,bx0=0,by0=0.
While∥uk-uk-1∥2∥uk∥2>ϵ,
uk+1=Gk by (22),
dxk+1=max(sk-1λ,0)·∇xuk+1+bxksk,
dyk+1=max(sk-1λ,0)·∇yuk+1+byksk,
bxk+1=bxk+(∇xuk+1-dxk+1),
byk+1=byk+(∇yuk+1-dyk+1).
End


## 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 ofu 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.Algorithm 3
Step1: Decomposing the denoising model of ROF to three subproblems (16)–(18) using the split
Bregman method.
Step2: To get u by solving (22) using the AMG method.
Step3: To get dx,dy, using the shrinkage operator (23).
Step4: To directly discretize bx,by, using difference scheme.Moreover, we use Krylov technique to accelerate subspace extrapolation.


## 2.3. The Split Bregman Method for Anisotropic TV Model

The anisotropic total variation ofu [27] is represented by
(24)∥∇xu∥1+∥∇yu∥1.
The anisotropic TV model for denoising can be formulated as the following minimization problem:
(25)minu,dx,dy∥∇xu∥1+∥∇yu∥1+μ2∥u-z∥22,
where μ is an appropriately chosen positive parameter. Define an operator [27]
(26)cut(c,1λ)=c-max(|c|-1λ,0)·c|c|={1λ,c>1λ,c,-1λ≤c≤1λ,-1λ,c<-1λ.
With the previous operator, we have Algorithm  4 to solve the anisotropic TV denoising problem.Algorithm 4
Initialize:u1=z,bx0=0,by0=0.
While∥uk-uk-1∥2∥uk∥2<ϵ,
bxk=cut(∇xuk+bxk-1,1λ),
byk=cut(∇yuk+byk-1,1λ),
uk+1=z-λμ(∇xTbk+∇yTbk).
EndIt should be noticed that there is no need to solve any equations like the isotropic model. This algorithm need not solve any elliptic equation ofu in [27], so the procedure is rapid, simple, and not costly.


## 3. Convergence Analysis

We consider the discretized system with general noise and blur:(27)μ(K*Kuk+1)ij+λ(4uijk+1-ui+1,jk+1-ui-1,jk+1-ui,j+1k+1-ui,j-1k+1)=μ(K*z)ij+λ[(dxk-bxk)i-1,j-(dxk-bxk)i,j=(μK*z)ij+λ+(dyk-byk)i,j-1-(dyk-byk)i,j]=Fij.
To simplify the notation, we abbreviate (27) as
(28)-λΔu=-μK*Ku+F. In general, the matrix K*K would be dense. There is one way to completely avoid the matrix inversion. A natural approach would then be
(29)-λΔu(s+1)=-μK*Ku(s)+F.
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)T(1,1,4,1,1) with σ=10. To overcome this instability, we propose to add a linear term B+ηI on both sides of (29) and get
(30)(-λΔ+B+ηI)u(s+1)=-(μK*Ku(s)-B-ηI)u(s)+F,
where η is arbitrary small positive number and B 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 matrixB satisfies
(31)2B-μK*K≥0,
the equation (30) is locally convergent.Letu∞ be the solution of (29); that is,
(32)-λΔu(∞)=-μK*Ku(∞)+F.
We then define e(∞)=u(s)-u(∞), and the error equation is
(33)(-λΔ+B+ηI)e(s+1)=-(μK*Ku(s)-B-ηI)e(s),
which is equivalent to
(34)η(e(s+1)-e(s))=-(Pe(s+1)-Qe(s))
with
(35)P=(-λΔ+B),Q=(B-μK*K).
It is easy to show that the matrices P and Q are symmetric. Taking the inner product on both sides of (34) with Pe(s+1)-Qe(s) gives
(36)η(e(s+1)-e(s),Pe(s+1)-Qe(s))=-∥Pe(s+1)-Qe(s)∥22≤0,
and we will now rearrange the formula (36) as
(37)η(e(s+1)-e(s))T×(P-Q2(e(s+1)+e(s))+P+Q2(e(s+1)-e(s)))≤0.
Rewrite the previous formula as
(38)(e(s+1))T(P-Q)e(s+1)-(e(s))T(P-Q)e(s)+(e(s+1)-e(s))T(P+Q)(e(s+1)-e(s))≤0.
If (P-Q) and (P+Q) are symmetric and positive definite, it follows from (38) that
(39)∥(e(s+1))∥P-Q2-∥(e(s))∥P-Q2+∥e(s+1)-e(s)∥P+Q2≤0,
where ∥v∥A2=vTAv with positive definite matrix A. The behavior of (39) leads to
(40)∥(e(s+1))∥P-Q2≤∥(e(s))∥P-Q2.
In other words, the iteration (30) is locally convergent if both
(41)P-Q=-λΔ+μK*K,P+Q=-λΔ+2B-μK*K
are symmetric and positive definite which means
(42)P-Q=-λΔ+μK*K>0,P+Q=-λΔ+2B-μK*K>0.
Since -λΔ and μK*K 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 u. We will elaborate on the choice of B 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 Figure1. 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].Original images. Image I, Image II, and Image III.
(a)
(b)
(c)Blur operator I:(43)1.5784(1,2,3,16,3,2,1)T(1,2,3,16,3,2,1)=1.5784[123163212463264236948963163248256483216369489632463264212316321].Blur operator II: an out-of-focus blur with the kernel of convolution given by the kernel(44)vk,l={c13,ifk2+l2≤4,0,otherwise,
where c=1.2188.Blur operator III: a truncated Gaussian blur given by(45)vk,l={ce-τ(k2+l2)if|k|,|l|≤5,0,otherwise,
where c=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 B to be compared with a diagonal matrix B [23].


#### 4.2.1. Diagonal MatrixB

In order to solve (30) efficiently, an obvious candidate is diagonal matrices of the form
(46)B=εI.
The condition (31) then reads
(47)ε≥ε*=12λmax(μK*K)=12λmax(μK2).To proceed with the estimate ofλmax(K), we notice that our sample blurring operators I–III can be represented by a mask of the form
(48)V=[v-r,-rv-r,-r+1⋯v-r,rv-r+1,-rv-r+1,-r+1⋯v-r+1,r⋮⋮⋮⋮vr,-rvr,-r+1⋯vr,r](2r+1)*(2r+1),
with vk,l=vl,k=v-k,l=vk,-l. It is easy to see that, under the Neumann boundary conditions, the largest eigenvalue for such K is given by
(49)λmax(K)=∑i=1L2‍Ki,j=∑k,l=-rr‍vk,l,
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):
(50)B=εI,ε≥ε*=12μ(∑k,l=-rr‍vk,l)2.Similar conditions can be derived for variants of (50), such as
(51)B=12diag(K*K)+12δI,δ≥δ*=maxi∑j≠i‍μ(K*K)i,j=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,δ≥δ*=maxi(∑j≠i‍μ(K*K)i,j-(K*K)i,i)=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,
or
(52)B=diag(K*K)+γI,γ≥γ*=12maxi(∑j≠i‍μ(K*K)i,j-2(K*K)i,i)=12μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2.In our examples, the diagonal entries ofK*K 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 MatrixB

In addition to diagonal matrices, we have also explored quinta-diagonal matricesBa,b given by the mask
(53)[0b0bab0b0].B
a
,
b has the same support as laplace operator Δ; therefore the computational cost in the AMG step is comparable to diagonal B. In view of (33) and (31), it is clear that the optimal B(a,b)* is the one that solves
(54)mina,bρ((-λΔ+B(a,b))-1(B(a,b)-μK*K))subjectto2B(a,b)-μK*K≥0,
where ρ(A) is the spectral radius of A. The actual minimizer of (54) depends on the true image u(∞) and there is no simple way of finding it. An accessible approximation is given by the following modified minimization problem:
(55)mina,bρ(B(a,b)-μK*K)subjectto2B(a,b)-μK*K≥0.
In our examples, the common eigenbasis of Ba,b and K is given by
(56)em,n(xk,yl)=cos(mπxk)cos(nπyl),m,n=1,…,L.
It is straightforward to compute the corresponding eigenvalues for (B(a,b)-μK*K) and the numerical solution of (55) by varying a and b over a suitable range. The optimal (a,b)*, 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 B(a,b)* indeed performs no worse than diagonal B 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).
Blur operator
b
*
δ
*
γ
*
(
a
,
b
)
*
I
1.13
1.95
0.83
(0.72, 0.17)
II
0.74
1.37
0.63
(0.49, 0.09)
III
1.06
1.98
0.93
(0.65, 0.15)


### 4.3. The Stopping Criterion and Acceleration Technique

The outer iteration is stopped upon a relative decrease of the (normalized) residual by a factor of10-4 for the blurring operators, namely, when the condition
(57)∥r(N)∥l2∥r(1)∥l2≤10-4
is reached. Here the normalization of the residual is given by
(58)r(s+1)=(D(s+1))-1×((-λΔ+B)v(s+1)+(μK*K-B)v(s+1)-F),
where v(s+1) is an approximate solution of (28) obtained with one V-cycle iteration (see also Section 5) and D(s+1)=diag(-λΔ+B). Note that the larger entries of the normalizing factor (D(s+1))-1 correspond to regions where v(s+1) is less smooth. It amplifies the (unnormalized) residual on regions where v(s+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 p steps on a subspace of dimension M≤p. To be more precise, we take two fixed integers 0≤M≤p, and for any integer n>0, let
(59)u~(c1,…,cM)=u(4n)+∑m=1M‍cm(u(4n+1-m)-u(4n-m)).
The residual of u~(c1,…,cM) can be approximated by
(60)r~(c1,…,cM)=defr(4n)+∑m=1M‍cm(r(4n+1-m)-r(4n-m)).One can minimizer~(c1,…,cM) with respect to (c1,…,cM) to get
(61)minc1,…,cM∥r~(c1,…,cM)∥l2=∥r~(c1*,…,cM*)∥l2
and reset u(pn) to u~(c1*,…,cM*) before going to (pn+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.


## 4.1. Images and Blurring Operators

The numerical examples given in this paper are based on the three images shown in Figure1. 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].Original images. Image I, Image II, and Image III.
(a)
(b)
(c)Blur operator I:(43)1.5784(1,2,3,16,3,2,1)T(1,2,3,16,3,2,1)=1.5784[123163212463264236948963163248256483216369489632463264212316321].Blur operator II: an out-of-focus blur with the kernel of convolution given by the kernel(44)vk,l={c13,ifk2+l2≤4,0,otherwise,
where c=1.2188.Blur operator III: a truncated Gaussian blur given by(45)vk,l={ce-τ(k2+l2)if|k|,|l|≤5,0,otherwise,
where c=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 B to be compared with a diagonal matrix B [23].


### 4.2.1. Diagonal MatrixB

In order to solve (30) efficiently, an obvious candidate is diagonal matrices of the form
(46)B=εI.
The condition (31) then reads
(47)ε≥ε*=12λmax(μK*K)=12λmax(μK2).To proceed with the estimate ofλmax(K), we notice that our sample blurring operators I–III can be represented by a mask of the form
(48)V=[v-r,-rv-r,-r+1⋯v-r,rv-r+1,-rv-r+1,-r+1⋯v-r+1,r⋮⋮⋮⋮vr,-rvr,-r+1⋯vr,r](2r+1)*(2r+1),
with vk,l=vl,k=v-k,l=vk,-l. It is easy to see that, under the Neumann boundary conditions, the largest eigenvalue for such K is given by
(49)λmax(K)=∑i=1L2‍Ki,j=∑k,l=-rr‍vk,l,
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):
(50)B=εI,ε≥ε*=12μ(∑k,l=-rr‍vk,l)2.Similar conditions can be derived for variants of (50), such as
(51)B=12diag(K*K)+12δI,δ≥δ*=maxi∑j≠i‍μ(K*K)i,j=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,δ≥δ*=maxi(∑j≠i‍μ(K*K)i,j-(K*K)i,i)=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,
or
(52)B=diag(K*K)+γI,γ≥γ*=12maxi(∑j≠i‍μ(K*K)i,j-2(K*K)i,i)=12μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2.In our examples, the diagonal entries ofK*K 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 MatrixB

In addition to diagonal matrices, we have also explored quinta-diagonal matricesBa,b given by the mask
(53)[0b0bab0b0].B
a
,
b has the same support as laplace operator Δ; therefore the computational cost in the AMG step is comparable to diagonal B. In view of (33) and (31), it is clear that the optimal B(a,b)* is the one that solves
(54)mina,bρ((-λΔ+B(a,b))-1(B(a,b)-μK*K))subjectto2B(a,b)-μK*K≥0,
where ρ(A) is the spectral radius of A. The actual minimizer of (54) depends on the true image u(∞) and there is no simple way of finding it. An accessible approximation is given by the following modified minimization problem:
(55)mina,bρ(B(a,b)-μK*K)subjectto2B(a,b)-μK*K≥0.
In our examples, the common eigenbasis of Ba,b and K is given by
(56)em,n(xk,yl)=cos(mπxk)cos(nπyl),m,n=1,…,L.
It is straightforward to compute the corresponding eigenvalues for (B(a,b)-μK*K) and the numerical solution of (55) by varying a and b over a suitable range. The optimal (a,b)*, 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 B(a,b)* indeed performs no worse than diagonal B 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).
Blur operator
b
*
δ
*
γ
*
(
a
,
b
)
*
I
1.13
1.95
0.83
(0.72, 0.17)
II
0.74
1.37
0.63
(0.49, 0.09)
III
1.06
1.98
0.93
(0.65, 0.15)


## 4.2.1. Diagonal MatrixB

In order to solve (30) efficiently, an obvious candidate is diagonal matrices of the form
(46)B=εI.
The condition (31) then reads
(47)ε≥ε*=12λmax(μK*K)=12λmax(μK2).To proceed with the estimate ofλmax(K), we notice that our sample blurring operators I–III can be represented by a mask of the form
(48)V=[v-r,-rv-r,-r+1⋯v-r,rv-r+1,-rv-r+1,-r+1⋯v-r+1,r⋮⋮⋮⋮vr,-rvr,-r+1⋯vr,r](2r+1)*(2r+1),
with vk,l=vl,k=v-k,l=vk,-l. It is easy to see that, under the Neumann boundary conditions, the largest eigenvalue for such K is given by
(49)λmax(K)=∑i=1L2‍Ki,j=∑k,l=-rr‍vk,l,
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):
(50)B=εI,ε≥ε*=12μ(∑k,l=-rr‍vk,l)2.Similar conditions can be derived for variants of (50), such as
(51)B=12diag(K*K)+12δI,δ≥δ*=maxi∑j≠i‍μ(K*K)i,j=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,δ≥δ*=maxi(∑j≠i‍μ(K*K)i,j-(K*K)i,i)=μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2,
or
(52)B=diag(K*K)+γI,γ≥γ*=12maxi(∑j≠i‍μ(K*K)i,j-2(K*K)i,i)=12μ(∑k,l=-rr‍vk,l)2-∑k,l=-rr‍vk,l2.In our examples, the diagonal entries ofK*K 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 MatrixB

In addition to diagonal matrices, we have also explored quinta-diagonal matricesBa,b given by the mask
(53)[0b0bab0b0].B
a
,
b has the same support as laplace operator Δ; therefore the computational cost in the AMG step is comparable to diagonal B. In view of (33) and (31), it is clear that the optimal B(a,b)* is the one that solves
(54)mina,bρ((-λΔ+B(a,b))-1(B(a,b)-μK*K))subjectto2B(a,b)-μK*K≥0,
where ρ(A) is the spectral radius of A. The actual minimizer of (54) depends on the true image u(∞) and there is no simple way of finding it. An accessible approximation is given by the following modified minimization problem:
(55)mina,bρ(B(a,b)-μK*K)subjectto2B(a,b)-μK*K≥0.
In our examples, the common eigenbasis of Ba,b and K is given by
(56)em,n(xk,yl)=cos(mπxk)cos(nπyl),m,n=1,…,L.
It is straightforward to compute the corresponding eigenvalues for (B(a,b)-μK*K) and the numerical solution of (55) by varying a and b over a suitable range. The optimal (a,b)*, 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 B(a,b)* indeed performs no worse than diagonal B 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).
Blur operator
b
*
δ
*
γ
*
(
a
,
b
)
*
I
1.13
1.95
0.83
(0.72, 0.17)
II
0.74
1.37
0.63
(0.49, 0.09)
III
1.06
1.98
0.93
(0.65, 0.15)


## 4.3. The Stopping Criterion and Acceleration Technique

The outer iteration is stopped upon a relative decrease of the (normalized) residual by a factor of10-4 for the blurring operators, namely, when the condition
(57)∥r(N)∥l2∥r(1)∥l2≤10-4
is reached. Here the normalization of the residual is given by
(58)r(s+1)=(D(s+1))-1×((-λΔ+B)v(s+1)+(μK*K-B)v(s+1)-F),
where v(s+1) is an approximate solution of (28) obtained with one V-cycle iteration (see also Section 5) and D(s+1)=diag(-λΔ+B). Note that the larger entries of the normalizing factor (D(s+1))-1 correspond to regions where v(s+1) is less smooth. It amplifies the (unnormalized) residual on regions where v(s+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 p steps on a subspace of dimension M≤p. To be more precise, we take two fixed integers 0≤M≤p, and for any integer n>0, let
(59)u~(c1,…,cM)=u(4n)+∑m=1M‍cm(u(4n+1-m)-u(4n-m)).
The residual of u~(c1,…,cM) can be approximated by
(60)r~(c1,…,cM)=defr(4n)+∑m=1M‍cm(r(4n+1-m)-r(4n-m)).One can minimizer~(c1,…,cM) with respect to (c1,…,cM) to get
(61)minc1,…,cM∥r~(c1,…,cM)∥l2=∥r~(c1*,…,cM*)∥l2
and reset u(pn) to u~(c1*,…,cM*) before going to (pn+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 imagez, 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):
(62)SNR=∥Ku-z∥l2∥u∥l2.
The signal-to-blur ratio (SBR)
(63)SBR=∥Ku-u∥l2∥u∥l2
measures the strength of the blurring operator K. The two kinds of error e2=∥uk+1-uk∥L2 and e0=∥uk+1-uk∥L∞ 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 imageu(0)=z and all images are contaminated with small noise σ=10 and large noise σ=40. We choose B 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 V-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 M=2. The rate of convergence of the outer iteration can be significantly improved by the Krylov subspace technique.


### 5.1. Image Denoising

Table2 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 Figure2.
Algorithm
(
σ
,
μ
,
λ
)
Iteration step
e
2
e
0
Time
(
10,0.1,0.2
)
N
=
3
6.88
87.36
0.0625
Algorithm2
N
=
50
1.16
30.64
0.9218
(
40,0.05,0.1
)
N
=
3
11.98
93.46
0.0468
N
=
50
10.04
106.16
0.6562
Algorithm4
(
10,0.1,0.2
)
N
=
50
4.88
46.25
0.5156
(
40,0.1,0.2
)
N
=
50
19.24
156.89
0.5156
Algorithm3
(
10,0.2,0.2
)
N
=
3
4.40
73.23
0.3125
(
40,0.03,0.06
)
N
=
3
9.29
111.18
0.3125Table 3
Results of Figure3.
Algorithm
(
σ
,
μ
,
λ
)
Iteration step
e
2
e
0
Time
(
10,0.1,0.2
)
N
=
3
5.94
133.24
0.2500
Algorithm2
N
=
50
5.96
139.42
3.7187
(
40,0.05,0.1
)
N
=
3
12.04
146.55
0.1875
N
=
50
12.12
147.62
2.750
Algorithm4
(
10,0.1,0.2
)
N
=
50
8.27
56.69
2.0625
(
40,0.03,0.00375
)
N
=
50
19.71
139.02
2.0781
Algorithm3
(
10,0.2,0.2
)
N
=
2
5.64
112.07
1.156
(
40,0.03,0.06
)
N
=
2
10.51
109.00
1.156(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)(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.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.2. Image Denoising and Deblurring

The compared results are given in Tables4-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 Figure4.
Algorithm
(
σ
,
ite
,
α
)/(σ,ite,μ,λ)
Blur operator
e
0
e
2
Time
I
93.91
8.47
5.8906
(10, 10, 0.02)
II
155.95
12.22
5.1250
Algorithm1
III
149.45
12.09
8.2500
I
167.84
14.67
7.2500
(40, 20, 0.15)
II
167.72
15.61
6.8906
III
168.29
15.08
11.1093
I
161.52
13.34
1.7500
(10, 5, 0.04, 0.08)
II
172.19
15.61
1.2812
Algorithm3
III
172.62
14.96
2.9062
I
151.96
15.17
1.7343
(40, 5, 0.04, 0.08)
II
170.75
16.24
1.2968
III
170.94
15.62
2.9531Table 5
Results of Figure5.
Algorithm
(
σ
,
ite
,
α
)/(σ,ite,μ,λ)
Blur operator
e
0
e
2
Time
I
82.67
3.30
5.4375
(10, 10, 0.05)
II
104.79
4.36
5.2187
Algorithm1
III
104.68
4.99
8.3437
I
122.00
5.86
5.4843
(40, 20, 1.0)
II
126.02
6.99
6.9375
III
123.50
6.61
11.0468
I
100.84
6.51
1.7500
(10, 5, 0.04, 0.08)
II
113.03
8.29
1.2812
Algorithm3
III
111.03
8.01
3.0156
I
109.85
9.79
1.7656
(40, 5, 0.04, 0.08)
II
117.62
9.40
1.2812
III
117.85
9.09
2.9375Table 6
Number of outer nonlinear iterationsN for various γ. Here α=0.05, σ=29.72 for Image I and α=2.5, σ=89.16 for Image II.
Blur
Image
M
0.2
0.4
0.6
0.8
0.9
1.2
2
4
10
B
a
*
,
b
*
0
∞
∞
∞
∞
53
60
81
136
297
46
I
1
∞
∞
148
50
48
28
32
48
93
23
III
2
∞
∞
37
28
28
27
29
41
78
22
0
∞
∞
109
109
109
113
119
131
162
109
II
1
∞
39
37
36
36
37
40
44
65
36
2
∞
34
32
32
33
33
35
42
57
32(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 20 iterations of Algorithm1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.
(a)
(b)
(c)
(d)(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 10 iterations of Algorithm1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.
(a)
(b)
(c)
(d)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 M=2 gives the optimal result among possible choices of γ and M 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 M=2 in all other numerical examples. We have also included the results for B=B(a,b)* in Table 6 for comparison. In general, B=B(a,b)* performs no worse than (52) with γ=γ*. In some cases, B=B(a,b)* can be up to 25% faster. In addition, if the blur operator is mild, our method can be used to solve Ku=z 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 M=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.


## 5.1. Image Denoising

Table2 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 Figure2.
Algorithm
(
σ
,
μ
,
λ
)
Iteration step
e
2
e
0
Time
(
10,0.1,0.2
)
N
=
3
6.88
87.36
0.0625
Algorithm2
N
=
50
1.16
30.64
0.9218
(
40,0.05,0.1
)
N
=
3
11.98
93.46
0.0468
N
=
50
10.04
106.16
0.6562
Algorithm4
(
10,0.1,0.2
)
N
=
50
4.88
46.25
0.5156
(
40,0.1,0.2
)
N
=
50
19.24
156.89
0.5156
Algorithm3
(
10,0.2,0.2
)
N
=
3
4.40
73.23
0.3125
(
40,0.03,0.06
)
N
=
3
9.29
111.18
0.3125Table 3
Results of Figure3.
Algorithm
(
σ
,
μ
,
λ
)
Iteration step
e
2
e
0
Time
(
10,0.1,0.2
)
N
=
3
5.94
133.24
0.2500
Algorithm2
N
=
50
5.96
139.42
3.7187
(
40,0.05,0.1
)
N
=
3
12.04
146.55
0.1875
N
=
50
12.12
147.62
2.750
Algorithm4
(
10,0.1,0.2
)
N
=
50
8.27
56.69
2.0625
(
40,0.03,0.00375
)
N
=
50
19.71
139.02
2.0781
Algorithm3
(
10,0.2,0.2
)
N
=
2
5.64
112.07
1.156
(
40,0.03,0.06
)
N
=
2
10.51
109.00
1.156(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)(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.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.2. Image Denoising and Deblurring

The compared results are given in Tables4-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 Figure4.
Algorithm
(
σ
,
ite
,
α
)/(σ,ite,μ,λ)
Blur operator
e
0
e
2
Time
I
93.91
8.47
5.8906
(10, 10, 0.02)
II
155.95
12.22
5.1250
Algorithm1
III
149.45
12.09
8.2500
I
167.84
14.67
7.2500
(40, 20, 0.15)
II
167.72
15.61
6.8906
III
168.29
15.08
11.1093
I
161.52
13.34
1.7500
(10, 5, 0.04, 0.08)
II
172.19
15.61
1.2812
Algorithm3
III
172.62
14.96
2.9062
I
151.96
15.17
1.7343
(40, 5, 0.04, 0.08)
II
170.75
16.24
1.2968
III
170.94
15.62
2.9531Table 5
Results of Figure5.
Algorithm
(
σ
,
ite
,
α
)/(σ,ite,μ,λ)
Blur operator
e
0
e
2
Time
I
82.67
3.30
5.4375
(10, 10, 0.05)
II
104.79
4.36
5.2187
Algorithm1
III
104.68
4.99
8.3437
I
122.00
5.86
5.4843
(40, 20, 1.0)
II
126.02
6.99
6.9375
III
123.50
6.61
11.0468
I
100.84
6.51
1.7500
(10, 5, 0.04, 0.08)
II
113.03
8.29
1.2812
Algorithm3
III
111.03
8.01
3.0156
I
109.85
9.79
1.7656
(40, 5, 0.04, 0.08)
II
117.62
9.40
1.2812
III
117.85
9.09
2.9375Table 6
Number of outer nonlinear iterationsN for various γ. Here α=0.05, σ=29.72 for Image I and α=2.5, σ=89.16 for Image II.
Blur
Image
M
0.2
0.4
0.6
0.8
0.9
1.2
2
4
10
B
a
*
,
b
*
0
∞
∞
∞
∞
53
60
81
136
297
46
I
1
∞
∞
148
50
48
28
32
48
93
23
III
2
∞
∞
37
28
28
27
29
41
78
22
0
∞
∞
109
109
109
113
119
131
162
109
II
1
∞
39
37
36
36
37
40
44
65
36
2
∞
34
32
32
33
33
35
42
57
32(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 20 iterations of Algorithm1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.
(a)
(b)
(c)
(d)(a) Original image. (b) Blurred and noisy image. (c) Denoised and deblurred image with 10 iterations of Algorithm1. (d) Denoised and deblurred image with 5 iterations of Algorithm 3.
(a)
(b)
(c)
(d)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 M=2 gives the optimal result among possible choices of γ and M 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 M=2 in all other numerical examples. We have also included the results for B=B(a,b)* in Table 6 for comparison. In general, B=B(a,b)* performs no worse than (52) with γ=γ*. In some cases, B=B(a,b)* can be up to 25% faster. In addition, if the blur operator is mild, our method can be used to solve Ku=z 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 M=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.



---

*Source: 797239-2013-03-17.xml*