# An Adaptive Total Generalized Variation Model with Augmented Lagrangian Method for Image Denoising

Authors: Chuan He; Changhua Hu; Xiaogang Yang; Huafeng He; Qi Zhang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We propose an adaptive total generalized variation (TGV) based model, aiming at achieving a balance between edge preservation and region smoothness for image denoising. The variable splitting (VS) and the classical augmented Lagrangian method (ALM) are used to solve the proposed model. With the proposed adaptive model and ALM, the regularization parameter, which balances the data fidelity and the regularizer, is refreshed with a closed form in each iterate, and the image denoising can be accomplished without manual interference. Numerical results indicate that our method is effective in staircasing effect suppression and holds superiority over some other state-of-the-art methods both in quantitative and in qualitative assessment.


---

## Body

## 1. Introduction

In the past few decades, many variation or partial differential equation (PDE) based restoration models [1–7] have been proposed to recover images from degraded observations, due to the ability of preserving significant image features such as edges or textures. Among these models, the total variation (TV) model, also named the Rudin-Osher-Fatemi (ROF) model [1], is distinguished for excellent edge preserving ability and becomes one of the most widely used regularizers in image restoration [1, 2, 8–10]. In particular, the TV denoising problem is in the following form:
(1)argminu∫Ω|∇u|dx+λ2∫Ω|u-u0|2dx,
where Ω is an open bounded domain in two dimensions, u is the image to be restored, u0 is the observation containing Gaussian white noise, and λ is the regularization parameter which balances the regularization term and the data fidelity term. ∫Ω|∇u|dx is the TV seminorm of the bounded variation (BV) space BV(Ω). The TV model is highly effective in preserving edges and corners, compared with the quadratic Tikhonov model. However, only when the original image is piecewise constant, the TV model is proved to be optimal. In fact, staircasing effect usually appears because most of natural images are not piecewise constant. Staircasing effect cannot meet the demands of human vision, due to the new artificial edges which do not exist in original images.To overcome the drawback of the TV model, researchers suggest introducing the higher-order derivatives of image functions [3–7, 12–16]. In order to eliminate the staircasing effect of TV model, Chambolle and Lions [14] proposed the following infimal-convolution minimization functional:
(2)argminu,v∫Ω|∇u-∇v|dx+α∫Ω|∂2v|dxK+λ∫Ω|u-u0|2dx,
where discontinuous components of the image are allotted to u-v while regions of moderate slopes are assigned to v. The above model was proved to be practically efficient. Later, a modified form of (2) was proposed in [5] and its regularizer is of the following form:
(3)argminu,v∫Ω|∇u-∇v|dx+α∫Ω|Δv|dx.
That is, the second-order derivative in (2) is substituted by the Laplacian in (3). The similar use of Laplacian operator can also be seen in some PDE-based methods [3].Since the classical TV model could not distinguish jumps from smooth transitions, Chan et al. [12] considered an additional penalization of the discontinuities in images. Precisely, they adopt
(4)∫Ω|∇u|dx+α∫Ωψ(|∇u|)h(Δu)dx
as the regularization term, where ψ is a real-valued function whose value approaches 0 while |∇u| approaches infinity. The absence of the staircasing effect for this choice was verified in [15].Bredies et al. [17] proposed the concept of total generalized variation (TGV), which is considered to be the generalization of TV. The TGV model is defined as
(5)TGVαk(u)=sup⁡{∫Ωudivkvdx∣v∈Cck(Ω,Symk(Rd)),kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkikll∥divlv∥∞≤αl,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkl=0,…,k-1∫Ωudivkvdx}withSymk(Rd)={ζ:Rd×Rd⋯×Rd︸ktimes⟶Rkkkkkkkkkkkk|ζisk-linearandsymmetricζ:Rd×Rd⋯×Rd︸ktimes},
where d≥1 denotes the image dimension, and, throughout this paper, we assume d=2; Symk(Rd) is the space of symmetric k-tensors on Rd; Cck(Ω,Symk(Rd)) is the space of compactly supported symmetric tensor field; αl is fixed positive parameter. From the definition of TGVαk, we learn that it involves the derivatives of u of order one to k. When k=1 and α0=1, TGVαk degenerates to the classical TV. Thus TGV can be seen as a generalization of TV.TGV involves and balances higher-order derivatives ofu. Image reconstruction with TGV regularization usually leads to result with piecewise polynomial intensities and sharp edges. Therefore, TGV can effectively suppress the staircasing effect. In [17], an accelerated first-order method of Nesterov [18] was proposed to solve the TGV-regularized denoising problem.In this paper, we propose an adaptive second-order TGV-regularized model for denoising and derive an augmented Lagrangian approach to handle the suggested model. Our denoising model is as follows:(6)argminuTGVα2(u)s.t.∫Ω|u-u0|2dx≤c.
According to the standard Lagrange duality, for a given c, there exists a nonnegative λ such that
(7)argminuTGVα2(u)+λ2∫Ω|u-u0|2dx
is equivalent to (6). However, with (6), we can automatically estimate the regularization parameter λ. We first utilize an indicator function of the feasible set to transform problem (6) into an unconstrained one; then the variable splitting technique is applied to transform the resulting unconstrained problem into a problem with linear penalizing constraints; finally, the obtained constrained problem is solved by the alternating direction method of multipliers (ADMM) [19–22], which is an instance of the classical ALM. The resulting image denoising algorithm is effective in staircasing effect suppression compared with some TV-based denoising methods, due to the second-order TGV regularizer. Besides, it achieves the adaptive estimation of the regularization parameter without inner iterative scheme. It is worth noting that the idea of this paper can be extended to TGV models with higher order than two. However, for simplicity, we only treat the second-order model and this is adequate for a large class of natural images.Our method differs from the previous works on at least two aspects. On one hand, compared with [16], which adopted the accelerated first-order method of Nesterov [18] to handle the unconstrained TGV-based denoising problem (7), we apply ALM to the constrained TGV-based denoising problem (6) and achieve the automatic estimation of the regularization parameter λ. Our strategy avoids the extra cost on the manual selection of λ by try-and-error. On the other hand, compared with the existing TV-based adaptive methods [10, 23, 24], we propose a more complicated adaptive method based on TGV, and it is apt to achieve more attractive results than the TV-based methods.The outline of the rest of the paper is organized as follows. Section2 provides the description of the adaptive second-order TGV-based model for image denoising. Based on the Lagrange duality, an equivalent form of TGVα2 is suggested. The derivation of the proposed method is presented in Section 3. Section 4 gives the numerical results that demonstrate the effectiveness of the proposed method. At last, Section 5 ends this paper with a brief conclusion.


