# Image Restoration Based on the Hybrid Total-Variation-Type Model

Authors: Baoli Shi; Zhi-Feng Pang; Yu-Fei Yang
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

We propose a hybrid total-variation-type model for the image restoration problem based on combining advantages of the ROF model with the LLT model. Since twoL1-norm terms in the proposed model make it difficultly solved by using some classically numerical methods directly, we first employ the alternating direction method of multipliers (ADMM) to solve a general form of the proposed model. Then, based on the ADMM and the Moreau-Yosida decomposition theory, a more efficient method called the proximal point method (PPM) is proposed and the convergence of the proposed method is proved. Some numerical results demonstrate the viability and efficiency of the proposed model and methods.


---

## Body

## 1. Introduction

Image restoration is of momentous significance in coherent imaging systems and various image processing applications. The goal is to recover the real image from the deteriorated image, for example, image denoising, image deblurring, image inpainting, and so forth; see [1–3] for details.For the additive noisy image, many denoising models have been proposed based on PDEs or variational methods over the last decades [1–3]. The essential idea for this class of models is to filter out the noise in an image while preserving significant features such as edges and textures. However, due to the ill-posedness of the restoration problem, we have to employ some regularization methods [4] to overcome it. The general form of regularization methods consists in minimizing an energy functional of the following form:
(1.1)F(u-)=λ2∥u--u-0∥𝒳2+ℛ(u-)
in Banach space 𝒳, where λ is the regularization parameter, ℛ, is a regularization term, u-0 is the observed image, and u- is the image to be restored. The development of the energy functional (1.1) actually profits from the ROF model [5] which is of the following form:
(1.2)minu-λ12∥u--u-0∥L22+|u-|BV,
where |u-|BV:=∫Ω‍|∇u-|dx. In the case without confusion, for simplification we omit the open set Ω with Lipschitz boundary for ∥·∥L2(Ω) and |·|L1(Ω). Due to edge-preserving property of the term |u-|BV, this model has been extended to other sorts of image processing problems such as
(1.3)minu-λ22∥Ku--u-0∥L22+|u-|BV
for image deblurring [6] and
(1.4)minu-λ32∥u--u-0∥L2(Ω∖D)2+|u-|BV
for image inpainting [2, 7], where K is a blurring operator and D is the inpainting domain. Furthermore, this model was applied to restore the multiplicative noisy image which usually appears in various image processing applications such as in laser images, ultrasound images [8], synthetic aperture radar (SAR) [9], and medical ultrasonic images [10]. One of the models based on BV was proposed by Huang et al. (HNW) [11] with the form
(1.5)minzα∫Ω‍(z+fe-z)dx+|z|BV,
where z is the logarithmic transformation of u and |z|BV keeps the total variation property related to |ez|BV and α>0. Here u satisfies f=uη for the noise η.However, as we all know the|·|BV term usually reduces the computational solution of the above models to be piecewise constant, which is also called the staircasing effect in smooth regions of the image. The staircase effect implies to produce new edges that do not exist in the true image so that the restored image is unsatisfactory to the eye. To overcome this drawback, some high-order models [12–15] have been proposed such as a model proposed by Lysaker, Lundervold, and Tai (the LLT model) with the following form:
(1.6)minu-β12∥u--u-0∥L22+|u-|BV2,
where |u-|BV2:=∫Ω‍|∇2u-|dx for the Hessian operator ∇2. However, these classes of models can blur the edges of the image in the course of restoration. Therefore, it is a natural choice to combine advantages of the ROF model and the LLT model if we want to preserve edges while avoiding the staircase effect in smooth regions. One of convex combinations between the BV and BV2 was proposed by Lysaker et al. [13] to restore the image with additive noise. But their model is not quite intuitive due to lack of gradient information in the weighting function. Since the edge detector function g:=g(u0):=1/(1+ϱ|∇(u0*Gσ(x))|) with Gσ(x)=(1/2πσ2)exp(-(x2/2σ2)) can depict the information of edges, we can employ it as a balance function; that is, we can apply the following model:
(1.7)minuβ2∥u-u-0∥L22+|(1-g)u|BV+|gu|BV2
to restore the noisy image [16]. Obviously, |(1-g)u|BV tends to be predominant where edges most likely appear and |gu|BV2 tends to be predominant at the locations with smooth signals. Based on the advantages of the hybrid model (1.7), we also extend it to the image restoration models (1.3)–(1.5) in this paper.Another topic for image restoration is to find some efficient methods to solve the above proposed models. In fact, there are many different methods based on PDE or convex optimization to solve the minimization problem (1.1) by means of the specific models (1.2)–(1.6). For example, for the purpose of solving the ROF model (1.2) or the LLT model (1.6), we can use the gradient descent method [5, 13], the Chambolle’s dual method [13, 17, 18], the primal and dual method [19–21], the second order cone programming method [22], the multigrid method [23], operator splitting method [24–26], the inverse scale method [27], and so forth. However, different from the ROF model (1.2) and the LLT model (1.6), the model (1.7) includes two L1-norm terms which make it solved more difficultly. More generally, the model (1.7) can be fell into the following framework:
(1.8)minuθ2∥u-f∥L22+h1(Λ1u)+h2(Λ2u),
where h1,h2:𝒳→R are proper, convex, and lower semicontinuous (l.s.c.), Λ1 and Λ2 are bounded linear operators, and θ is a parameter. A specific form of (1.8) was considered by Afonso and Bioucas-Dias in [28] where a Bregman iterative method was proposed to solve the model with the combination of the L1-norm term and the total variation (TV) term. Actually this splitting Bregman method is formally equivalent to the alternating direction method of multipliers (ADMM) [24, 29–34]. However, the ADMM ineluctably tends to solve some subproblems which correspond to the related modified problems. Furthermore, these make us obtain the numerical results by requiring much more computational cost. In order to obtain an efficient numerical method, it is a natural choice to avoid solving the related subproblems. In this paper, we propose a proximal point method (PPM) which can be deduced from the ADMM. This deduction is based on the connection that the sum of the projection operator and the shrinkage operator is equal to the identity operator; it is known as the Moreau-Yosida decomposition Theorem  31.5 in [35]. Then the PPM not only keeps the advantages of the ADMM but also requires much less computational cost. This implies that the PPM is much more fast and efficient, especially for the larger scale images. Furthermore, using the monotone operator theory, we give the convergence analysis of the proposed method. Moreover, we extend the PPM to solve the model (1.7) to image deblurring, image inpainting, and the multiplicative noisy image restoration. Experimental results show that the restored images generated by the proposed models and methods are desirable.The paper is organized as follows. In Section2, we recall some knowledge related to convex analysis. In Section 3, we first propose the ADMM to solve the problem (1.8) and then give the PPM to improve this method. In Section 4, we give some applications by using the proposed algorithms and also compare the related models and the proposed methods. Some concluding remarks are given in Section 5.


## 2. Notations and Definitions

Let us describe some notations and definitions used in this paper. For simplifying, we useX=Rn. Usually, we can set n=2 for the gray-scale images. The related contents can be referred to [1, 27, 36–43].Definition 2.1.
The operatorA:X→R is monotone if it satisfies
(2.1)(y1-y2,x1-x2)≥0
for all y1∈A(x1) and y2∈A(x2), and A is maximal if it is not strictly contained in any other monotone operator on X.Definition 2.2.
LetG:X→R∪{+∞} be a convex function. The subdifferential of G at a point y∈X is defined by
(2.2)∂G(y)={η:G(y)-G(x)≤(η,y-x),∀x∈X},
where η∈∂G(y) is called a subgradient. It reduces to the classical gradient if G(x) is differentiable.Definition 2.3.
Assume thatA is a maximal monotone operator. Denote by Hc its Yosida approximation:
(2.3)HcA=1c(I-JcA),
where I denotes the identity operator [38, 39] and JcA is the resolvent of A with the form as
(2.4)JcA=(I+cA)-1.Definition 2.4.
LetG:X→R be a convex proper function. The conjugate function of G is defined by
(2.5)G*(x*)=supx∈X{(x,x*)-G(x)}
for all x*∈X.Definition 2.5.
LetG:X→R∪{∞} be convex and t>0. The proximal mapping to G:X→X is defined by
(2.6)ProxtG(x):=argminy{G(y)+12t∥y-x∥L22}
for y∈X.Definition 2.6.
The projection operator𝒫ℬτ(·):X→X onto the closed disc ℬτ:={x∈X:|x|L1≤1/τ} is defined by
(2.7)𝒫ℬτ(x)=x|x|L1min(|x|L1,1τ),
where x∈X and τ>0.Definition 2.7.
The shrinkage operatorSτ(·):X→X is defined by
(2.8)𝒮τ(x)=x|x|L1max(|x|L1-1τ,0),
where we use the convention 0/0=0.Remark 2.8.
It is obvious that the functionG:X→R and its conjugate function G* satisfy the following relationship:
(2.9)ProxtG(x)+Prox1/tG*(x)=x,
for t>0. Especially, the projection operator 𝒫ℬτ(·) and the shrinkage operator Sτ(·) satisfy
(2.10)𝒫ℬτ(x)+Sτ(x)=x
for any x∈X. In fact, this corresponds to the classic Moreau-Yosida decomposition Theorem  31.5 in [35].


## 3. The Alternating Direction Method of Multipliers (ADMM) and the Proximal Point Method (PPM)

Variable splitting methods such as the ADMM [29–32, 44] and the operator splitting methods [42, 45, 46] have been recently used in the image, signal, and data processing community. The key of this class of methods is to transform the original problem into some subproblems so that we can easily solve these subproblems by employing some numerical methods. In this section we first consider to use the ADMM to solve the general minimization problem (1.8). However, the computational cost of the ADMM is tediously increased due to its looser form. In order to overcome this drawback, we thus change this method to a compacter form called the proximal method based on the relationship (2.9) in Remark 2.8.


### 3.1. The Alternating Direction Method of Multipliers (ADMM)

We now consider the following constrained problem:(3.1)minu,v,wθ2∥u-f∥L22+h1(v)+h2(z),s.t.v=Λ1u,z=Λ2u,
which is clearly equivalent to the constrained problem (1.8) in the feasible set {(u,v,z):v=Λ1u,andz=Λ2u}. Throughout the following subsections, we always assume that Λ1 and Λ2 are a surjective map. It seems that the problem (3.1) including three variables looks more complex than the original unconstrained problem (1.8). In fact, this problem can be solved more easily under the condition that h1 and h2 are nondifferentiable. In the augmented Lagrangian framework [33, 34, 47, 48], the problem (3.1) equivalently solves the following minimization Lagrangian function:
(3.2)minu,v,z,ζ1,ζ2,μ1,μ2L(u,v,z,ζ1,ζ2,μ1,μ2)=θ2∥u-f∥L22+h1(v)+(ζ1,Λ1u-v)+μ12∥Λ1u-v∥L22+h2(z)+(ζ2,Λ2u-z)+μ22∥Λ2u-z∥L22,
where ζi is the Lagrangian multiplier and μi is the penalty parameter for i=1, 2. Then we can use the following augmented Lagrangian method (ALM): (3.3a)(un,vn,zn)=argminu,v,zL(u,v,z,ζ1n-1,ζ2n-1,μ1,μ2),(3.3b)ζ1n=ζ1n-1+μ1(Λ1u-v),(3.3c)ζ2n=ζ2n-1+μ2(Λ2u-z)with choosing the original values ζ10 and ζ20 to solve (3.2). If we set din=ζin/μi for i=1, 2 and omit the terms which are independent of (un,vn,zn) in (3.3a), the above strategy (3.3a)–(3.3c) can be written as (3.4a)(un,vn,zn)=argminu,v,zθ2∥u-f∥L22+h1(v)+μ12∥Λ1u-v+d1n-1∥L22+h2(z),+μ22∥Λ2u-z+d2n-1∥L22,(3.4b)d1n=d1n-1+Λ1un-vn,(3.4c)d2n=d2n-1+Λ2un-znfor the original values d10 and d20. Then we can use the following ADMM to solve (3.4a)–(3.4c).Algorithm 3.1 (ADMM for solving (3.4a)–(3.4c)).
(1) Choose the original values:v0, z0, d10, and d20. Set θ, μ1, μ2>0 and n=1.
(2) Compute(un,vn,zn,d1n,d2n) by(3.5a)
u
n
=
arg
min
u
θ
2
∥
u
-
f
∥
L
2
2
+
μ
1
2
∥
Λ
1
u
-
v
n
-
1
+
d
1
n
-
1
∥
L
2
2
+
μ
2
2
∥
Λ
2
u
-
z
n
-
1
+
d
2
n
-
1
∥
L
2
2
︸
=
h
(
u
)
,
(3.5b)
v
n
=
arg
min
v
h
1
(
v
)
+
μ
1
2
∥
Λ
1
u
n
-
v
+
d
1
n
-
1
∥
L
2
2
,
(3.5c)
z
n
=
arg
min
z
h
2
(
z
)
+
μ
2
2
∥
Λ
2
u
n
-
z
+
d
2
n
-
1
∥
L
2
2
,
(3.5d)
d
1
n
=
d
1
n
-
1
+
Λ
1
u
n
-
v
n
,
(3.5e)
d
2
n
=
d
2
n
-
1
+
Λ
2
u
n
-
z
n
.(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Sinceh(u) is differentiable and strictly convex, we can get the unique solution of (3.5a) which satisfies (3.6a)θ(un-f)+μ1Λ1*(Λ1un-vn-1+d1n-1)+μ2Λ2*(Λ2un-zn-1+d2n-1)=0,(3.6b)0=w1n+μ1(vn-d1n-1-Λ1un),(3.6c)0=w2n+μ2(zn-d2n-1-Λ2un),(3.6d)d1n=d1n-1+Λ1un-vn,(3.6e)d2n=d2n-1+Λ2un-zn,where Λ1* and Λ2* are the adjoint operators of Λ1 and Λ2, w1n∈∂h1(vn) and w2n∈∂h2(zn), respectively. It follows that the solution un in (3.6a) can be directly obtained by the following explicit formulation:
(3.7)un=(θI+μ1Λ1*Λ1+μ2Λ2*Λ2)︸M-1(θf+μ1Λ1*(vn-1-d1n-1)+μ2Λ2*(zn-1-d2n-1)).
However, when the operator M is ill-posed, the solution is unsuitable or unavailable. Hence we have to go back to (3.6a) and to employ some iteration strategy such as the Gauss-Seidel method to solve this equation. On the other hand, it is obvious that (3.5b) and (3.5c) can be looked at as the proximal mapping, so the solutions of the minimization problems (3.5b) and (3.5c) can be obviously written as
(3.8)vn=Prox1/μ1h1(Λ1un+d1n-1),zn=Prox1/μ2h2(Λ2un+d2n-1).Theorem 3.2.
Assume that(u*,v*,z*,w1*,w2*) is the saddle point of the Lagrange function
(3.9)ℒ(u,v,z,w1,w2)=θ2∥u-f∥L22+h1(v)+(w1,Λ1u-v)+h2(z)+(w2,Λ2u-z).
Then u* is the solution of the minimization problem (1.8). Furthermore, the sequence {(un,vn,zn)} generated by Algorithm 3.1 converges to (u*,v*,z*).Notice that Algorithm3.1 can be actually looked at as the split Bregman method [25]. The based idea of this method is to introduce some intermediate variables so as to transform the original problem into some subproblems which are easily solved. The connection between the split Bregman method and the ADMM has been shown in [29, 49]. However, our algorithm considers the sum of three convex functions, which is more general than the related algorithms in [25, 49]. Furthermore, it must be noted that v and z are completely separated in (3.4a), so the two subproblems (3.5b) and (3.5c) are parallel. Therefore the convergence results of the ADMM can be applied here.


### 3.2. The Proximal Point Method

Though the ADMM in Algorithm3.1 can effectively solve the original problem (3.1), we have to solve five subproblems. This actually makes this method suffer from a looser form as in [25, 45] so that it can badly affect its numerical computation efficiency. In this subsection, we propose a compacter form comparing to the ADMM. This formation called the PPM by using the relationship (2.9) in Remark 2.8 can reduce the original five subproblems of the ADMM in Algorithm 3.1 to solve three subproblems, thus it can improve computation cost of the ADMM. Now we have rewritten (3.5a)–(3.5e) with a little variation as the following form: (3.10a)vn=argminvh1(v)+μ12∥v-d1n-1-Λ1un-1∥L22,(3.10b)zn=argminzh2(z)+μ22∥z-d2n-1-Λ2un-1∥L22,(3.10c)d1n=d1n-1+Λ1un-1-vn,(3.10d)d2n=d2n-1+Λ2un-1-zn,(3.10e)un=argminuθ2∥u-f∥L22+μ12∥Λ1u-vn+d1n∥L22+μ22∥Λ2u-zn+d2n∥L22with the first order optimality conditions given by (3.11a)vn=prox1/μ1h1(Λ1un-1+d1n-1),(3.11b)zn=prox1/μ1h2(Λ2un-1+d2n-1),(3.11c)d1n=d1n-1+Λ1un-1-vn,(3.11d)d2n=d2n-1+Λ2un-1-zn,(3.11e)θ(un-f)+μ1Λ1*(Λ1un-vn+d1n)+μ2Λ2*(Λ2un-zn+d2n)=0.If (3.11e) is replaced by
(3.12)θ(un-f)+μ1Λ1*(Λ1un-1-vn+d1n-1)+μ2Λ2*(Λ2un-1-zn+d2n-1)=0,
it follows from (3.11a)–(3.11e) and Moreau-Yosida decomposition Theorem  31.5 in [35] that
(3.13)d1n=proxμ1h1*(Λ1un-1+d1n-1),d2n=proxμ2h2*(Λ2un-1+d2n-1),un=f-(μ1θΛ1*d1n+μ2θΛ2*d2n).
So we propose the following algorithm to solve (3.1).Algorithm 3.3 (PPM for solving (3.1)).
(1) Choose the original values:d10, d20, and u0. Set θ, μ1, μ2>0 and n=1.
(2) Compute(d1n,d2n,un) by (3.13).
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Lemma 3.4.
Setx1, x2∈X and A is a maximal monotone operator, then the operators HcA and JcA satisfy
(3.14)1c2∥JcA(x1)-JcA(x2)∥L22+∥HcA(x1)-HcA(x2)∥L22≤1c2∥x1-x2∥L22.Theorem 3.5.
Assume thath1(x) and h2(u) are convex and proper. If μ1∈(0,1/θ∥Λ1∥2) and μ2∈(0,1/θ∥Λ2∥2), here ∥·∥:=max{∥Kx∥L2:x∈X with ∥x∥L2≤1} for a continuous linear operator K, then the sequence {(un,d1n,d2n)} generated by Algorithm 3.3 converges to the limit point (d-1,d-2,u-). Furthermore, the limit point u- is the solution of (1.8).


## 3.1. The Alternating Direction Method of Multipliers (ADMM)

We now consider the following constrained problem:(3.1)minu,v,wθ2∥u-f∥L22+h1(v)+h2(z),s.t.v=Λ1u,z=Λ2u,
which is clearly equivalent to the constrained problem (1.8) in the feasible set {(u,v,z):v=Λ1u,andz=Λ2u}. Throughout the following subsections, we always assume that Λ1 and Λ2 are a surjective map. It seems that the problem (3.1) including three variables looks more complex than the original unconstrained problem (1.8). In fact, this problem can be solved more easily under the condition that h1 and h2 are nondifferentiable. In the augmented Lagrangian framework [33, 34, 47, 48], the problem (3.1) equivalently solves the following minimization Lagrangian function:
(3.2)minu,v,z,ζ1,ζ2,μ1,μ2L(u,v,z,ζ1,ζ2,μ1,μ2)=θ2∥u-f∥L22+h1(v)+(ζ1,Λ1u-v)+μ12∥Λ1u-v∥L22+h2(z)+(ζ2,Λ2u-z)+μ22∥Λ2u-z∥L22,
where ζi is the Lagrangian multiplier and μi is the penalty parameter for i=1, 2. Then we can use the following augmented Lagrangian method (ALM): (3.3a)(un,vn,zn)=argminu,v,zL(u,v,z,ζ1n-1,ζ2n-1,μ1,μ2),(3.3b)ζ1n=ζ1n-1+μ1(Λ1u-v),(3.3c)ζ2n=ζ2n-1+μ2(Λ2u-z)with choosing the original values ζ10 and ζ20 to solve (3.2). If we set din=ζin/μi for i=1, 2 and omit the terms which are independent of (un,vn,zn) in (3.3a), the above strategy (3.3a)–(3.3c) can be written as (3.4a)(un,vn,zn)=argminu,v,zθ2∥u-f∥L22+h1(v)+μ12∥Λ1u-v+d1n-1∥L22+h2(z),+μ22∥Λ2u-z+d2n-1∥L22,(3.4b)d1n=d1n-1+Λ1un-vn,(3.4c)d2n=d2n-1+Λ2un-znfor the original values d10 and d20. Then we can use the following ADMM to solve (3.4a)–(3.4c).Algorithm 3.1 (ADMM for solving (3.4a)–(3.4c)).
(1) Choose the original values:v0, z0, d10, and d20. Set θ, μ1, μ2>0 and n=1.
(2) Compute(un,vn,zn,d1n,d2n) by(3.5a)
u
n
=
arg
min
u
θ
2
∥
u
-
f
∥
L
2
2
+
μ
1
2
∥
Λ
1
u
-
v
n
-
1
+
d
1
n
-
1
∥
L
2
2
+
μ
2
2
∥
Λ
2
u
-
z
n
-
1
+
d
2
n
-
1
∥
L
2
2
︸
=
h
(
u
)
,
(3.5b)
v
n
=
arg
min
v
h
1
(
v
)
+
μ
1
2
∥
Λ
1
u
n
-
v
+
d
1
n
-
1
∥
L
2
2
,
(3.5c)
z
n
=
arg
min
z
h
2
(
z
)
+
μ
2
2
∥
Λ
2
u
n
-
z
+
d
2
n
-
1
∥
L
2
2
,
(3.5d)
d
1
n
=
d
1
n
-
1
+
Λ
1
u
n
-
v
n
,
(3.5e)
d
2
n
=
d
2
n
-
1
+
Λ
2
u
n
-
z
n
.(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Sinceh(u) is differentiable and strictly convex, we can get the unique solution of (3.5a) which satisfies (3.6a)θ(un-f)+μ1Λ1*(Λ1un-vn-1+d1n-1)+μ2Λ2*(Λ2un-zn-1+d2n-1)=0,(3.6b)0=w1n+μ1(vn-d1n-1-Λ1un),(3.6c)0=w2n+μ2(zn-d2n-1-Λ2un),(3.6d)d1n=d1n-1+Λ1un-vn,(3.6e)d2n=d2n-1+Λ2un-zn,where Λ1* and Λ2* are the adjoint operators of Λ1 and Λ2, w1n∈∂h1(vn) and w2n∈∂h2(zn), respectively. It follows that the solution un in (3.6a) can be directly obtained by the following explicit formulation:
(3.7)un=(θI+μ1Λ1*Λ1+μ2Λ2*Λ2)︸M-1(θf+μ1Λ1*(vn-1-d1n-1)+μ2Λ2*(zn-1-d2n-1)).
However, when the operator M is ill-posed, the solution is unsuitable or unavailable. Hence we have to go back to (3.6a) and to employ some iteration strategy such as the Gauss-Seidel method to solve this equation. On the other hand, it is obvious that (3.5b) and (3.5c) can be looked at as the proximal mapping, so the solutions of the minimization problems (3.5b) and (3.5c) can be obviously written as
(3.8)vn=Prox1/μ1h1(Λ1un+d1n-1),zn=Prox1/μ2h2(Λ2un+d2n-1).Theorem 3.2.
Assume that(u*,v*,z*,w1*,w2*) is the saddle point of the Lagrange function
(3.9)ℒ(u,v,z,w1,w2)=θ2∥u-f∥L22+h1(v)+(w1,Λ1u-v)+h2(z)+(w2,Λ2u-z).
Then u* is the solution of the minimization problem (1.8). Furthermore, the sequence {(un,vn,zn)} generated by Algorithm 3.1 converges to (u*,v*,z*).Notice that Algorithm3.1 can be actually looked at as the split Bregman method [25]. The based idea of this method is to introduce some intermediate variables so as to transform the original problem into some subproblems which are easily solved. The connection between the split Bregman method and the ADMM has been shown in [29, 49]. However, our algorithm considers the sum of three convex functions, which is more general than the related algorithms in [25, 49]. Furthermore, it must be noted that v and z are completely separated in (3.4a), so the two subproblems (3.5b) and (3.5c) are parallel. Therefore the convergence results of the ADMM can be applied here.


## 3.2. The Proximal Point Method

Though the ADMM in Algorithm3.1 can effectively solve the original problem (3.1), we have to solve five subproblems. This actually makes this method suffer from a looser form as in [25, 45] so that it can badly affect its numerical computation efficiency. In this subsection, we propose a compacter form comparing to the ADMM. This formation called the PPM by using the relationship (2.9) in Remark 2.8 can reduce the original five subproblems of the ADMM in Algorithm 3.1 to solve three subproblems, thus it can improve computation cost of the ADMM. Now we have rewritten (3.5a)–(3.5e) with a little variation as the following form: (3.10a)vn=argminvh1(v)+μ12∥v-d1n-1-Λ1un-1∥L22,(3.10b)zn=argminzh2(z)+μ22∥z-d2n-1-Λ2un-1∥L22,(3.10c)d1n=d1n-1+Λ1un-1-vn,(3.10d)d2n=d2n-1+Λ2un-1-zn,(3.10e)un=argminuθ2∥u-f∥L22+μ12∥Λ1u-vn+d1n∥L22+μ22∥Λ2u-zn+d2n∥L22with the first order optimality conditions given by (3.11a)vn=prox1/μ1h1(Λ1un-1+d1n-1),(3.11b)zn=prox1/μ1h2(Λ2un-1+d2n-1),(3.11c)d1n=d1n-1+Λ1un-1-vn,(3.11d)d2n=d2n-1+Λ2un-1-zn,(3.11e)θ(un-f)+μ1Λ1*(Λ1un-vn+d1n)+μ2Λ2*(Λ2un-zn+d2n)=0.If (3.11e) is replaced by
(3.12)θ(un-f)+μ1Λ1*(Λ1un-1-vn+d1n-1)+μ2Λ2*(Λ2un-1-zn+d2n-1)=0,
it follows from (3.11a)–(3.11e) and Moreau-Yosida decomposition Theorem  31.5 in [35] that
(3.13)d1n=proxμ1h1*(Λ1un-1+d1n-1),d2n=proxμ2h2*(Λ2un-1+d2n-1),un=f-(μ1θΛ1*d1n+μ2θΛ2*d2n).
So we propose the following algorithm to solve (3.1).Algorithm 3.3 (PPM for solving (3.1)).
(1) Choose the original values:d10, d20, and u0. Set θ, μ1, μ2>0 and n=1.
(2) Compute(d1n,d2n,un) by (3.13).
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Lemma 3.4.
Setx1, x2∈X and A is a maximal monotone operator, then the operators HcA and JcA satisfy
(3.14)1c2∥JcA(x1)-JcA(x2)∥L22+∥HcA(x1)-HcA(x2)∥L22≤1c2∥x1-x2∥L22.Theorem 3.5.
Assume thath1(x) and h2(u) are convex and proper. If μ1∈(0,1/θ∥Λ1∥2) and μ2∈(0,1/θ∥Λ2∥2), here ∥·∥:=max{∥Kx∥L2:x∈X with ∥x∥L2≤1} for a continuous linear operator K, then the sequence {(un,d1n,d2n)} generated by Algorithm 3.3 converges to the limit point (d-1,d-2,u-). Furthermore, the limit point u- is the solution of (1.8).


## 4. Some Applications in Image Restoration

In Section4.1, we apply the ADMM and the above PPM to the image denoising problem. Here we also compare the proposed hybrid model with the ROF model and the LLT model. Then, based on the proposed hybrid model, we set the PPM as a basic method to solve image deblurring, image inpainting, and image denoising for the multiplicative noise in the last three subsections. For simplicity, we assume that the image region Ω is squared with the size M×M and set S=RM×M, T=S×S, and Z=T×T as in [17]. The usual scalar product can be denoted as 〈p1,p2〉T:=∑i=1M‍∑j=1M‍pi,j1pi,j2 for p1,p2∈T and (p,q)Z=∑i,j=1M‍(pi,j11qi,j11+pi,j12qi,j12+pi,j21qi,j21+pi,j22qi,j22) for p,q∈Z. The L1 norm of p=(p1,p2)∈T is defined by |p|ℓ1=p12+p22 and the ℓ1 norm of q=(q1,q2;q3,q4)∈Z is defined by |q|ℓ1=q12+q22+q32+q42. If u∈S, we use ∇+=(∇x+,∇y+)∈T to denote the first order forward difference operator with
(4.1)∇x+ui,j={ui+1,j-ui,jfor1≤i<N,0fori=N,∇y+ui,j={ui,j+1-ui,jfor1≤j<N,0forj=N,
and use ∇-=(∇x-,∇y-) to denote the first order backward difference operator with
(4.2)∇x-ui,j={-ui,j,fori=1,ui,j-ui-1,jfor1<i<N,u(i,j)fori=N,∇y-ui,j={-ui,jforj=1,ui,j-ui,j-1for1<j<N,ui,jforj=N,
for i,j=1,…,N. Based on the first order difference operators, we can give the second order difference operator as follows:
(4.3)∇2ui,j=(∇x-(∇x+ui,j)∇x+(∇y+ui,j)∇y-(∇x-ui,j)∇y-(∇y+ui,j))∈Z.
Using the same approach, we can define some other second order operators such as ∇x+(∇x-ui,j), ∇x-(∇y-ui,j), ∇y+(∇x+ui,j), and ∇y+(∇y-ui,j). Then we can give the first order and the second divergence operators as
(4.4)divui,j=∇x-ui,j+∇y-ui,j,div2ui,j=∇x+(∇x-ui,j)+∇x-(∇y-ui,j)+∇y+(∇x+ui,j)+∇y+(∇y-ui,j).
Furthermore, if we set p∈T and q∈Z, it is easy to deduce that
(4.5)∥divp∥22≤8∥p∥T2,∥div2q∥22≤64∥q∥Z2.Remark 4.1.
The related examples in the following subsections are performed using Windows 7 and Matlab 2009(a) on a desktop with Intel Core i5 processor at 2.4 GHz and 4 GB memory. All of the parameters for related models are chosen by trial and empirically which can yield better restored images. On the other hand, we should notice that it is not very expensive when we use the ADMM and the PPM to getun, but the total computational effort of one outer iteration requiring many inner steps can be very huge. In order to reduce the computational effort and keep fair comparison of these two methods, we so simplify the inner-outer iterative framework by performing only one-step in inner iteration. It is obvious that these sets are very efficient from the following numerical experiences.


### 4.1. Image Denoising for the Additive Noise

In this subsection, we consider to use the ADMM and the PPM to solve (1.7) for restoring the additive noisy image. If we set Λ1=∇ and Λ2=∇2, then the algorithms are proposed as follows.Algorithm 4.2 (ADMM to solve (1.7)).
(1) Choose the originald10=v0=0∈T and d20=z0=0∈Z. Set θ, μ1, μ2>0 and n=1.
(2) Compute(un,vn,zn,d1n,d2n) by (4.6a)(θI-μ1Δ+μ2Δ2)un=θf+μ1div(vn-1-d1n-1)-μ2div2(zn-1-d2n-1)(4.6b)vn=∇un+d1n-1|∇un+d1n-1|ℓ1·max(|∇un+d1n-1|ℓ1-(1-g(x))μ1,0),(4.6c)zn=∇2un+d2n-1|∇2un+d2n-1|ℓ1·max(|∇2un+d2n-1|ℓ1-g(x)μ2,0),(4.6d)d1n=d1n-1+∇un-vn,(4.6e)d2n=d2n-1+∇2un-zn,where Δ=div∘∇ and Δ2=div2∘∇2.
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Remark 4.3.
For the first subproblem (4.6a), we can use the Gauss-Seidel method as shown in [25] to get the solution. However, in this paper, we use the following strategy:
(4.7)ui,jn=θfi,j+μ1div(vi,jn-1-d1,i,jn-1)+4μ1ui,jn-1-μ2div2(zi,jn-1-d2,i,jn-1)+28μ2ui,jn-1θ+4μ1+28μ2,
where some information of operators Δ and Δ2 related to u is used. The formulas (4.6b) and (4.6c) of Algorithm 4.2 can be easily deduced from
(4.8)minv|(1-g)v|ℓ1+μ12∥v-∇un-d1n-1∥22,minz|gv|ℓ1+μ22∥z-∇2un-d2n-1∥22.
Furthermore, following form Theorem 3.2, we can also deduce that the sequence {un} generated by Algorithm 4.2 converges to the solution of (1.7).Algorithm 4.4 (PPM to solve (1.7)).
(1) Choose the originald10=0∈S, d20=0∈Z and u0=f. Set θ, μ1, μ2>0 and n=1.
(2) Compute(d1n,d2n,un) by
(4.9)d1n=∇un-1+d1n-1|∇un-1+d1n-1|ℓ1·min(|∇un-1+d1n-1|ℓ1,1),d2n=∇2un-1+d2n-1|∇2un-1+d2n-1|ℓ1·min(|∇2un-1+d2n-1|ℓ1,1),un=f-μ1θ(1-g)divd1n+μ2θgdiv2d2n.
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).For Algorithm4.4, based on the relations in (4.5) and Theorem 3.5, we have the following result.Theorem 4.5.
Ifμ1∈(0,1/8θ) and μ2∈(0,1/64θ), then the sequence {un} generated by Algorithm 4.4 converges to the solution of (1.7).Remark 4.6.
For the above two algorithms, we can also setg as a constant between 0 and 1. In fact, it is easy to find that the algorithms correspond to solving the ROF model or the LLT model, respectively, when g=0 or g=1. At this time, the iteration strategy can be simplified as
(4.10)ui,jn=θfi,j+μ1div(vi,jn-1-d1,i,jn-1)+4μ1ui,jn-1θ+4μ1,ui,jn=θfi,j-μ2div2(zi,jn-1-d2,i,jn-1)+28μ2ui,jn-1θ+28μ2
for these two models, respectively. Furthermore, when g∈(0,1), these two algorithms correspond to solve the model which is the convex combination of the ROF model and the LLT model.Example 4.7.
In this example, we compare the ADMM with the PPM for solving the ROF model (1.4), the LLT model (1.6), and the hybrid model (1.7). The original images with three different sizes shown in Figure 1 are added to the Gaussian white noise with the standard deviation σ=15.3. Iterations were terminated when the stop conditions ∥un+1-un∥2/∥un∥2≤ε are met. It is easy to find the related results from Table 1 that the PPM is faster than the ADMM. Especially, the average CPU time of the PPM compared with that of the ADMM can save about 50% for the ROF model and the LLT model. It saves about 40% for the hybrid model.Table 1
The related results in Example4.7.
g
=
0
The ROF model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
1
28
×
128
1.0
×
1
0
-
4
0.2184
86
24.7949
0.1248
83
24.7043
256
×
256
1.0
×
1
0
-
4
0.9204
77
16.2934
0.2652
38
16.3540
512
×
512
1.0
×
1
0
-
4
4.3368
64
15.7146
1.9656
53
15.9353
g
=
1
The LLT model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
128
×
128
1.0
×
1
0
-
4
0.4368
60
24.5688
0.3432
80
25.0287
256
×
256
1.0
×
1
0
-
4
1.7004
55
16.9065
0.9360
53
16.9083
512
×
512
1.0
×
1
0
-
4
15.4441
55
15.9432
7.2384
50
15.9512
g
∈
(
0,1
)
The hybrid model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
128
×
128
1.0
×
1
0
-
4
0.7176
69
25.7224
0.5304
70
25.4982
256
×
256
1.0
×
1
0
-
4
2.5584
55
16.9699
1.4664
52
16.9701
512
×
512
1.0
×
1
0
-
4
19.6561
56
15.9503
9.3445
50
15.9512The original images and the noisy images with three different sizes in Example4.7.
(a)
(b)
(c)
(d)
(e)
(f)Example 4.8.
In this example, the noisy image is added to the Gaussian white noisy with the standard deviationσ=12. The algorithms will be stopped after 100 iterations. We compare the results generated by the ROF model, the LLT model, the convex combination of the ROF model and the hybrid model. As we can see from Table 2, the hybrid model get the lowest MSE and the highest SNR; these imply that the hybrid model can give the best restored image. On the other hand, it is easy to find that the ROF model makes staircasing effect appear and the LLT model leads to edge blurring. In fact, they are based on the fact that the restored model by the ROF model is piecewise constant on large areas and the LLT model as a higher model damps oscillations much faster in the region of edges. For the convex combined model and the hybrid model, they can efficiently suppress these two drawbacks. Furthermore, the hybrid model is more efficient than the convex combined model, because the hybrid model uses the edge detector function which can efficiently coordinate edge information. It should be noticed that here we use the Chambolle’s strategy [17] to solve the convex combined model so that it is slower than the strategy by using the PPM to solve the hybrid model. To focus on these facts, we present some zoomed-in local results and select a slice of the images which meets contours and the smooth regions shown in Figures 2 and 3.Table 2
The related data in Example4.8 (here R.P. is the regularization parameter).
Model
R.P.
Time (s)
SNR
MSE
ROF
4.5
0.8268
17.7677
44.0400
LLT
3.0
2.0904
17.9998
42.7826
Convex
3.5
2.6832
18.2399
39.7756
Hybrid
4.5
2.5272
18.3694
39.1447The original and the noisy image in Example4.8.
(a)
(b)The related restored images in Example4.8.
(a)
(b)
(c)
(d)


### 4.2. Other Applications

In this subsection, we extend the hybrid model to other classes of restoration problems. As we can see in Section4.1, the hybrid model has some advantages compared with the ROF model and the LLT model. Since the PPM is faster and more efficient than the ADMM, we only employ the PPM to solve the related image restoration model. Now we first consider the following iteration strategy: (4.11a)zn=ℱ(un)=:argminzβ2∥un-z∥L22+h(z),(4.11b)un+1=𝒢(zn):=argminuβ2∥u-zn∥L22+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx︸𝒟(u),where h(z)∈C1(Ω). It is easy to find that the minimization problem 𝒟(u) in (4.11b) is coercive and strictly convex, so the subproblem (4.11b) has a unique solution. Based on [50, Lemma 2.4], we also deduce that the operator 𝒢 is (1/2)-averaged nonexpansive.Theorem 4.9.
Assume that the functional(4.12)ℳ(u,z)=β2∥u-z∥L22+h(z)+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx
is convex, coercive, and bounded below, then the sequence {(un,zn)} generated by (4.11a)-(4.11b) converges to a point (u*,z*). Furthermore, when u*=z*, u* is the solution of
(4.13)minuh(u)+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx.Remark 4.10.
It should be noticed that the following three models satisfy the conditions of Theorem4.9, so we can use the strategy (4.11a)-(4.11b) to solve these models.


#### 4.2.1. Image Deblurring

Now we apply the hybrid model to the image deblurring problems with the following formula:(4.14)minuλ2∥Ku-u0∥L22+|u|BV(1-g)+|u|BVg2.
Because of the compactness, the blurring operator K is not continuously invertible, which implies that we cannot directly use the proximal method to solve minimization problem (4.14). However, based on the penalty method [48], we can transform the problem (4.14) to solve the following problem: (4.15a)zn=(μλKTK+I)-1(μλKTf+un),(4.15b)un+1:=argminu12μ∥u-zn∥L22+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx.The key of (4.15a)-(4.15b) is to separate the blurring operator K from the original problem (4.14) so that we can avoid the ill-posed operator K. Furthermore, it is obvious that the problem (4.15b) satisfie Theorem 4.9, so we can use the proximal method to solve it.Example 4.11.
In this experiment, we use the image Lena, which is blurred with a Gaussian kernel of “hsize = 3” and in which is added the Gaussian white noise with the standard deviationσ=0.02. The related images and SNRs are arranged in Figure 4. As we can see in Example 4.7, the proximal method tends to be stable before 100 iterations, so here we fix the algorithm to be stopped when the iteration attains 100. It is easy to find that the hybrid model can give a better image quality than the other two models. Especially, we also observe some staircasing effect in Figure 4(b) and edge blurring in Figure 4(c). However, all of the drawbacks can be efficiently suppressed in Figure 1(d).The related SNRs in Example4.11. (a) 11.2166; (b) 14.3378; (c) 14.4778; (d) 14.5118.
(a)
Original image
(b)
ROF model
(c)
LLT model
(d)
Hybrid model


#### 4.2.2. Image Inpainting

Here we consider to use the hybrid model for the image inpainting problem with the following form:(4.16)minuλ-2∥u-u0∥L2(Ω∖D)2+|(1-g)u|BV+|gu|BV2,
where D is the inpainting domain and λ- satisfies
(4.17)λ-={λ,ifu∈Ω∖D,0,ifu∈D.
If we introduce an auxiliary variable z, based on the penalty method [48] again, then the solution of the problem (4.15a)-(4.15b) can be approximated by solving the following problem: (4.18a)un:=zn+λ-μu01+λ-μ,(4.18b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.Example 4.12.
In this example, we show the results of real image inpainting in Figures6 and 7. We consider to use these three models to restore images shown in Figures 5(c) and 5(d), which are, respectively, contaminated by the mask image and the noise with the standard deviation σ=12. The parameters of these three models are chosen for getting better restored images and the algorithms will be stopped after 500 iterations (see also Table 3). It is obvious that these three models can efficiently restore the original deteriorated images. However, there are much more inappropriate information restored by the ROF model than the LLT model and the hybrid model. These advantages is based on the fact that the fourth-order linear diffusion (the LLT model and the hybrid model) damps oscillations much faster than second-order diffusion (the ROF model). Especially, these unsuitable effects can be easily seen from the local zooming images in Figures 6 and 7. Furthermore, the restored image by the hybrid model looks more natural than the LLT model.Table 3
The related results in Example4.12.
Model
Image (c) in Figure5
Image (d) in Figure5
λ
μ
SNR
λ
μ
SNR
ROF
220
0.005
18.5533
32
0.008
15.4055
LLT
150
0.005
16.2934
60
0.005
16.4336
Hybrid
220
0.005
20.2713
70
0.0025
16.6034The related images in Example4.12.
(a)
Original image
(b)
Mask image
(c)
Painting image
(d)
Noisy imageThe related restored images corresponding to the painting image in Example4.12.The related restored images corresponding to the noisy image in Example4.12.


#### 4.2.3. Image Denoising for the Multiplicative Noise

Based on the hybrid model (1.7), the multiplicative noise model (1.5) can be naturally extended to the following minimization:
(4.19)minzα∫Ω‍(z+fe-z)dx+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.
It is obvious that the problem (4.19) can be approximated by(4.20a)un:=argminu12μ∥zn-u∥L22+α∫Ω‍(u+fe-u)dx,(4.20b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.For the first subproblem (4.20a), its solution un can be approximately obtained by using the Newton method to solve the following nonlinear equation:
(4.21)μ(u-zn)+α(1-fe-u)=0.Example 4.13.
In this example, we restore the multiplicative noisy image. The noisy Lena image shown in Figure7 is contaminated by the Gamma noise (L=33) with mean one, in which probability density function p(u) is given by
(4.22)p(u)={LLuL-1Γ(L)e-Lu,ifu>0,0,ifu≤0,
where L is an integer and Γ(·) is a Gamma function. Based on the iteration formula (4.20a)-(4.20b), we should notice that there are two interior iterations. For employing the Newton method to solve the problem (4.21), we set the stepsize t=1 and the Newton method will be stopped when ∥un+1-un∥ℓ2/∥un∥ℓ2≤1.0×10-5. For solving the second subproblem (4.20b), we set the fixed iterations with 40. In fact for using the PPM to solve the ROF model, the LLT model and the hybrid model, the solutions of these models will tend to the steady. The outer iteration will be stopped after 200 iterations. The related restored images are shown in Figure 8, and it is easy to find that the hybrid model gives a better restored image than the two other models.Figure 8
The related images and the local zooming images in Example4.13. SNRs: (a) 9.2017; (b) 15.6059; (c) 15.1754; (d) 16.0546.


## 4.1. Image Denoising for the Additive Noise

In this subsection, we consider to use the ADMM and the PPM to solve (1.7) for restoring the additive noisy image. If we set Λ1=∇ and Λ2=∇2, then the algorithms are proposed as follows.Algorithm 4.2 (ADMM to solve (1.7)).
(1) Choose the originald10=v0=0∈T and d20=z0=0∈Z. Set θ, μ1, μ2>0 and n=1.
(2) Compute(un,vn,zn,d1n,d2n) by (4.6a)(θI-μ1Δ+μ2Δ2)un=θf+μ1div(vn-1-d1n-1)-μ2div2(zn-1-d2n-1)(4.6b)vn=∇un+d1n-1|∇un+d1n-1|ℓ1·max(|∇un+d1n-1|ℓ1-(1-g(x))μ1,0),(4.6c)zn=∇2un+d2n-1|∇2un+d2n-1|ℓ1·max(|∇2un+d2n-1|ℓ1-g(x)μ2,0),(4.6d)d1n=d1n-1+∇un-vn,(4.6e)d2n=d2n-1+∇2un-zn,where Δ=div∘∇ and Δ2=div2∘∇2.
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).Remark 4.3.
For the first subproblem (4.6a), we can use the Gauss-Seidel method as shown in [25] to get the solution. However, in this paper, we use the following strategy:
(4.7)ui,jn=θfi,j+μ1div(vi,jn-1-d1,i,jn-1)+4μ1ui,jn-1-μ2div2(zi,jn-1-d2,i,jn-1)+28μ2ui,jn-1θ+4μ1+28μ2,
where some information of operators Δ and Δ2 related to u is used. The formulas (4.6b) and (4.6c) of Algorithm 4.2 can be easily deduced from
(4.8)minv|(1-g)v|ℓ1+μ12∥v-∇un-d1n-1∥22,minz|gv|ℓ1+μ22∥z-∇2un-d2n-1∥22.
Furthermore, following form Theorem 3.2, we can also deduce that the sequence {un} generated by Algorithm 4.2 converges to the solution of (1.7).Algorithm 4.4 (PPM to solve (1.7)).
(1) Choose the originald10=0∈S, d20=0∈Z and u0=f. Set θ, μ1, μ2>0 and n=1.
(2) Compute(d1n,d2n,un) by
(4.9)d1n=∇un-1+d1n-1|∇un-1+d1n-1|ℓ1·min(|∇un-1+d1n-1|ℓ1,1),d2n=∇2un-1+d2n-1|∇2un-1+d2n-1|ℓ1·min(|∇2un-1+d2n-1|ℓ1,1),un=f-μ1θ(1-g)divd1n+μ2θgdiv2d2n.
(3) If the stop criterion is not satisfied, setn:=n+1 and go to step (2).For Algorithm4.4, based on the relations in (4.5) and Theorem 3.5, we have the following result.Theorem 4.5.
Ifμ1∈(0,1/8θ) and μ2∈(0,1/64θ), then the sequence {un} generated by Algorithm 4.4 converges to the solution of (1.7).Remark 4.6.
For the above two algorithms, we can also setg as a constant between 0 and 1. In fact, it is easy to find that the algorithms correspond to solving the ROF model or the LLT model, respectively, when g=0 or g=1. At this time, the iteration strategy can be simplified as
(4.10)ui,jn=θfi,j+μ1div(vi,jn-1-d1,i,jn-1)+4μ1ui,jn-1θ+4μ1,ui,jn=θfi,j-μ2div2(zi,jn-1-d2,i,jn-1)+28μ2ui,jn-1θ+28μ2
for these two models, respectively. Furthermore, when g∈(0,1), these two algorithms correspond to solve the model which is the convex combination of the ROF model and the LLT model.Example 4.7.
In this example, we compare the ADMM with the PPM for solving the ROF model (1.4), the LLT model (1.6), and the hybrid model (1.7). The original images with three different sizes shown in Figure 1 are added to the Gaussian white noise with the standard deviation σ=15.3. Iterations were terminated when the stop conditions ∥un+1-un∥2/∥un∥2≤ε are met. It is easy to find the related results from Table 1 that the PPM is faster than the ADMM. Especially, the average CPU time of the PPM compared with that of the ADMM can save about 50% for the ROF model and the LLT model. It saves about 40% for the hybrid model.Table 1
The related results in Example4.7.
g
=
0
The ROF model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
1
28
×
128
1.0
×
1
0
-
4
0.2184
86
24.7949
0.1248
83
24.7043
256
×
256
1.0
×
1
0
-
4
0.9204
77
16.2934
0.2652
38
16.3540
512
×
512
1.0
×
1
0
-
4
4.3368
64
15.7146
1.9656
53
15.9353
g
=
1
The LLT model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
128
×
128
1.0
×
1
0
-
4
0.4368
60
24.5688
0.3432
80
25.0287
256
×
256
1.0
×
1
0
-
4
1.7004
55
16.9065
0.9360
53
16.9083
512
×
512
1.0
×
1
0
-
4
15.4441
55
15.9432
7.2384
50
15.9512
g
∈
(
0,1
)
The hybrid model
Stopping
ADMM
PPM
Size
Conditionsε
Time (s)
Ite.
SNR
Time (s)
Ite.
SNR
128
×
128
1.0
×
1
0
-
4
0.7176
69
25.7224
0.5304
70
25.4982
256
×
256
1.0
×
1
0
-
4
2.5584
55
16.9699
1.4664
52
16.9701
512
×
512
1.0
×
1
0
-
4
19.6561
56
15.9503
9.3445
50
15.9512The original images and the noisy images with three different sizes in Example4.7.
(a)
(b)
(c)
(d)
(e)
(f)Example 4.8.
In this example, the noisy image is added to the Gaussian white noisy with the standard deviationσ=12. The algorithms will be stopped after 100 iterations. We compare the results generated by the ROF model, the LLT model, the convex combination of the ROF model and the hybrid model. As we can see from Table 2, the hybrid model get the lowest MSE and the highest SNR; these imply that the hybrid model can give the best restored image. On the other hand, it is easy to find that the ROF model makes staircasing effect appear and the LLT model leads to edge blurring. In fact, they are based on the fact that the restored model by the ROF model is piecewise constant on large areas and the LLT model as a higher model damps oscillations much faster in the region of edges. For the convex combined model and the hybrid model, they can efficiently suppress these two drawbacks. Furthermore, the hybrid model is more efficient than the convex combined model, because the hybrid model uses the edge detector function which can efficiently coordinate edge information. It should be noticed that here we use the Chambolle’s strategy [17] to solve the convex combined model so that it is slower than the strategy by using the PPM to solve the hybrid model. To focus on these facts, we present some zoomed-in local results and select a slice of the images which meets contours and the smooth regions shown in Figures 2 and 3.Table 2
The related data in Example4.8 (here R.P. is the regularization parameter).
Model
R.P.
Time (s)
SNR
MSE
ROF
4.5
0.8268
17.7677
44.0400
LLT
3.0
2.0904
17.9998
42.7826
Convex
3.5
2.6832
18.2399
39.7756
Hybrid
4.5
2.5272
18.3694
39.1447The original and the noisy image in Example4.8.
(a)
(b)The related restored images in Example4.8.
(a)
(b)
(c)
(d)


## 4.2. Other Applications

In this subsection, we extend the hybrid model to other classes of restoration problems. As we can see in Section4.1, the hybrid model has some advantages compared with the ROF model and the LLT model. Since the PPM is faster and more efficient than the ADMM, we only employ the PPM to solve the related image restoration model. Now we first consider the following iteration strategy: (4.11a)zn=ℱ(un)=:argminzβ2∥un-z∥L22+h(z),(4.11b)un+1=𝒢(zn):=argminuβ2∥u-zn∥L22+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx︸𝒟(u),where h(z)∈C1(Ω). It is easy to find that the minimization problem 𝒟(u) in (4.11b) is coercive and strictly convex, so the subproblem (4.11b) has a unique solution. Based on [50, Lemma 2.4], we also deduce that the operator 𝒢 is (1/2)-averaged nonexpansive.Theorem 4.9.
Assume that the functional(4.12)ℳ(u,z)=β2∥u-z∥L22+h(z)+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx
is convex, coercive, and bounded below, then the sequence {(un,zn)} generated by (4.11a)-(4.11b) converges to a point (u*,z*). Furthermore, when u*=z*, u* is the solution of
(4.13)minuh(u)+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx.Remark 4.10.
It should be noticed that the following three models satisfy the conditions of Theorem4.9, so we can use the strategy (4.11a)-(4.11b) to solve these models.


### 4.2.1. Image Deblurring

Now we apply the hybrid model to the image deblurring problems with the following formula:(4.14)minuλ2∥Ku-u0∥L22+|u|BV(1-g)+|u|BVg2.
Because of the compactness, the blurring operator K is not continuously invertible, which implies that we cannot directly use the proximal method to solve minimization problem (4.14). However, based on the penalty method [48], we can transform the problem (4.14) to solve the following problem: (4.15a)zn=(μλKTK+I)-1(μλKTf+un),(4.15b)un+1:=argminu12μ∥u-zn∥L22+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx.The key of (4.15a)-(4.15b) is to separate the blurring operator K from the original problem (4.14) so that we can avoid the ill-posed operator K. Furthermore, it is obvious that the problem (4.15b) satisfie Theorem 4.9, so we can use the proximal method to solve it.Example 4.11.
In this experiment, we use the image Lena, which is blurred with a Gaussian kernel of “hsize = 3” and in which is added the Gaussian white noise with the standard deviationσ=0.02. The related images and SNRs are arranged in Figure 4. As we can see in Example 4.7, the proximal method tends to be stable before 100 iterations, so here we fix the algorithm to be stopped when the iteration attains 100. It is easy to find that the hybrid model can give a better image quality than the other two models. Especially, we also observe some staircasing effect in Figure 4(b) and edge blurring in Figure 4(c). However, all of the drawbacks can be efficiently suppressed in Figure 1(d).The related SNRs in Example4.11. (a) 11.2166; (b) 14.3378; (c) 14.4778; (d) 14.5118.
(a)
Original image
(b)
ROF model
(c)
LLT model
(d)
Hybrid model


### 4.2.2. Image Inpainting

Here we consider to use the hybrid model for the image inpainting problem with the following form:(4.16)minuλ-2∥u-u0∥L2(Ω∖D)2+|(1-g)u|BV+|gu|BV2,
where D is the inpainting domain and λ- satisfies
(4.17)λ-={λ,ifu∈Ω∖D,0,ifu∈D.
If we introduce an auxiliary variable z, based on the penalty method [48] again, then the solution of the problem (4.15a)-(4.15b) can be approximated by solving the following problem: (4.18a)un:=zn+λ-μu01+λ-μ,(4.18b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.Example 4.12.
In this example, we show the results of real image inpainting in Figures6 and 7. We consider to use these three models to restore images shown in Figures 5(c) and 5(d), which are, respectively, contaminated by the mask image and the noise with the standard deviation σ=12. The parameters of these three models are chosen for getting better restored images and the algorithms will be stopped after 500 iterations (see also Table 3). It is obvious that these three models can efficiently restore the original deteriorated images. However, there are much more inappropriate information restored by the ROF model than the LLT model and the hybrid model. These advantages is based on the fact that the fourth-order linear diffusion (the LLT model and the hybrid model) damps oscillations much faster than second-order diffusion (the ROF model). Especially, these unsuitable effects can be easily seen from the local zooming images in Figures 6 and 7. Furthermore, the restored image by the hybrid model looks more natural than the LLT model.Table 3
The related results in Example4.12.
Model
Image (c) in Figure5
Image (d) in Figure5
λ
μ
SNR
λ
μ
SNR
ROF
220
0.005
18.5533
32
0.008
15.4055
LLT
150
0.005
16.2934
60
0.005
16.4336
Hybrid
220
0.005
20.2713
70
0.0025
16.6034The related images in Example4.12.
(a)
Original image
(b)
Mask image
(c)
Painting image
(d)
Noisy imageThe related restored images corresponding to the painting image in Example4.12.The related restored images corresponding to the noisy image in Example4.12.


### 4.2.3. Image Denoising for the Multiplicative Noise

Based on the hybrid model (1.7), the multiplicative noise model (1.5) can be naturally extended to the following minimization:
(4.19)minzα∫Ω‍(z+fe-z)dx+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.
It is obvious that the problem (4.19) can be approximated by(4.20a)un:=argminu12μ∥zn-u∥L22+α∫Ω‍(u+fe-u)dx,(4.20b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.For the first subproblem (4.20a), its solution un can be approximately obtained by using the Newton method to solve the following nonlinear equation:
(4.21)μ(u-zn)+α(1-fe-u)=0.Example 4.13.
In this example, we restore the multiplicative noisy image. The noisy Lena image shown in Figure7 is contaminated by the Gamma noise (L=33) with mean one, in which probability density function p(u) is given by
(4.22)p(u)={LLuL-1Γ(L)e-Lu,ifu>0,0,ifu≤0,
where L is an integer and Γ(·) is a Gamma function. Based on the iteration formula (4.20a)-(4.20b), we should notice that there are two interior iterations. For employing the Newton method to solve the problem (4.21), we set the stepsize t=1 and the Newton method will be stopped when ∥un+1-un∥ℓ2/∥un∥ℓ2≤1.0×10-5. For solving the second subproblem (4.20b), we set the fixed iterations with 40. In fact for using the PPM to solve the ROF model, the LLT model and the hybrid model, the solutions of these models will tend to the steady. The outer iteration will be stopped after 200 iterations. The related restored images are shown in Figure 8, and it is easy to find that the hybrid model gives a better restored image than the two other models.Figure 8
The related images and the local zooming images in Example4.13. SNRs: (a) 9.2017; (b) 15.6059; (c) 15.1754; (d) 16.0546.


## 4.2.1. Image Deblurring

Now we apply the hybrid model to the image deblurring problems with the following formula:(4.14)minuλ2∥Ku-u0∥L22+|u|BV(1-g)+|u|BVg2.
Because of the compactness, the blurring operator K is not continuously invertible, which implies that we cannot directly use the proximal method to solve minimization problem (4.14). However, based on the penalty method [48], we can transform the problem (4.14) to solve the following problem: (4.15a)zn=(μλKTK+I)-1(μλKTf+un),(4.15b)un+1:=argminu12μ∥u-zn∥L22+∫Ω‍|(1-g)∇u|dx+∫Ω‍|g∇2u|dx.The key of (4.15a)-(4.15b) is to separate the blurring operator K from the original problem (4.14) so that we can avoid the ill-posed operator K. Furthermore, it is obvious that the problem (4.15b) satisfie Theorem 4.9, so we can use the proximal method to solve it.Example 4.11.
In this experiment, we use the image Lena, which is blurred with a Gaussian kernel of “hsize = 3” and in which is added the Gaussian white noise with the standard deviationσ=0.02. The related images and SNRs are arranged in Figure 4. As we can see in Example 4.7, the proximal method tends to be stable before 100 iterations, so here we fix the algorithm to be stopped when the iteration attains 100. It is easy to find that the hybrid model can give a better image quality than the other two models. Especially, we also observe some staircasing effect in Figure 4(b) and edge blurring in Figure 4(c). However, all of the drawbacks can be efficiently suppressed in Figure 1(d).The related SNRs in Example4.11. (a) 11.2166; (b) 14.3378; (c) 14.4778; (d) 14.5118.
(a)
Original image
(b)
ROF model
(c)
LLT model
(d)
Hybrid model


## 4.2.2. Image Inpainting

Here we consider to use the hybrid model for the image inpainting problem with the following form:(4.16)minuλ-2∥u-u0∥L2(Ω∖D)2+|(1-g)u|BV+|gu|BV2,
where D is the inpainting domain and λ- satisfies
(4.17)λ-={λ,ifu∈Ω∖D,0,ifu∈D.
If we introduce an auxiliary variable z, based on the penalty method [48] again, then the solution of the problem (4.15a)-(4.15b) can be approximated by solving the following problem: (4.18a)un:=zn+λ-μu01+λ-μ,(4.18b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.Example 4.12.
In this example, we show the results of real image inpainting in Figures6 and 7. We consider to use these three models to restore images shown in Figures 5(c) and 5(d), which are, respectively, contaminated by the mask image and the noise with the standard deviation σ=12. The parameters of these three models are chosen for getting better restored images and the algorithms will be stopped after 500 iterations (see also Table 3). It is obvious that these three models can efficiently restore the original deteriorated images. However, there are much more inappropriate information restored by the ROF model than the LLT model and the hybrid model. These advantages is based on the fact that the fourth-order linear diffusion (the LLT model and the hybrid model) damps oscillations much faster than second-order diffusion (the ROF model). Especially, these unsuitable effects can be easily seen from the local zooming images in Figures 6 and 7. Furthermore, the restored image by the hybrid model looks more natural than the LLT model.Table 3
The related results in Example4.12.
Model
Image (c) in Figure5
Image (d) in Figure5
λ
μ
SNR
λ
μ
SNR
ROF
220
0.005
18.5533
32
0.008
15.4055
LLT
150
0.005
16.2934
60
0.005
16.4336
Hybrid
220
0.005
20.2713
70
0.0025
16.6034The related images in Example4.12.
(a)
Original image
(b)
Mask image
(c)
Painting image
(d)
Noisy imageThe related restored images corresponding to the painting image in Example4.12.The related restored images corresponding to the noisy image in Example4.12.


## 4.2.3. Image Denoising for the Multiplicative Noise

Based on the hybrid model (1.7), the multiplicative noise model (1.5) can be naturally extended to the following minimization:
(4.19)minzα∫Ω‍(z+fe-z)dx+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.
It is obvious that the problem (4.19) can be approximated by(4.20a)un:=argminu12μ∥zn-u∥L22+α∫Ω‍(u+fe-u)dx,(4.20b)zn+1:=argminz12μ∥z-un∥L22+∫Ω‍|(1-g)∇z|dx+∫Ω‍|g∇2z|dx.For the first subproblem (4.20a), its solution un can be approximately obtained by using the Newton method to solve the following nonlinear equation:
(4.21)μ(u-zn)+α(1-fe-u)=0.Example 4.13.
In this example, we restore the multiplicative noisy image. The noisy Lena image shown in Figure7 is contaminated by the Gamma noise (L=33) with mean one, in which probability density function p(u) is given by
(4.22)p(u)={LLuL-1Γ(L)e-Lu,ifu>0,0,ifu≤0,
where L is an integer and Γ(·) is a Gamma function. Based on the iteration formula (4.20a)-(4.20b), we should notice that there are two interior iterations. For employing the Newton method to solve the problem (4.21), we set the stepsize t=1 and the Newton method will be stopped when ∥un+1-un∥ℓ2/∥un∥ℓ2≤1.0×10-5. For solving the second subproblem (4.20b), we set the fixed iterations with 40. In fact for using the PPM to solve the ROF model, the LLT model and the hybrid model, the solutions of these models will tend to the steady. The outer iteration will be stopped after 200 iterations. The related restored images are shown in Figure 8, and it is easy to find that the hybrid model gives a better restored image than the two other models.Figure 8
The related images and the local zooming images in Example4.13. SNRs: (a) 9.2017; (b) 15.6059; (c) 15.1754; (d) 16.0546.


## 5. Concluding Remarks

In this paper, based on the edge detector function, we proposed a hybrid model to overcome some drawbacks of the ROF model and the LLT model. Following the augmented Lagrangian method, we can employ the ADMM of multipliers to solve this hybrid model. In this paper, we mainly proposed the PPM to solve this model due to the fact that the PPM unnecessarily solves a PDE compared with the ADMM so that it is more effective than the ADMM. The convergence of the proposed method was also given. However, the convergence rate of the proposed method is onlyO(1/k), so our future work is to improve our method with convergence rate O(1/k2) based on the same strategies in [51, 52].



---

*Source: 376802-2012-11-29.xml*