## 2. Adaptive Second-Order TGV-Based Model for Image Denoising

The space of bounded generalized variation (BGV) functions of orderk with weight α is defined as
(8)BGVαk(Ω)={u∈L1(Ω)|TGVαk(u)<∞}.
Correspondingly, the BGV norm is defined as
(9)∥u∥BGVαk=∥u∥1+TGVαk(u).
The TGV seminorm rather than the BGV norm is usually used as a regularizer.In this paper, we just takek=2 into consideration for simplicity. The second-order TGV can be written as
(10)TGVα2(u)=sup⁡{∫Ωudiv2vdx|v∈Cc2(Ω,Sym2(Rd)),kkkkkkkkkkkkkkkkkl∥v∥∞≤α0,∥div⁡v∥∞≤α1∫Ωudiv2vdx|v∈Cc2(Ω,Sym2(Rd)),},
where the divergences are defined as
(11)(div⁡v)i=∑j=1d∂vij∂xj,1≤i≤d;div2v=∑i,j=1d∂2vij∂xi∂xj.
In fact, Sym2(Rd) is equivalent to the space of all symmetric d×d matrices. The infinite norms in (10) are given by
(12)kklk∥v∥∞=sup⁡x∈Ω(∑i,j=1d|vij(x)|2)1/2,∥div⁡v∥∞=sup⁡x∈Ω(∑i=1d|(div⁡v)i|2)1/2.For the convenience of the derivation of our algorithm, we apply the discrete form in the following and the tensors and vectors are denoted in bold type font. In order to make use of ADMM, we apply an equivalent definition ofTGVα2 [17, 22] based on the Lagrange duality. With this definition, we have
(13)TGVα2(u)=min⁡pα0∥ε(p)∥1+α1∥∇u-p∥1,
where u∈Rmn denotes an m×n image, p∈Rmn×Rmn belongs to the two-dimensional 1-tensor field, and ε denotes the symmetrized derivative operator. Suppose that ui,j and pi,j denote the (i,j)th components of u and p, respectively. Then we have pi,j=[pi,j,1,pi,j,2]∈Sym1(R2) and the (i,j)th component of ε(p) is given by
(14)ε(p)i,j=[ε(p)i,j,1ε(p)i,j,3ε(p)i,j,3ε(p)i,j,2]=[∇x1pi,j,1(∇x2pi,j,1+∇x1pi,j,2)2(∇x2pi,j,1+∇x1pi,j,2)2∇x2pi,j,2]kkkkkkkkkk∈Sym2(R2)0≤i≤m,0≤j≤n,
where ∇x1 and ∇x2 denote the difference operators in directions x1 and x2. According to the definition of operators ∇ and ε, (∇u)i,j and ε(p)i,j are two-dimensional 1-tensor and symmetric 2-tensor, respectively. Besides, the ∥∥1s of p and ε(p) are defined as ∥p∥1=∑i,jm,n|pi,j|=∑i,jm,npi,j,12+pi,j,22 and ∥ε(p)∥1=∑i,jm,n|ε(p)i,j|=∑i,jm,nε(p)i,j,12+ε(p)i,j,22+2ε(p)i,j,32, respectively. The deduction of (13) is given in the Appendix.Then the constrained second-order TGV-regularized denoising problem (6) can be rewritten as
(15)argminu,pα0∥ε(p)∥1+α1∥∇u-p∥1s.t.kkkk∥u-u0∥22≤c.


## 3. Methodology

### 3.1. The Augmented Lagrangian Model of Adaptive TGV-Based Denoising

Problem (15) can be transformed into an unconstrained problem, with the following discontinuous objective functional:
(16)argminu,pα0∥ε(p)∥1+α1∥∇u-p∥1+IΦ(u),
where IΦ(u) is the indicator function of the feasible set defined by
(17)IΦ(u)={0,ifu∈Φ≜{u:∥u-u0∥22≤c};+∞otherwise.
Note that Φ is a closed Euclidean ball centered at u0 with radius c.The solution of problem (16) suffers from its nonlinearity and nondifferentiability. Referring to the variable splitting, we introduce three auxiliary variables to simplify the solution process of (16): a variable w for liberating u out from the constraint of the feasible set; a variable y and a variable z for liberating ε(p) and ∇u-p out from the nondifferentiable 1-norms, respectively. Then problem (16) can be transformed into the following equivalent constrained problem:
(18)argminu,pα0∥y∥1+α1∥z∥1+IΦ(w)s.t.kkklu=w,ε(p)=y,∇u-p=z.In order to liberateu out from the feasible set constraint, we introduce auxiliary variable w. Similar operation can also be found in [10]. Without this operation, we should resort to an inner iterative scheme to update the regularization parameter.The corresponding augmented Lagrangian functional of (18) is defined as
(19)LA(u,p,w,y,z;μ,ξ,η)≜IΦ(w)-〈μ,w-u〉+β12∥w-u∥22+α0∥y∥1-〈ξ,y-ε(p)〉+β22∥y-ε(p)∥22+α1∥z∥1-〈η,z-∇u+p〉+β32∥z-∇u+p∥22,
where μ, ξ, and η are Lagrange multipliers and β1, β2, and β3 are penalty parameters which should be positive. According to the classical ADMM, we should solve the following iterative scheme:
(20)uk+1=argminuLA(u,pk,wk,yk,zk;μk,ξk,ηk),pk+1=argminpLA(uk+1,p,wk,yk,zk;μk,ξk,ηk),wk+1=argminwLA(uk+1,pk+1,w,yk,zk;μk,ξk,ηk),yk+1=argminyLA(uk+1,pk+1,wk+1,y,zk;μk,ξk,ηk),zk+1=argminzLA(uk+1,pk+1,wk+1,yk+1,z;μk,ξk,ηk),μk+1=μk-β1(wk+1-uk+1),ξk+1=ξk-β2(yk+1-ε(uk+1)),ηk+1=ηk-β3(zk+1-∇uk+1+pk+1).


### 3.2. Solution of the Subproblems

With the auxiliaryw, the u subproblem becomes quadratic and irrelevant to the constraint of the feasible set. It allows the following objective:
(21)uk+1=argminuβ12∥wk-u-μkβ1∥22+β32∥zk-∇u+pk-ηkβ3∥22.
The minimization problem (21) can be solved by the following equation:
(22)uk+1=(β1β3I+∇T∇)-1(+∇T(zk+pk-ηkβ3)(zk+pk-ηkβ3)β1β3(wk-μkβ1)kkkkkkkkkkkkkkkk+∇T(zk+pk-ηkβ3)).
With the circulant boundary condition of images, we can solve (22) with several FFTs and IFFTs [8, 10].Following the same way, the subproblem with respect top is also quadratic and we have the objective functional as follows:
(23)pk+1=argminpβ22∥yk-ε(p)-ξkβ2∥22+β32∥zk-∇uk+1+p-ηkβ3∥22.
Then, for p1k+1, we have
(24)(β2∇x1T∇x1+β22∇x2T∇x2+β3I)p1k+1=β2[∇x1T(y1k-ξ1kβ2)+∇x2T(y3k-ξ3kβ2-12∇x1p2k)]+β3(∇x1uk+1+η1kβ3-z1k),
and for p2k+1, we have
(25)(β2∇x2T∇x2+β22∇x1T∇x1+β3I)p2k+1=β2[∇x2T(y2k-ξ2kβ2)+∇x1T(y3k-ξ3kβ2-12∇x2p1k+1)]+β3(∇x2uk+1+η2kβ3-z2k),
where p1 and p2 are the combinations of pi,j,1 and pi,j,2 (0≤i≤m, 0≤j≤n), respectively. Similar to the solution of (22), problems (24) and (25) can also be solved conveniently through several FFTs and IFFTs under the assumption of the circulant boundary condition.The subproblem fory can be written as
(26)yk+1=argminy∥y∥1+β22α0∥y-ε(pk+1)-ξkβ2∥22.
Problem (26) can be solved component-wisely through the following 4-dimensional shrinkage operation:
(27)yi,jk+1=max⁡⁡{∥ε(pk+1)i,j+ξi,jkβ2∥2-α0β2,0}×ε(pk+1)i,j+(ξi,jk/β2)∥ε(pk+1)i,j+(ξi,jk/β2)∥2lllkkill0≤i≤m,0≤j≤n.Thez subproblem is given by
(28)zk+1=argminz∥z∥1+β32α1∥z-∇uk+1+pk+1-ηkβ3∥22,
and it can be solved component-wisely through the following 2-dimensional shrinkage operation:
(29)zi,jk+1=max⁡⁡{∥(∇uk+1)i,j+ηi,jkβ3-pi,jk+1∥2-α1β3,0}×(∇uk+1)i,j+(ηi,jk/β3)-pi,jk+1∥(∇uk+1)i,j+(ηi,jk/β3)-pi,jk+1∥2,kkkkkkkkkkkk0≤i≤m,0≤j≤n.The subproblem with respect tow can be written as
(30)wk+1=argminwIΦ(w)+β12∥w-(uk+1+μkβ1)∥22=argminwλk+12∥w-u0∥22+β12∥w-(uk+1+μkβ1)∥22.
Consequently, the solution of problem (30) is
(31)wk+1=λk+1u0+β1(uk+1+(μk/β1))λk+1+β1.
The solution of λk+1 falls into two cases according to the range of uk+1+(μk/β1). On one hand, if
(32)∥uk+1+μkβ1-u0∥22≤c,
we can set λk+1=0, and, obviously, wk+1=uk+1+μk/β1 satisfies the feasible set constraint. On the other hand, if (32) is not true, wk+1 should fulfill the following equation:
(33)∥wk+1-u0∥22=c.
Substituting (31) into (33), we get
(34)λk+1=β1∥uk+1+(μk/β1)-u0∥2c-β1.
The resulting image denoising algorithm is summarized in Algorithm 1 TGV2ID-ADMM.Algorithm 1:
Algorithm TGV
2
ID-ADMM:  Second-order TGV-regularized image denoising with ADMM.
Input:  u0,c.
(1)Initialize u0,p0,w0,y0,z0,μ0,ξ0,η0. Set k=0 and β1,β2,β3>0.
(2)while  stopping criterion is not satisfied,do
(3) Compute  uk+1 according to (22);
(4) Compute  pk+1 according to (24) and (25);
(5) Compute  yk+1 according to (27);
(6) Compute  zk+1 according to (29);
(5)if(32) holds,then
(6) λk+1=0 and wk+1=uk+1+(μk/β1);
(7)else
(8) Update λk+1 and wk+1 according to (34) and (31);
(9)end  if
(10)   Updateμk+1,ξk+1, and  ηk+1 according to (20);
(11)k=k+1;
(12)end  while
(13)returnuk+1.The adoption of the variablew is essential for the adaptive estimate of the regularization parameter λ. With the assistance of w, u is liberated out from the constraint of the feasible set. Thus, the update of λ is free from the disturbance of the update of u, and a closed form for updating λ is achieved in each step without inner iteration. From functionals (16) and (18) we learn that, by setting (α0,α1)=(0,1), (y,p)=(0,0), and β2=0, Algorithm TGV2ID-ADMM will degenerate to a TV-based denoising algorithm, and we denote this case as TGV1ID-ADMM.The convergence of Algorithm TGV2ID-ADMM follows from the convergence analysis for the TV-based ADMM in [11, 25], due to the convex property of TGVα2. In this paper, we do not repeat the lengthy analysis procedure. However, we have the following essential convergence theorem for the proposed method.Theorem 1.
For fixedβ1,β2,β3>0, the sequence {uk,pk,wk,yk,zk,μk,ξk,ηk,λk} generated by Algorithm TGV2ID-ADMM from any initial point (u0,p0,w0,y0,z0,μ0,ξ0,η0) converges to (u*,p*,w*,y*,z*,μ*,ξ*,η*,λ*), where u* is the solution of functional (15) and λ* is the regularization parameter corresponding to the feasible set constraint u∈Φ.


## 3.1. The Augmented Lagrangian Model of Adaptive TGV-Based Denoising

Problem (15) can be transformed into an unconstrained problem, with the following discontinuous objective functional:
(16)argminu,pα0∥ε(p)∥1+α1∥∇u-p∥1+IΦ(u),
where IΦ(u) is the indicator function of the feasible set defined by
(17)IΦ(u)={0,ifu∈Φ≜{u:∥u-u0∥22≤c};+∞otherwise.
Note that Φ is a closed Euclidean ball centered at u0 with radius c.The solution of problem (16) suffers from its nonlinearity and nondifferentiability. Referring to the variable splitting, we introduce three auxiliary variables to simplify the solution process of (16): a variable w for liberating u out from the constraint of the feasible set; a variable y and a variable z for liberating ε(p) and ∇u-p out from the nondifferentiable 1-norms, respectively. Then problem (16) can be transformed into the following equivalent constrained problem:
(18)argminu,pα0∥y∥1+α1∥z∥1+IΦ(w)s.t.kkklu=w,ε(p)=y,∇u-p=z.In order to liberateu out from the feasible set constraint, we introduce auxiliary variable w. Similar operation can also be found in [10]. Without this operation, we should resort to an inner iterative scheme to update the regularization parameter.The corresponding augmented Lagrangian functional of (18) is defined as
(19)LA(u,p,w,y,z;μ,ξ,η)≜IΦ(w)-〈μ,w-u〉+β12∥w-u∥22+α0∥y∥1-〈ξ,y-ε(p)〉+β22∥y-ε(p)∥22+α1∥z∥1-〈η,z-∇u+p〉+β32∥z-∇u+p∥22,
where μ, ξ, and η are Lagrange multipliers and β1, β2, and β3 are penalty parameters which should be positive. According to the classical ADMM, we should solve the following iterative scheme:
(20)uk+1=argminuLA(u,pk,wk,yk,zk;μk,ξk,ηk),pk+1=argminpLA(uk+1,p,wk,yk,zk;μk,ξk,ηk),wk+1=argminwLA(uk+1,pk+1,w,yk,zk;μk,ξk,ηk),yk+1=argminyLA(uk+1,pk+1,wk+1,y,zk;μk,ξk,ηk),zk+1=argminzLA(uk+1,pk+1,wk+1,yk+1,z;μk,ξk,ηk),μk+1=μk-β1(wk+1-uk+1),ξk+1=ξk-β2(yk+1-ε(uk+1)),ηk+1=ηk-β3(zk+1-∇uk+1+pk+1).


## 3.2. Solution of the Subproblems

With the auxiliaryw, the u subproblem becomes quadratic and irrelevant to the constraint of the feasible set. It allows the following objective:
(21)uk+1=argminuβ12∥wk-u-μkβ1∥22+β32∥zk-∇u+pk-ηkβ3∥22.
The minimization problem (21) can be solved by the following equation:
(22)uk+1=(β1β3I+∇T∇)-1(+∇T(zk+pk-ηkβ3)(zk+pk-ηkβ3)β1β3(wk-μkβ1)kkkkkkkkkkkkkkkk+∇T(zk+pk-ηkβ3)).
With the circulant boundary condition of images, we can solve (22) with several FFTs and IFFTs [8, 10].Following the same way, the subproblem with respect top is also quadratic and we have the objective functional as follows:
(23)pk+1=argminpβ22∥yk-ε(p)-ξkβ2∥22+β32∥zk-∇uk+1+p-ηkβ3∥22.
Then, for p1k+1, we have
(24)(β2∇x1T∇x1+β22∇x2T∇x2+β3I)p1k+1=β2[∇x1T(y1k-ξ1kβ2)+∇x2T(y3k-ξ3kβ2-12∇x1p2k)]+β3(∇x1uk+1+η1kβ3-z1k),
and for p2k+1, we have
(25)(β2∇x2T∇x2+β22∇x1T∇x1+β3I)p2k+1=β2[∇x2T(y2k-ξ2kβ2)+∇x1T(y3k-ξ3kβ2-12∇x2p1k+1)]+β3(∇x2uk+1+η2kβ3-z2k),
where p1 and p2 are the combinations of pi,j,1 and pi,j,2 (0≤i≤m, 0≤j≤n), respectively. Similar to the solution of (22), problems (24) and (25) can also be solved conveniently through several FFTs and IFFTs under the assumption of the circulant boundary condition.The subproblem fory can be written as
(26)yk+1=argminy∥y∥1+β22α0∥y-ε(pk+1)-ξkβ2∥22.
Problem (26) can be solved component-wisely through the following 4-dimensional shrinkage operation:
(27)yi,jk+1=max⁡⁡{∥ε(pk+1)i,j+ξi,jkβ2∥2-α0β2,0}×ε(pk+1)i,j+(ξi,jk/β2)∥ε(pk+1)i,j+(ξi,jk/β2)∥2lllkkill0≤i≤m,0≤j≤n.Thez subproblem is given by
(28)zk+1=argminz∥z∥1+β32α1∥z-∇uk+1+pk+1-ηkβ3∥22,
and it can be solved component-wisely through the following 2-dimensional shrinkage operation:
(29)zi,jk+1=max⁡⁡{∥(∇uk+1)i,j+ηi,jkβ3-pi,jk+1∥2-α1β3,0}×(∇uk+1)i,j+(ηi,jk/β3)-pi,jk+1∥(∇uk+1)i,j+(ηi,jk/β3)-pi,jk+1∥2,kkkkkkkkkkkk0≤i≤m,0≤j≤n.The subproblem with respect tow can be written as
(30)wk+1=argminwIΦ(w)+β12∥w-(uk+1+μkβ1)∥22=argminwλk+12∥w-u0∥22+β12∥w-(uk+1+μkβ1)∥22.
Consequently, the solution of problem (30) is
(31)wk+1=λk+1u0+β1(uk+1+(μk/β1))λk+1+β1.
The solution of λk+1 falls into two cases according to the range of uk+1+(μk/β1). On one hand, if
(32)∥uk+1+μkβ1-u0∥22≤c,
we can set λk+1=0, and, obviously, wk+1=uk+1+μk/β1 satisfies the feasible set constraint. On the other hand, if (32) is not true, wk+1 should fulfill the following equation:
(33)∥wk+1-u0∥22=c.
Substituting (31) into (33), we get
(34)λk+1=β1∥uk+1+(μk/β1)-u0∥2c-β1.
The resulting image denoising algorithm is summarized in Algorithm 1 TGV2ID-ADMM.Algorithm 1:
Algorithm TGV
2
ID-ADMM:  Second-order TGV-regularized image denoising with ADMM.
Input:  u0,c.
(1)Initialize u0,p0,w0,y0,z0,μ0,ξ0,η0. Set k=0 and β1,β2,β3>0.
(2)while  stopping criterion is not satisfied,do
(3) Compute  uk+1 according to (22);
(4) Compute  pk+1 according to (24) and (25);
(5) Compute  yk+1 according to (27);
(6) Compute  zk+1 according to (29);
(5)if(32) holds,then
(6) λk+1=0 and wk+1=uk+1+(μk/β1);
(7)else
(8) Update λk+1 and wk+1 according to (34) and (31);
(9)end  if
(10)   Updateμk+1,ξk+1, and  ηk+1 according to (20);
(11)k=k+1;
(12)end  while
(13)returnuk+1.The adoption of the variablew is essential for the adaptive estimate of the regularization parameter λ. With the assistance of w, u is liberated out from the constraint of the feasible set. Thus, the update of λ is free from the disturbance of the update of u, and a closed form for updating λ is achieved in each step without inner iteration. From functionals (16) and (18) we learn that, by setting (α0,α1)=(0,1), (y,p)=(0,0), and β2=0, Algorithm TGV2ID-ADMM will degenerate to a TV-based denoising algorithm, and we denote this case as TGV1ID-ADMM.The convergence of Algorithm TGV2ID-ADMM follows from the convergence analysis for the TV-based ADMM in [11, 25], due to the convex property of TGVα2. In this paper, we do not repeat the lengthy analysis procedure. However, we have the following essential convergence theorem for the proposed method.Theorem 1.
For fixedβ1,β2,β3>0, the sequence {uk,pk,wk,yk,zk,μk,ξk,ηk,λk} generated by Algorithm TGV2ID-ADMM from any initial point (u0,p0,w0,y0,z0,μ0,ξ0,η0) converges to (u*,p*,w*,y*,z*,μ*,ξ*,η*,λ*), where u* is the solution of functional (15) and λ* is the regularization parameter corresponding to the feasible set constraint u∈Φ.


## 4. Experiment Results

In this section, we illustrate the effectiveness of the proposed algorithm on suppressing staircasing effect and removing Gaussian noise in image. Besides, we also show the robustness of the results with respect to the penalty parameters. We performed our algorithm under MATLAB v7.8.0 and Windows 7 on a PC with an Intel Core (TM) i5 CPU at 3.20 GHz and 8 GB of RAM.The root mean squared error (RMSE) and the peak signal-to-noise ratio (PSNR) used in comparison are defined as(35)RMSE=∥u-uclean∥2mn,PSNR=10log10(2552·mn∥u-uclean∥22),
where uclean is the original image that contains no noise. Besides, in subsections 4.1 and 4.2, we set the penalty parameters as β1=10(BSNR/10-1)×β and β2=β3=β=0.3 for TGV2ID-ADMM (β2=0 for TGV1ID-ADMM) to achieve consistently promising result with fast speed, where BSNR is the blurred signal-to-noise ratio defined by BSNR = 10 log10 (var(u0)/σ2) (var(u0) denotes the variance of u0).


### 4.1. Staircasing Effect Reduction by the Proposed Method

We first compare Algorithm TGV2ID-ADMM with Algorithm TGV1ID-ADMM to illustrate the effectiveness of TGVα2 model in staircasing effect reduction. We use ∥uk+1-uk∥22/∥uk∥22≤10-6 as the stopping criteria for these two algorithms, where uk denotes the restored result in the kth iteration. For the second-order case, we set (α0,α1)=(3,1), whereas for the one-order case, we set (α0,α1)=(0,1).In this experiment, we use a synthetic piecewise affine image shown in Figure1 as the test image. The original image is contaminated by Gaussian noise of standard variance σ=15 at first. Then we imply TGV2ID-ADMM and TGV1ID-ADMM to remove the noise. Table 1 shows the results in terms of RMSE, PSNR, total iterations, and CPU time. The ground truth, noised, and restored images by the two algorithms are displayed in Figure 1. Furthermore, for better visualization, we additionally provide the three-dimensional close-ups of the marked regions of the two restored images in Figure 1. From Table 1 we observe that TGV2ID-ADMM does better than TGV1ID-ADMM in terms of both RMSE and PSNR. Figure 1 shows that the denoised image of TGV2ID-ADMM almost contains no artificial edges in affine regions. In contrast, the restored result of TGV1ID-ADMM contains obvious staircasing effect in affine regions. The three-dimensional closed-ups vividly demonstrate this phenomenon. This illustrates that our TGV-based algorithm is effective in staircasing effect reduction.Table 1
Results of the staircasing effect reduction experiment.
Piecewise affine image denoising withσ=15
Model
RMSE
PSNR
Iterations
Time (s)
Noised
14.99
24.61
—
—
TGV2
1.89
42.61
155
18.21
TGV1
2.40
39.81
128
7.76First row: ground truth and Gaussian noised (σ=15) piecewise affine images; second row: restored images by TGV2ID-ADMM and TGV1ID-ADMM; third row: three-dimensional close-ups of the marked regions in the two restored images.
(a)
(b)
(c)Table1 also shows that, to accomplish the denoising task, TGV2ID-ADMM usually costs more CPU time than TGV1ID-ADMM, since TGVα2 model involves much more calculation. However, the cost is worthy due to the impressive improvement on both quantitative and qualitative restoration quality. Figure 2 displays the evolutions of λs and PSNRs achieved by the two algorithms. It is learnt that, the regularization parameters of both converge to the optimal points at last, which guarantees the automatic implementation of the two algorithms.Evolutions ofλs (a) and PSNRs (b) achieved by TGV2ID-ADMM (TGV2) and TGV1ID-ADMM (TGV1).
(a)
(b)


### 4.2. Comparison in Accuracy

In this subsection, we compare TGV2ID-ADMM with the other two famous adaptive TV-based denoising algorithms: Chambolle’s projection algorithm [23] and Split Bregman algorithm [24], both possessing public online implementations at “http://www.ipol.im/”. Two natural images, Lena and Peppers both of size 512 × 512 shown in Figure 3, are used for comparison. The parameter setting for TGV2ID-ADMM is the same as that in the previous subsection. We obtain the test results of the two competitors through online experimental operation.Test images: Lena and Peppers.
(a)
(b)We add Gaussian noise of standard variances of 20, 30, and 40 to Lena and Peppers to obtain the noised observations, respectively. Then we apply these three algorithms to restore the noisy images. Table2 shows the comparison results in terms of RMSE and PSNR. The best result for each comparison item is highlighted in bold type font. Table 2 shows that TGV2ID-ADMM holds superiority on both RMSE and PSNR for all the tested cases. Figure 4 displays the noised Lena under Gaussian noise of σ=30 and the restorations by the three algorithms, whereas Figure 5 exhibits the noised Peppers under Gaussian noise of σ=40 and the corresponding restorations. Figures 4 and 5 demonstrate that TGV2ID-ADMM obtains results with better visual impression and efficiently suppresses the staircasing effect. In contrast, both TV-based Chambolle’s projection algorithm and TV-based Split Bregman algorithm achieve results with obvious staircasing effect. Since we apply test images with different levels of noise, the robustness of our algorithm towards the noise level is verified to a certain extent.Table 2
Comparison results in accuracy.
σ
Noised
TGV2ID-ADMM
Chambolle
Split Bregman
RMSE
PSNR
RMSE
PSNR
RMSE
PSNR
RMSE
PSNR
Lena512×512
20
19.99
22.11
7.09
31.12
7.40
30.75
8.05
30.01
30
30.00
18.59
8.67
29.37
9.45
28.62
9.95
28.17
40
39.97
16.10
10.06
28.07
10.54
27.67
11.66
26.80
Peppers512×512
20
19.96
22.13
7.05
31.17
7.38
30.77
7.82
30.27
30
29.98
18.59
8.68
29.36
9.43
28.64
9.92
28.20
40
40.10
16.07
9.74
28.36
10.74
27.51
12.38
26.28Noised Lena under Gaussian noise ofσ=30 and the restorations by TGV2ID-ADMM, Chambolle, and Split Bregman, respectively.
(a)
(b)
(c)
(d)Noised Peppers under Gaussian noise ofσ=40 and the restorations by TGV2ID-ADMM, Chambolle, and Split Bregman, respectively.
(a)
(b)
(c)
(d)


### 4.3. Solution Robustness with Respect to the Penalty Parameters

Although the positive assumption of penalty parameters is sufficient for the convergence of ADMM, the results of ADMM are commonly influenced by the choice of the penalty parameters to a certain extent in practice. As suggested by a referee, we add an experiment to show the robustness of the results of TGV2ID-ADMM with respect to the penalty parameters, under the two denoising background problems mentioned above, that is, the Lena denoising problem under Gaussian noise of σ=30 and the Peppers denoising problem under Gaussian noise of σ=40. We still set β1=10(BSNR/10-1)×β and β2=β3=β but change β from 0.01 to 1 with a step size of 0.01. In Figure 6, we plot PSNR versus β for the denoised Lena and Peppers. Figure 6 demonstrates that the optimal β should be focalized in [0.2,0.3] and its location is robust towards the variation of image and noise level. The results of our method possess sufficient robustness with respect to the variation of penalty parameters to a certain extent, since the absolute error between the maximum and the minimum of PSNR is less than 0.18 dB in the experiment, and this error could not introduce obvious distinction in visual quality. In the former two experiments, the setting of β=0.3 is approximately optimal for the proposed algorithm.PSNR versusβ for (a) Lena image under Gaussian noise of σ=30 and (b) Peppers image under Gaussian noise of σ=40, respectively. The absolute error between the maximum and the minimum of PSNR for each image is less than 0.18 dB.
(a)
Lena image under Gaussian noise ofσ=30
(b)
Peppers image under Gaussian noise ofσ=40


## 4.1. Staircasing Effect Reduction by the Proposed Method

We first compare Algorithm TGV2ID-ADMM with Algorithm TGV1ID-ADMM to illustrate the effectiveness of TGVα2 model in staircasing effect reduction. We use ∥uk+1-uk∥22/∥uk∥22≤10-6 as the stopping criteria for these two algorithms, where uk denotes the restored result in the kth iteration. For the second-order case, we set (α0,α1)=(3,1), whereas for the one-order case, we set (α0,α1)=(0,1).In this experiment, we use a synthetic piecewise affine image shown in Figure1 as the test image. The original image is contaminated by Gaussian noise of standard variance σ=15 at first. Then we imply TGV2ID-ADMM and TGV1ID-ADMM to remove the noise. Table 1 shows the results in terms of RMSE, PSNR, total iterations, and CPU time. The ground truth, noised, and restored images by the two algorithms are displayed in Figure 1. Furthermore, for better visualization, we additionally provide the three-dimensional close-ups of the marked regions of the two restored images in Figure 1. From Table 1 we observe that TGV2ID-ADMM does better than TGV1ID-ADMM in terms of both RMSE and PSNR. Figure 1 shows that the denoised image of TGV2ID-ADMM almost contains no artificial edges in affine regions. In contrast, the restored result of TGV1ID-ADMM contains obvious staircasing effect in affine regions. The three-dimensional closed-ups vividly demonstrate this phenomenon. This illustrates that our TGV-based algorithm is effective in staircasing effect reduction.Table 1
Results of the staircasing effect reduction experiment.
Piecewise affine image denoising withσ=15
Model
RMSE
PSNR
Iterations
Time (s)
Noised
14.99
24.61
—
—
TGV2
1.89
42.61
155
18.21
TGV1
2.40
39.81
128
7.76First row: ground truth and Gaussian noised (σ=15) piecewise affine images; second row: restored images by TGV2ID-ADMM and TGV1ID-ADMM; third row: three-dimensional close-ups of the marked regions in the two restored images.
(a)
(b)
(c)Table1 also shows that, to accomplish the denoising task, TGV2ID-ADMM usually costs more CPU time than TGV1ID-ADMM, since TGVα2 model involves much more calculation. However, the cost is worthy due to the impressive improvement on both quantitative and qualitative restoration quality. Figure 2 displays the evolutions of λs and PSNRs achieved by the two algorithms. It is learnt that, the regularization parameters of both converge to the optimal points at last, which guarantees the automatic implementation of the two algorithms.Evolutions ofλs (a) and PSNRs (b) achieved by TGV2ID-ADMM (TGV2) and TGV1ID-ADMM (TGV1).
(a)
(b)


## 4.2. Comparison in Accuracy

In this subsection, we compare TGV2ID-ADMM with the other two famous adaptive TV-based denoising algorithms: Chambolle’s projection algorithm [23] and Split Bregman algorithm [24], both possessing public online implementations at “http://www.ipol.im/”. Two natural images, Lena and Peppers both of size 512 × 512 shown in Figure 3, are used for comparison. The parameter setting for TGV2ID-ADMM is the same as that in the previous subsection. We obtain the test results of the two competitors through online experimental operation.Test images: Lena and Peppers.
(a)
(b)We add Gaussian noise of standard variances of 20, 30, and 40 to Lena and Peppers to obtain the noised observations, respectively. Then we apply these three algorithms to restore the noisy images. Table2 shows the comparison results in terms of RMSE and PSNR. The best result for each comparison item is highlighted in bold type font. Table 2 shows that TGV2ID-ADMM holds superiority on both RMSE and PSNR for all the tested cases. Figure 4 displays the noised Lena under Gaussian noise of σ=30 and the restorations by the three algorithms, whereas Figure 5 exhibits the noised Peppers under Gaussian noise of σ=40 and the corresponding restorations. Figures 4 and 5 demonstrate that TGV2ID-ADMM obtains results with better visual impression and efficiently suppresses the staircasing effect. In contrast, both TV-based Chambolle’s projection algorithm and TV-based Split Bregman algorithm achieve results with obvious staircasing effect. Since we apply test images with different levels of noise, the robustness of our algorithm towards the noise level is verified to a certain extent.Table 2
Comparison results in accuracy.
σ
Noised
TGV2ID-ADMM
Chambolle
Split Bregman
RMSE
PSNR
RMSE
PSNR
RMSE
PSNR
RMSE
PSNR
Lena512×512
20
19.99
22.11
7.09
31.12
7.40
30.75
8.05
30.01
30
30.00
18.59
8.67
29.37
9.45
28.62
9.95
28.17
40
39.97
16.10
10.06
28.07
10.54
27.67
11.66
26.80
Peppers512×512
20
19.96
22.13
7.05
31.17
7.38
30.77
7.82
30.27
30
29.98
18.59
8.68
29.36
9.43
28.64
9.92
28.20
40
40.10
16.07
9.74
28.36
10.74
27.51
12.38
26.28Noised Lena under Gaussian noise ofσ=30 and the restorations by TGV2ID-ADMM, Chambolle, and Split Bregman, respectively.
(a)
(b)
(c)
(d)Noised Peppers under Gaussian noise ofσ=40 and the restorations by TGV2ID-ADMM, Chambolle, and Split Bregman, respectively.
(a)
(b)
(c)
(d)


## 4.3. Solution Robustness with Respect to the Penalty Parameters

Although the positive assumption of penalty parameters is sufficient for the convergence of ADMM, the results of ADMM are commonly influenced by the choice of the penalty parameters to a certain extent in practice. As suggested by a referee, we add an experiment to show the robustness of the results of TGV2ID-ADMM with respect to the penalty parameters, under the two denoising background problems mentioned above, that is, the Lena denoising problem under Gaussian noise of σ=30 and the Peppers denoising problem under Gaussian noise of σ=40. We still set β1=10(BSNR/10-1)×β and β2=β3=β but change β from 0.01 to 1 with a step size of 0.01. In Figure 6, we plot PSNR versus β for the denoised Lena and Peppers. Figure 6 demonstrates that the optimal β should be focalized in [0.2,0.3] and its location is robust towards the variation of image and noise level. The results of our method possess sufficient robustness with respect to the variation of penalty parameters to a certain extent, since the absolute error between the maximum and the minimum of PSNR is less than 0.18 dB in the experiment, and this error could not introduce obvious distinction in visual quality. In the former two experiments, the setting of β=0.3 is approximately optimal for the proposed algorithm.PSNR versusβ for (a) Lena image under Gaussian noise of σ=30 and (b) Peppers image under Gaussian noise of σ=40, respectively. The absolute error between the maximum and the minimum of PSNR for each image is less than 0.18 dB.
(a)
Lena image under Gaussian noise ofσ=30
(b)
Peppers image under Gaussian noise ofσ=40


## 5. Concluding Remarks

We propose an adaptive TGV-based model for noise removal in this paper. The variable splitting (VS) and the classical augmented Lagrangian method are used to handle the proposed model. From the experimental results, we observe that the proposed algorithm is effective in suppressing staircasing effect and preserving edges in images, and it is superior to some other famed adaptive denoising methods both in quantitative and in qualitative assessment. Besides, our work can be smoothly generalized to image deblurring problems.



---

*Source: 157893-2014-07-10.xml*