# A New Fast Algorithm for Constrained Four-Directional Total Variation Image Denoising Problem

Authors: Fan Liao; Jean Louis Coatrieux; Jiasong Wu; Huazhong Shu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A new four-directional total variation (4-TV) model, applicable to isotropic and anisotropic TV functions, is proposed for image denoising. A dual based fast gradient projection algorithm for the constrained 4-TV image denoising problem is also reported which combines the well-known gradient projection and the fast gradient projection methods. Experimental results show that this model provides in most cases a better signal to noise ratio when compared to previous models like the reference TV, the total generalized variation, and the nonlocal total variation.


---

## Body

## 1. Introduction

Variational models have found a wide variety of applications in image processing and computer vision, in particular in restoration tasks such as denoising, deblurring, and blind deconvolution. One of the major concerns in machine vision remains to preserve important image features (edges, lines, and also textures) while removing noise. Total variation (TV) based image restoration models were first introduced by Rudin et al. (ROF) in their pioneering work [1] for edge preserving image denoising. It has been extended to many other problems and modified in a variety of ways to improve its performance. The unconstrained TV based denoising model has the following form:
(1)
min
⁡
x
x
-
b
2
+
λ
x
TV
,
where · is the Euclidian norm, x is the image to be recovered, b is the observed image, and ·
TV stands for the discrete TV norm. The positive parameter λ balances the measurements and the noise sensitivity.Many methods were proposed to solve (1). They include partial differential equation (PDE) [2, 3], semismooth Newton [4], multilevel optimization [5], primal-dual active set [6], second-order cone programming [7], high-order total variation minimization [8], split Bregman [9, 10], total generalized variation (TGV) [11], nonlocal total variation (NLTV) [12, 13], compressed sensing based [14], and dual methods [15–25].The methods of much interest to this paper are the dual approach proposed by Chambolle [20–22], the methods on the dual approach for the constrained denoising problem reported by Beck and Teboulle [23, 24], and the four-directional total variation (4-TV) proposed by Sakurai et al. [16]. Chambolle developed for the denoising problem a globally convergent first-order primal-dual algorithm, which is much faster than the conventional gradient descent algorithm. Beck and Teboulle used the gradient projection (GP) and the fast gradient projection (FGP) methods. Sakurai et al. [16] first proposed the four-directional total variation for anisotropic case, but a complete mathematical proof was not provided (it was admitted as such in their paper [16]). The FGP method relies on papers published by Nesterov [26, 27] where a fast first-order method is derived by coupling the gradient-based method with smoothing techniques. The rate of convergence of FGP is of the order O
(
1
/
k
2
) as opposed to the slower O
(
1
/
k
) rate of convergence of GP, where k is the number of iterations.In this paper, we provide the complete mathematical proof for both anisotropic and isotropic cases together with a proper definition of the 4-TV model. Moreover, we propose a fast constrained 4-TV algorithm for image denoising problem. To our knowledge, this is the first time that the double information in both time and space domains is jointly used.The paper is organized as follows. In Section2, the theoretical derivation of the 4-TV model is presented, and a new fast algorithm for this new model is described. Section 3 reports our experimental results including a comparison with the total generalized variation (TGV) and the nonlocal total variation (NLTV). A discussion is provided in Section 4. Finally, a short conclusion is given in Section 5.


## 2. Methods

### 2.1. The Discrete Four-Directional TV Model

We first consider the discrete ROF model which is a convex but nonsmooth minimization problem(2)
min
⁡
x
x
-
b
2
+
λ
TV
x
,
λ
>
0
,
where the discrete total variation in (1) is represented by TV
(
·
).Here, we only consider images defined on a rectangular domain, so we havex
∈
R
m
×
n and b
∈
R
m
×
n.The isotropic TV is defined as(3)
T
V
I
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
2
+
x
i
,
j
-
x
i
,
j
+
1
2
+
∑
i
=
1
m
-
1
x
i
,
n
-
x
i
+
1
,
n
+
∑
j
=
1
n
-
1
x
m
,
j
-
x
m
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
,
and the anisotropic TV as
(4)
T
V
A
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
+
x
i
,
j
-
x
i
,
j
+
1
+
∑
i
=
1
m
-
1
x
i
,
n
-
x
i
+
1
,
n
+
∑
j
=
1
n
-
1
x
m
,
j
-
x
m
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
.The discrete 4-TV model is an extended version of the conventional TV proposed by Sakurai et al. [16]. It relies on four different directional components (horizontal, vertical, and two different diagonal directional components) while the various conventional 2-directional TV just adopt 2 different components (horizontal and vertical directional components) as shown in Figure 1 (extracted from [16]).Total variation components [16].
(a)
Two-directional TV
(b)
Four-directional TVThe anisotropic 4-TV is defined by Sakurai et al. [16] as
(5)
T
V
4
A
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
+
x
i
,
j
-
x
i
,
j
+
1
hhhhhhhhh
+
x
i
,
j
-
x
i
+
1
,
j
+
1
hhhhhhhhh
+
x
i
+
1
,
j
-
x
i
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
h
h
x
∈
R
m
×
n
,
where, in the above formula, Sakurai et al. did not consider the boundary conditions. Moreover, the isotropic case was not taken into consideration in their paper.In order to get a more accurate 4-TV based model, we construct a new imagex
′.Letx
∈
R
m
×
n
(
1
≤
i
≤
m
,
1
≤
j
≤
n
) be the corrupted image and x
′
∈
R
(
m
+
2
)
×
(
n
+
2
)
(
0
≤
i
≤
m
+
1
,
0
≤
j
≤
n
+
1
) the new image. These images are shown as follows.(a)
The corrupted imagex is
(6)
x
1,1
…
x
1
,
n
⋮
⋱
⋮
x
m
,
1
…
x
m
,
n
.
(b)
The new imagex
′ is
(7)
x
0,0
′
x
0,1
′
⋯
x
0
,
n
′
x
0
,
n
+
1
′
x
1,0
′
x
1,1
′
⋯
x
1
,
n
′
x
1
,
n
+
1
′
⋮
⋮
⋱
⋮
⋮
x
m
,
0
′
x
m
,
1
′
⋯
x
m
,
n
′
x
m
,
n
+
1
′
x
m
+
1,0
′
x
m
+
1,1
′
⋯
x
m
+
1
,
n
′
x
m
+
1
,
n
+
1
′
.The images in (a) and (b) share the following relations:(8)
x
i
,
j
′
=
x
i
,
j
1
≤
i
≤
m
,
1
≤
j
≤
n
,
x
0
,
j
′
=
x
1
,
j
i
=
0
,
1
≤
j
≤
n
,
x
m
+
1
,
j
′
=
x
m
,
j
i
=
m
+
1
,
1
≤
j
≤
n
,
x
i
,
0
′
=
x
i
,
1
1
≤
i
≤
m
,
j
=
0
,
x
i
,
n
+
1
′
=
x
i
,
n
1
≤
i
≤
m
,
j
=
n
+
1
,
x
0,0
′
=
x
1,1
,
x
0
,
n
+
1
′
=
x
1
,
n
,
x
m
+
1,0
′
=
x
m
,
1
,
x
m
+
1
,
n
+
1
′
=
x
m
,
n
.Therefore, the anisotropic 4-TV can be defined as(9)
T
V
4
A
x
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
+
x
i
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
l
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
,
x
′
∈
R
(
m
+
2
)
×
(
n
+
2
)
and the isotropic 4-directional as
(10)
T
V
4
I
x
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
2
+
x
i
,
j
′
-
x
i
,
j
+
1
′
2
h
h
h
h
h
h
h
h
h
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
2
h
h
h
h
h
h
h
h
h
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
2
1
/
2
,
h
h
h
h
h
h
l
x
∈
R
m
×
n
,
x
′
∈
R
m
+
2
×
n
+
2
,
where the above boundary is expanded by constructing the new image x
′
∈
R
(
m
+
2
)
×
(
n
+
2
). Moreover, the domain of x
′
∈
R
(
m
+
2
)
×
(
n
+
2
) plays an important role in the 4-directional TV based dual method.Consequently, the discrete 4-directional TV model is defined by(11)
min
⁡
x
x
-
b
2
+
λ
T
V
4
x
.By adopting the above 4-directional TV model, each iteration step uses double information in the space domain. For this reason, the new model can be expected to have more general and effective properties than the standard one in image denoising problem. In addition, both anisotropic and isotropic TV cases can be dealt with.


### 2.2. Constrained 4-Directional Total Variation Based Denoising

#### 2.2.1. The Dual Approach

We consider the constrained 4-directional TV based denoising problem which corresponds to(12)
min
⁡
x
∈
C
x
-
b
2
+
λ
T
V
4
x
,
λ
>
0
,
where C is a closed convex subset of E
≡
R
m
×
n and the nonsmooth functional 4-directional TV is either anisotropic TV or isotropic TV. In order to avoid repetition, we will mainly consider the isotropic TV, and the results for the anisotropic TV will be briefly presented.The TV function is characterized by the nonsmoothness. The characteristic of the nonsmoothness is the key difficulty in problem equation (2). Chambolle [21] proposed a dual approach to surmount this shortcoming. Beck and Teboulle [23] and Sakurai et al. [16] followed the same approach and we also follow it for constructing our constrained 4-directional dual method.First, we definep
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
) as follows:
(13)
p
i
,
j
1
=
x
i
,
j
′
-
x
i
+
1
,
j
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
+
1
,
p
i
,
j
2
=
x
i
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
+
1
,
0
≤
j
≤
n
,
p
i
,
j
3
=
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
,
p
i
,
j
4
=
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
.Here, we do not assume reflexive boundary conditions owing to constructing the new image spacex
′.Now, letP be the set of matrix-group (
p
1
,
p
2
,
p
3
,
p
4
) that satisfies
(14)
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
≤
1
,
hhhhhhhhlhh
(
0
≤
i
≤
m
,
0
≤
j
≤
n
)
,
p
i
,
n
+
1
1
≤
1
,
0
≤
i
≤
m
,
p
m
+
1
,
j
2
≤
1
,
0
≤
j
≤
n
.Then, we also have the relation:(15)
x
1
2
+
x
2
2
+
x
3
2
+
x
4
2
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
p
1
x
1
+
p
2
x
2
+
p
3
x
3
+
p
4
x
4
:
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
h
h
h
h
h
h
h
h
h
h
h
l
h
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
.We introduce the linear operationΛ
∈
R
m
×
n, which is defined by
(16)
Λ
p
1
,
p
2
,
p
3
,
p
4
i
,
j
=
(
p
i
,
j
1
-
p
i
-
1
,
j
1
)
+
(
p
i
,
j
2
-
p
i
,
j
-
1
2
)
+
(
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
)
+
(
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
)
,
hhhhhhhhhhh
(
1
≤
i
≤
m
,
1
≤
j
≤
n
)
.So, we can write(17)
T
V
4
I
x
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
x
i
,
j
′
-
x
i
+
1
,
j
′
+
p
i
,
j
2
x
i
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
+
p
i
,
j
3
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
p
i
,
j
4
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
+
∑
i
=
1
m
p
i
,
0
1
2
+
p
i
,
n
+
1
1
2
+
∑
j
=
1
n
p
0
,
j
2
2
+
p
m
+
1
,
j
2
2
,
where 1
≤
i
≤
m
,
1
≤
j
≤
n. The term ∑
i
=
1
m
(
p
i
,
0
1
)
2
+
(
p
i
,
n
+
1
1
)
2
+
∑
j
=
1
n
(
p
0
,
j
2
)
2
+
(
p
m
+
1
,
j
2
)
2 is constant and can thus be omitted.So, (12) becomes
(18)
min
⁡
x
∈
C
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
x
-
b
2
+
λ
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
=
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
min
⁡
x
∈
C
x
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
+
b
2
.Because the above function is concave inp
1
,
p
2
,
p
3
,
p
4 and convex in x, we can exchange max and min [28].LetP
[
min
⁡
,
max
⁡
] be the orthogonal projection operator on the set [
min
⁡
,
max
⁡
]. So, P
[
min
⁡
,
max
⁡
] is given by
(19)
P
[
min
⁡
,
max
⁡
]
x
i
,
j
=
min
⁡
,
x
i
,
j
<
min
⁡
,
x
i
,
j
,
min
⁡
≤
x
i
,
j
≤
max
⁡
,
max
⁡
,
x
i
,
j
>
max
⁡
.The optimal solution of the constrained 4-directional TV based denoising model in (12) is
(20)
x
=
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
.By neglecting the constant termb
2 in (18), we obtain
(21)
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.Letp
1
,
p
2
,
p
3
,
p
4
∈
P be the optimal solution of
(22)
min
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
h
(
p
1
,
p
2
,
p
3
,
p
4
)
≡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.The operatorΛ
T
∈
R
m
×
n which is the adjoint to Λ is given by
(23)
Λ
T
x
=
p
1
,
p
2
,
p
3
,
p
4
.We consider the functions
∈
R
m
×
n defined by
(24)
s
x
=
x
-
P
[
min
⁡
,
max
⁡
]
x
2
.Equation (22) can be rewritten as
(25)
h
p
1
,
p
2
,
p
3
,
p
4
=
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.And the gradient ofs
(
·
) is given by
(26)
∇
s
x
=
2
x
-
P
[
min
⁡
,
max
⁡
]
x
.Therefore,(27)
∇
h
p
1
,
p
2
,
p
3
,
p
4
=
∇
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
l
h
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
λ
2
Λ
T
∇
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
-
λ
Λ
T
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
.From the above derivation, we see that the dual problem expressed by (22) is a convex minimization problem where the function is also continuously differentiable in a constraint set. Thus, the first-order gradient-based algorithms can be applied.The only difference between the isotopic TV and anisotropic TV cases is contained in the relation:(28)
x
=
max
⁡
p
{
p
x
:
p
≤
1
}
replacing (15).The gradient of (24) has been obtained. In order to solve the dual problem equation (22), we also need to calculate the Lipschitz constant of the gradient objective function [23]. Let L
(
h
) be the Lipschitz constant of the gradient objective function given by (22).The Euclidian norm of the matrix-pairsp
1
,
p
2
,
p
3
,
p
4, where p
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
), is
(29)
p
1
,
p
2
,
p
3
,
p
4
=
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
.For every two groups of matricesp
1
1
,
p
1
2
,
p
1
3
,
p
1
4, p
2
1
,
p
2
2
,
p
2
3
,
p
2
4, we have
(30)
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
=
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
Λ
T
·
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
h
l
h
h
h
h
h
h
h
h
l
h
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
T
·
Λ
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
Λ
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
≤
λ
2
2
Λ
T
·
Λ
·
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
.Now,(31)
Λ
T
x
2
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
2
+
x
i
,
j
′
-
x
i
,
j
+
1
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
+
1
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
x
i
,
j
′
2
=
16
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
2
,
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
-
p
i
-
1
,
j
1
2
+
p
i
,
j
2
-
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
+
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
2
+
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
2
+
p
i
-
1
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
h
+
p
i
,
j
3
2
+
p
i
-
1
,
j
-
1
3
2
h
h
h
h
h
h
h
h
h
+
p
i
-
1
,
j
4
2
+
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
p
1
,
p
2
,
p
3
,
p
4
2
=
16
∑
i
=
1
m
∑
j
=
1
n
p
1
,
p
2
,
p
3
,
p
4
2
.Thus, we haveΛ
(
p
1
,
p
2
,
p
3
,
p
4
)
≤
4
p
1
,
p
2
,
p
3
,
p
4 and Λ
T
(
x
)
≤
4
x
i
,
j. So, we can obtain that Λ
≤
4, Λ
T
≤
4, and
(32)
L
(
h
)
=
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
·
Λ
T
≤
8
λ
2
.The overall procedure to implement this constrained 4-directional gradient projection (4-GP) algorithm can be summarized as shown in Algorithm1.Algorithm 1:
4-GP  (
b
,
λ
,
N
).
Input:
b—Observed image.
λ—Regularization parameter.
N—Number of iterations.
Output:
x
*—Obtained optimal solution of (12) up to a tolerance.
Step  0.  Take (
p
0
1
,
p
0
2
,
p
0
3
,
p
0
4
)
=
(
0,0
,
0,0
).
Step  k.
(
k
=
1,2
,
…
,
N
)
Calculate
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
=
P
P
{
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
+
1
8
λ
Λ
T
(
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
)
)
}
Set
x
*
=
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
N
1
,
p
N
2
,
p
N
3
,
p
N
4
)
)In the constrained case, a groupp
1
,
p
2
,
p
3
,
p
4 is constrained by P
P
p
1
,
p
2
,
p
3
,
p
4
=
r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1, 0
≤
i
≤
m
,
0
≤
j
≤
n. The group r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1 is given by
(33)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.Note that the difference between the isotopic TV and anisotropic TV is as follows:(34)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
h
h
h
h
h
h
h
h
h
l
h
h
l
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.


### 2.3. The Fast Dual Approach

It has been shown from the above derivations that the dual problem equation (22) is a convex minimization problem where the function is also continuously differentiable and in a constraint set.The original fast gradient projection algorithm can be traced back to the gradient mapping approach proposed by Nesterov [26]. Since then, a number of new algorithms, inspired by Nesterov’s work, have been reported [17, 20, 23, 24, 27, 29].Here, we used the constrained 4-directional fast gradient projection (4-FGP) algorithm for the denoising problem. The 4-FGP algorithm has a convergence rate inO
(
1
/
k
2
) by utilizing double information (most recent two steps) in the time space better than the convergence rate O
(
1
/
k
) of the 4-GP algorithm. Following the FGP algorithm described in [23], the 4-FGP algorithm can be described in Algorithm 2.Algorithm 2:4-FGP (b, λ, N).
Input:
b—Observed image.
λ—Regularization parameter.
N—Number of iterations.
Output:
x
*—Obtained optimal solution of (12) up to a tolerance.
Step  0.  Take (
r
1
1
,
r
1
2
,
r
1
4
,
r
1
1
)
=
(
p
0
1
,
p
0
2
,
p
0
3
,
p
0
4
)
=
(
0,0
,
0,0
), t
1
=
1.
Step  k.
(
k
=
1,2
,
…
,
N
) Calculate
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
=
P
P
{
(
r
k
1
,
r
k
2
,
r
k
3
,
r
k
4
)
+
1
8
λ
Λ
T
(
P
[
min,max
]
(
b
-
λ
2
Λ
(
r
k
1
,
r
k
2
,
r
k
3
,
r
k
4
)
)
)
},
t
k
+
1
=
1
+
1
+
(
t
k
)
2
2,
(
r
k
+
1
1
,
r
k
+
1
2
,
r
k
+
1
3
,
r
k
+
1
4
)
=
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
+
(
t
k
-
1
t
k
+
1
)
(
p
k
1
-
p
k
-
1
1
,
p
k
2
-
p
k
-
1
2
,
p
k
3
-
p
k
-
1
3
,
p
k
4
-
p
k
-
1
4
)
Set  x
*
=
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
N
1
,
p
N
2
,
p
N
3
,
p
N
4
)
)
∑
12
12


## 2.1. The Discrete Four-Directional TV Model

We first consider the discrete ROF model which is a convex but nonsmooth minimization problem(2)
min
⁡
x
x
-
b
2
+
λ
TV
x
,
λ
>
0
,
where the discrete total variation in (1) is represented by TV
(
·
).Here, we only consider images defined on a rectangular domain, so we havex
∈
R
m
×
n and b
∈
R
m
×
n.The isotropic TV is defined as(3)
T
V
I
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
2
+
x
i
,
j
-
x
i
,
j
+
1
2
+
∑
i
=
1
m
-
1
x
i
,
n
-
x
i
+
1
,
n
+
∑
j
=
1
n
-
1
x
m
,
j
-
x
m
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
,
and the anisotropic TV as
(4)
T
V
A
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
+
x
i
,
j
-
x
i
,
j
+
1
+
∑
i
=
1
m
-
1
x
i
,
n
-
x
i
+
1
,
n
+
∑
j
=
1
n
-
1
x
m
,
j
-
x
m
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
.The discrete 4-TV model is an extended version of the conventional TV proposed by Sakurai et al. [16]. It relies on four different directional components (horizontal, vertical, and two different diagonal directional components) while the various conventional 2-directional TV just adopt 2 different components (horizontal and vertical directional components) as shown in Figure 1 (extracted from [16]).Total variation components [16].
(a)
Two-directional TV
(b)
Four-directional TVThe anisotropic 4-TV is defined by Sakurai et al. [16] as
(5)
T
V
4
A
x
=
∑
i
=
1
m
-
1
∑
j
=
1
n
-
1
x
i
,
j
-
x
i
+
1
,
j
+
x
i
,
j
-
x
i
,
j
+
1
hhhhhhhhh
+
x
i
,
j
-
x
i
+
1
,
j
+
1
hhhhhhhhh
+
x
i
+
1
,
j
-
x
i
,
j
+
1
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
h
h
x
∈
R
m
×
n
,
where, in the above formula, Sakurai et al. did not consider the boundary conditions. Moreover, the isotropic case was not taken into consideration in their paper.In order to get a more accurate 4-TV based model, we construct a new imagex
′.Letx
∈
R
m
×
n
(
1
≤
i
≤
m
,
1
≤
j
≤
n
) be the corrupted image and x
′
∈
R
(
m
+
2
)
×
(
n
+
2
)
(
0
≤
i
≤
m
+
1
,
0
≤
j
≤
n
+
1
) the new image. These images are shown as follows.(a)
The corrupted imagex is
(6)
x
1,1
…
x
1
,
n
⋮
⋱
⋮
x
m
,
1
…
x
m
,
n
.
(b)
The new imagex
′ is
(7)
x
0,0
′
x
0,1
′
⋯
x
0
,
n
′
x
0
,
n
+
1
′
x
1,0
′
x
1,1
′
⋯
x
1
,
n
′
x
1
,
n
+
1
′
⋮
⋮
⋱
⋮
⋮
x
m
,
0
′
x
m
,
1
′
⋯
x
m
,
n
′
x
m
,
n
+
1
′
x
m
+
1,0
′
x
m
+
1,1
′
⋯
x
m
+
1
,
n
′
x
m
+
1
,
n
+
1
′
.The images in (a) and (b) share the following relations:(8)
x
i
,
j
′
=
x
i
,
j
1
≤
i
≤
m
,
1
≤
j
≤
n
,
x
0
,
j
′
=
x
1
,
j
i
=
0
,
1
≤
j
≤
n
,
x
m
+
1
,
j
′
=
x
m
,
j
i
=
m
+
1
,
1
≤
j
≤
n
,
x
i
,
0
′
=
x
i
,
1
1
≤
i
≤
m
,
j
=
0
,
x
i
,
n
+
1
′
=
x
i
,
n
1
≤
i
≤
m
,
j
=
n
+
1
,
x
0,0
′
=
x
1,1
,
x
0
,
n
+
1
′
=
x
1
,
n
,
x
m
+
1,0
′
=
x
m
,
1
,
x
m
+
1
,
n
+
1
′
=
x
m
,
n
.Therefore, the anisotropic 4-TV can be defined as(9)
T
V
4
A
x
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
+
x
i
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
l
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
R
m
×
n
,
x
′
∈
R
(
m
+
2
)
×
(
n
+
2
)
and the isotropic 4-directional as
(10)
T
V
4
I
x
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
2
+
x
i
,
j
′
-
x
i
,
j
+
1
′
2
h
h
h
h
h
h
h
h
h
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
2
h
h
h
h
h
h
h
h
h
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
2
1
/
2
,
h
h
h
h
h
h
l
x
∈
R
m
×
n
,
x
′
∈
R
m
+
2
×
n
+
2
,
where the above boundary is expanded by constructing the new image x
′
∈
R
(
m
+
2
)
×
(
n
+
2
). Moreover, the domain of x
′
∈
R
(
m
+
2
)
×
(
n
+
2
) plays an important role in the 4-directional TV based dual method.Consequently, the discrete 4-directional TV model is defined by(11)
min
⁡
x
x
-
b
2
+
λ
T
V
4
x
.By adopting the above 4-directional TV model, each iteration step uses double information in the space domain. For this reason, the new model can be expected to have more general and effective properties than the standard one in image denoising problem. In addition, both anisotropic and isotropic TV cases can be dealt with.


## 2.2. Constrained 4-Directional Total Variation Based Denoising

### 2.2.1. The Dual Approach

We consider the constrained 4-directional TV based denoising problem which corresponds to(12)
min
⁡
x
∈
C
x
-
b
2
+
λ
T
V
4
x
,
λ
>
0
,
where C is a closed convex subset of E
≡
R
m
×
n and the nonsmooth functional 4-directional TV is either anisotropic TV or isotropic TV. In order to avoid repetition, we will mainly consider the isotropic TV, and the results for the anisotropic TV will be briefly presented.The TV function is characterized by the nonsmoothness. The characteristic of the nonsmoothness is the key difficulty in problem equation (2). Chambolle [21] proposed a dual approach to surmount this shortcoming. Beck and Teboulle [23] and Sakurai et al. [16] followed the same approach and we also follow it for constructing our constrained 4-directional dual method.First, we definep
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
) as follows:
(13)
p
i
,
j
1
=
x
i
,
j
′
-
x
i
+
1
,
j
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
+
1
,
p
i
,
j
2
=
x
i
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
+
1
,
0
≤
j
≤
n
,
p
i
,
j
3
=
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
,
p
i
,
j
4
=
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
.Here, we do not assume reflexive boundary conditions owing to constructing the new image spacex
′.Now, letP be the set of matrix-group (
p
1
,
p
2
,
p
3
,
p
4
) that satisfies
(14)
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
≤
1
,
hhhhhhhhlhh
(
0
≤
i
≤
m
,
0
≤
j
≤
n
)
,
p
i
,
n
+
1
1
≤
1
,
0
≤
i
≤
m
,
p
m
+
1
,
j
2
≤
1
,
0
≤
j
≤
n
.Then, we also have the relation:(15)
x
1
2
+
x
2
2
+
x
3
2
+
x
4
2
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
p
1
x
1
+
p
2
x
2
+
p
3
x
3
+
p
4
x
4
:
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
h
h
h
h
h
h
h
h
h
h
h
l
h
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
.We introduce the linear operationΛ
∈
R
m
×
n, which is defined by
(16)
Λ
p
1
,
p
2
,
p
3
,
p
4
i
,
j
=
(
p
i
,
j
1
-
p
i
-
1
,
j
1
)
+
(
p
i
,
j
2
-
p
i
,
j
-
1
2
)
+
(
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
)
+
(
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
)
,
hhhhhhhhhhh
(
1
≤
i
≤
m
,
1
≤
j
≤
n
)
.So, we can write(17)
T
V
4
I
x
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
x
i
,
j
′
-
x
i
+
1
,
j
′
+
p
i
,
j
2
x
i
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
+
p
i
,
j
3
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
p
i
,
j
4
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
+
∑
i
=
1
m
p
i
,
0
1
2
+
p
i
,
n
+
1
1
2
+
∑
j
=
1
n
p
0
,
j
2
2
+
p
m
+
1
,
j
2
2
,
where 1
≤
i
≤
m
,
1
≤
j
≤
n. The term ∑
i
=
1
m
(
p
i
,
0
1
)
2
+
(
p
i
,
n
+
1
1
)
2
+
∑
j
=
1
n
(
p
0
,
j
2
)
2
+
(
p
m
+
1
,
j
2
)
2 is constant and can thus be omitted.So, (12) becomes
(18)
min
⁡
x
∈
C
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
x
-
b
2
+
λ
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
=
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
min
⁡
x
∈
C
x
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
+
b
2
.Because the above function is concave inp
1
,
p
2
,
p
3
,
p
4 and convex in x, we can exchange max and min [28].LetP
[
min
⁡
,
max
⁡
] be the orthogonal projection operator on the set [
min
⁡
,
max
⁡
]. So, P
[
min
⁡
,
max
⁡
] is given by
(19)
P
[
min
⁡
,
max
⁡
]
x
i
,
j
=
min
⁡
,
x
i
,
j
<
min
⁡
,
x
i
,
j
,
min
⁡
≤
x
i
,
j
≤
max
⁡
,
max
⁡
,
x
i
,
j
>
max
⁡
.The optimal solution of the constrained 4-directional TV based denoising model in (12) is
(20)
x
=
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
.By neglecting the constant termb
2 in (18), we obtain
(21)
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.Letp
1
,
p
2
,
p
3
,
p
4
∈
P be the optimal solution of
(22)
min
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
h
(
p
1
,
p
2
,
p
3
,
p
4
)
≡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.The operatorΛ
T
∈
R
m
×
n which is the adjoint to Λ is given by
(23)
Λ
T
x
=
p
1
,
p
2
,
p
3
,
p
4
.We consider the functions
∈
R
m
×
n defined by
(24)
s
x
=
x
-
P
[
min
⁡
,
max
⁡
]
x
2
.Equation (22) can be rewritten as
(25)
h
p
1
,
p
2
,
p
3
,
p
4
=
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.And the gradient ofs
(
·
) is given by
(26)
∇
s
x
=
2
x
-
P
[
min
⁡
,
max
⁡
]
x
.Therefore,(27)
∇
h
p
1
,
p
2
,
p
3
,
p
4
=
∇
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
l
h
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
λ
2
Λ
T
∇
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
-
λ
Λ
T
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
.From the above derivation, we see that the dual problem expressed by (22) is a convex minimization problem where the function is also continuously differentiable in a constraint set. Thus, the first-order gradient-based algorithms can be applied.The only difference between the isotopic TV and anisotropic TV cases is contained in the relation:(28)
x
=
max
⁡
p
{
p
x
:
p
≤
1
}
replacing (15).The gradient of (24) has been obtained. In order to solve the dual problem equation (22), we also need to calculate the Lipschitz constant of the gradient objective function [23]. Let L
(
h
) be the Lipschitz constant of the gradient objective function given by (22).The Euclidian norm of the matrix-pairsp
1
,
p
2
,
p
3
,
p
4, where p
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
), is
(29)
p
1
,
p
2
,
p
3
,
p
4
=
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
.For every two groups of matricesp
1
1
,
p
1
2
,
p
1
3
,
p
1
4, p
2
1
,
p
2
2
,
p
2
3
,
p
2
4, we have
(30)
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
=
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
Λ
T
·
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
h
l
h
h
h
h
h
h
h
h
l
h
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
T
·
Λ
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
Λ
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
≤
λ
2
2
Λ
T
·
Λ
·
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
.Now,(31)
Λ
T
x
2
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
2
+
x
i
,
j
′
-
x
i
,
j
+
1
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
+
1
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
x
i
,
j
′
2
=
16
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
2
,
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
-
p
i
-
1
,
j
1
2
+
p
i
,
j
2
-
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
+
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
2
+
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
2
+
p
i
-
1
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
h
+
p
i
,
j
3
2
+
p
i
-
1
,
j
-
1
3
2
h
h
h
h
h
h
h
h
h
+
p
i
-
1
,
j
4
2
+
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
p
1
,
p
2
,
p
3
,
p
4
2
=
16
∑
i
=
1
m
∑
j
=
1
n
p
1
,
p
2
,
p
3
,
p
4
2
.Thus, we haveΛ
(
p
1
,
p
2
,
p
3
,
p
4
)
≤
4
p
1
,
p
2
,
p
3
,
p
4 and Λ
T
(
x
)
≤
4
x
i
,
j. So, we can obtain that Λ
≤
4, Λ
T
≤
4, and
(32)
L
(
h
)
=
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
·
Λ
T
≤
8
λ
2
.The overall procedure to implement this constrained 4-directional gradient projection (4-GP) algorithm can be summarized as shown in Algorithm1.Algorithm 1:
4-GP  (
b
,
λ
,
N
).
Input:
b—Observed image.
λ—Regularization parameter.
N—Number of iterations.
Output:
x
*—Obtained optimal solution of (12) up to a tolerance.
Step  0.  Take (
p
0
1
,
p
0
2
,
p
0
3
,
p
0
4
)
=
(
0,0
,
0,0
).
Step  k.
(
k
=
1,2
,
…
,
N
)
Calculate
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
=
P
P
{
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
+
1
8
λ
Λ
T
(
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
)
)
}
Set
x
*
=
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
N
1
,
p
N
2
,
p
N
3
,
p
N
4
)
)In the constrained case, a groupp
1
,
p
2
,
p
3
,
p
4 is constrained by P
P
p
1
,
p
2
,
p
3
,
p
4
=
r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1, 0
≤
i
≤
m
,
0
≤
j
≤
n. The group r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1 is given by
(33)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.Note that the difference between the isotopic TV and anisotropic TV is as follows:(34)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
h
h
h
h
h
h
h
h
h
l
h
h
l
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.


## 2.2.1. The Dual Approach

We consider the constrained 4-directional TV based denoising problem which corresponds to(12)
min
⁡
x
∈
C
x
-
b
2
+
λ
T
V
4
x
,
λ
>
0
,
where C is a closed convex subset of E
≡
R
m
×
n and the nonsmooth functional 4-directional TV is either anisotropic TV or isotropic TV. In order to avoid repetition, we will mainly consider the isotropic TV, and the results for the anisotropic TV will be briefly presented.The TV function is characterized by the nonsmoothness. The characteristic of the nonsmoothness is the key difficulty in problem equation (2). Chambolle [21] proposed a dual approach to surmount this shortcoming. Beck and Teboulle [23] and Sakurai et al. [16] followed the same approach and we also follow it for constructing our constrained 4-directional dual method.First, we definep
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
) as follows:
(13)
p
i
,
j
1
=
x
i
,
j
′
-
x
i
+
1
,
j
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
+
1
,
p
i
,
j
2
=
x
i
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
+
1
,
0
≤
j
≤
n
,
p
i
,
j
3
=
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
,
p
i
,
j
4
=
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
,
0
≤
i
≤
m
,
0
≤
j
≤
n
.Here, we do not assume reflexive boundary conditions owing to constructing the new image spacex
′.Now, letP be the set of matrix-group (
p
1
,
p
2
,
p
3
,
p
4
) that satisfies
(14)
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
≤
1
,
hhhhhhhhlhh
(
0
≤
i
≤
m
,
0
≤
j
≤
n
)
,
p
i
,
n
+
1
1
≤
1
,
0
≤
i
≤
m
,
p
m
+
1
,
j
2
≤
1
,
0
≤
j
≤
n
.Then, we also have the relation:(15)
x
1
2
+
x
2
2
+
x
3
2
+
x
4
2
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
p
1
x
1
+
p
2
x
2
+
p
3
x
3
+
p
4
x
4
:
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
h
h
h
h
h
h
h
h
h
h
h
l
h
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
≤
1
.We introduce the linear operationΛ
∈
R
m
×
n, which is defined by
(16)
Λ
p
1
,
p
2
,
p
3
,
p
4
i
,
j
=
(
p
i
,
j
1
-
p
i
-
1
,
j
1
)
+
(
p
i
,
j
2
-
p
i
,
j
-
1
2
)
+
(
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
)
+
(
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
)
,
hhhhhhhhhhh
(
1
≤
i
≤
m
,
1
≤
j
≤
n
)
.So, we can write(17)
T
V
4
I
x
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
x
i
,
j
′
-
x
i
+
1
,
j
′
+
p
i
,
j
2
x
i
,
j
′
-
x
i
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
+
p
i
,
j
3
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
p
i
,
j
4
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
=
max
⁡
p
1
,
p
2
,
p
3
,
p
4
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
+
∑
i
=
1
m
p
i
,
0
1
2
+
p
i
,
n
+
1
1
2
+
∑
j
=
1
n
p
0
,
j
2
2
+
p
m
+
1
,
j
2
2
,
where 1
≤
i
≤
m
,
1
≤
j
≤
n. The term ∑
i
=
1
m
(
p
i
,
0
1
)
2
+
(
p
i
,
n
+
1
1
)
2
+
∑
j
=
1
n
(
p
0
,
j
2
)
2
+
(
p
m
+
1
,
j
2
)
2 is constant and can thus be omitted.So, (12) becomes
(18)
min
⁡
x
∈
C
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
x
-
b
2
+
λ
Tr
⁡
Λ
p
1
,
p
2
,
p
3
,
p
4
T
x
=
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
min
⁡
x
∈
C
x
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
+
b
2
.Because the above function is concave inp
1
,
p
2
,
p
3
,
p
4 and convex in x, we can exchange max and min [28].LetP
[
min
⁡
,
max
⁡
] be the orthogonal projection operator on the set [
min
⁡
,
max
⁡
]. So, P
[
min
⁡
,
max
⁡
] is given by
(19)
P
[
min
⁡
,
max
⁡
]
x
i
,
j
=
min
⁡
,
x
i
,
j
<
min
⁡
,
x
i
,
j
,
min
⁡
≤
x
i
,
j
≤
max
⁡
,
max
⁡
,
x
i
,
j
>
max
⁡
.The optimal solution of the constrained 4-directional TV based denoising model in (12) is
(20)
x
=
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
.By neglecting the constant termb
2 in (18), we obtain
(21)
max
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
h
h
h
h
h
h
h
h
h
h
h
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.Letp
1
,
p
2
,
p
3
,
p
4
∈
P be the optimal solution of
(22)
min
⁡
(
p
1
,
p
2
,
p
3
,
p
4
)
∈
P
h
(
p
1
,
p
2
,
p
3
,
p
4
)
≡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
,
p
2
,
p
3
,
p
4
)
-
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.The operatorΛ
T
∈
R
m
×
n which is the adjoint to Λ is given by
(23)
Λ
T
x
=
p
1
,
p
2
,
p
3
,
p
4
.We consider the functions
∈
R
m
×
n defined by
(24)
s
x
=
x
-
P
[
min
⁡
,
max
⁡
]
x
2
.Equation (22) can be rewritten as
(25)
h
p
1
,
p
2
,
p
3
,
p
4
=
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
.And the gradient ofs
(
·
) is given by
(26)
∇
s
x
=
2
x
-
P
[
min
⁡
,
max
⁡
]
x
.Therefore,(27)
∇
h
p
1
,
p
2
,
p
3
,
p
4
=
∇
-
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
h
h
h
l
h
+
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
λ
2
Λ
T
∇
s
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
-
λ
Λ
T
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
=
-
λ
Λ
T
P
min
⁡
,
max
⁡
b
-
λ
2
Λ
p
1
,
p
2
,
p
3
,
p
4
.From the above derivation, we see that the dual problem expressed by (22) is a convex minimization problem where the function is also continuously differentiable in a constraint set. Thus, the first-order gradient-based algorithms can be applied.The only difference between the isotopic TV and anisotropic TV cases is contained in the relation:(28)
x
=
max
⁡
p
{
p
x
:
p
≤
1
}
replacing (15).The gradient of (24) has been obtained. In order to solve the dual problem equation (22), we also need to calculate the Lipschitz constant of the gradient objective function [23]. Let L
(
h
) be the Lipschitz constant of the gradient objective function given by (22).The Euclidian norm of the matrix-pairsp
1
,
p
2
,
p
3
,
p
4, where p
1
∈
R
(
m
+
1
)
×
(
n
+
2
), p
2
∈
R
(
m
+
2
)
×
(
n
+
1
), p
3
∈
R
(
m
+
1
)
×
(
n
+
1
), p
4
∈
R
(
m
+
1
)
×
(
n
+
1
), is
(29)
p
1
,
p
2
,
p
3
,
p
4
=
p
1
2
+
p
2
2
+
p
3
2
+
p
4
2
.For every two groups of matricesp
1
1
,
p
1
2
,
p
1
3
,
p
1
4, p
2
1
,
p
2
2
,
p
2
3
,
p
2
4, we have
(30)
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
=
λ
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
Λ
T
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
Λ
T
·
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
h
l
h
h
h
h
h
h
h
h
l
h
-
P
[
min
⁡
,
max
⁡
]
b
-
λ
2
Λ
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
T
·
Λ
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
Λ
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
≤
λ
2
2
Λ
T
·
Λ
·
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
-
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
.Now,(31)
Λ
T
x
2
=
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
-
x
i
+
1
,
j
′
2
+
x
i
,
j
′
-
x
i
,
j
+
1
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
′
-
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
-
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
,
j
+
1
′
2
+
x
i
,
j
′
2
h
h
h
h
h
h
h
h
+
x
i
+
1
,
j
+
1
′
2
+
x
i
+
1
,
j
′
2
+
x
i
,
j
+
1
′
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
x
i
,
j
′
2
=
16
∑
i
=
1
m
∑
j
=
1
n
x
i
,
j
2
,
Λ
p
1
,
p
2
,
p
3
,
p
4
2
=
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
-
p
i
-
1
,
j
1
2
+
p
i
,
j
2
-
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
+
p
i
,
j
3
-
p
i
-
1
,
j
-
1
3
2
+
p
i
-
1
,
j
4
-
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
p
i
,
j
1
2
+
p
i
-
1
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
-
1
2
2
h
h
h
h
h
h
h
h
h
+
p
i
,
j
3
2
+
p
i
-
1
,
j
-
1
3
2
h
h
h
h
h
h
h
h
h
+
p
i
-
1
,
j
4
2
+
p
i
,
j
-
1
4
2
≤
2
∑
i
=
1
m
∑
j
=
1
n
8
p
1
,
p
2
,
p
3
,
p
4
2
=
16
∑
i
=
1
m
∑
j
=
1
n
p
1
,
p
2
,
p
3
,
p
4
2
.Thus, we haveΛ
(
p
1
,
p
2
,
p
3
,
p
4
)
≤
4
p
1
,
p
2
,
p
3
,
p
4 and Λ
T
(
x
)
≤
4
x
i
,
j. So, we can obtain that Λ
≤
4, Λ
T
≤
4, and
(32)
L
(
h
)
=
∇
h
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
∇
h
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
(
p
1
1
,
p
1
2
,
p
1
3
,
p
1
4
)
-
(
p
2
1
,
p
2
2
,
p
2
3
,
p
2
4
)
≤
λ
2
2
Λ
·
Λ
T
≤
8
λ
2
.The overall procedure to implement this constrained 4-directional gradient projection (4-GP) algorithm can be summarized as shown in Algorithm1.Algorithm 1:
4-GP  (
b
,
λ
,
N
).
Input:
b—Observed image.
λ—Regularization parameter.
N—Number of iterations.
Output:
x
*—Obtained optimal solution of (12) up to a tolerance.
Step  0.  Take (
p
0
1
,
p
0
2
,
p
0
3
,
p
0
4
)
=
(
0,0
,
0,0
).
Step  k.
(
k
=
1,2
,
…
,
N
)
Calculate
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
=
P
P
{
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
+
1
8
λ
Λ
T
(
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
k
-
1
1
,
p
k
-
1
2
,
p
k
-
1
3
,
p
k
-
1
4
)
)
)
}
Set
x
*
=
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
N
1
,
p
N
2
,
p
N
3
,
p
N
4
)
)In the constrained case, a groupp
1
,
p
2
,
p
3
,
p
4 is constrained by P
P
p
1
,
p
2
,
p
3
,
p
4
=
r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1, 0
≤
i
≤
m
,
0
≤
j
≤
n. The group r
1
,
r
2
,
r
3
,
r
4
∈
R
m
+
1
×
n
+
1 is given by
(33)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
2
+
p
i
,
j
2
2
+
p
i
,
j
3
2
+
p
i
,
j
4
2
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.Note that the difference between the isotopic TV and anisotropic TV is as follows:(34)
r
i
,
j
1
=
p
i
,
j
1
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
2
=
p
i
,
j
2
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
3
=
p
i
,
j
3
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
r
i
,
j
4
=
p
i
,
j
4
max
⁡
1
,
p
i
,
j
1
+
p
i
,
j
2
+
p
i
,
j
3
+
p
i
,
j
4
,
h
h
h
h
h
h
h
h
h
l
h
h
l
h
h
0
≤
i
≤
m
,
0
≤
j
≤
n
.


## 2.3. The Fast Dual Approach

It has been shown from the above derivations that the dual problem equation (22) is a convex minimization problem where the function is also continuously differentiable and in a constraint set.The original fast gradient projection algorithm can be traced back to the gradient mapping approach proposed by Nesterov [26]. Since then, a number of new algorithms, inspired by Nesterov’s work, have been reported [17, 20, 23, 24, 27, 29].Here, we used the constrained 4-directional fast gradient projection (4-FGP) algorithm for the denoising problem. The 4-FGP algorithm has a convergence rate inO
(
1
/
k
2
) by utilizing double information (most recent two steps) in the time space better than the convergence rate O
(
1
/
k
) of the 4-GP algorithm. Following the FGP algorithm described in [23], the 4-FGP algorithm can be described in Algorithm 2.Algorithm 2:4-FGP (b, λ, N).
Input:
b—Observed image.
λ—Regularization parameter.
N—Number of iterations.
Output:
x
*—Obtained optimal solution of (12) up to a tolerance.
Step  0.  Take (
r
1
1
,
r
1
2
,
r
1
4
,
r
1
1
)
=
(
p
0
1
,
p
0
2
,
p
0
3
,
p
0
4
)
=
(
0,0
,
0,0
), t
1
=
1.
Step  k.
(
k
=
1,2
,
…
,
N
) Calculate
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
=
P
P
{
(
r
k
1
,
r
k
2
,
r
k
3
,
r
k
4
)
+
1
8
λ
Λ
T
(
P
[
min,max
]
(
b
-
λ
2
Λ
(
r
k
1
,
r
k
2
,
r
k
3
,
r
k
4
)
)
)
},
t
k
+
1
=
1
+
1
+
(
t
k
)
2
2,
(
r
k
+
1
1
,
r
k
+
1
2
,
r
k
+
1
3
,
r
k
+
1
4
)
=
(
p
k
1
,
p
k
2
,
p
k
3
,
p
k
4
)
+
(
t
k
-
1
t
k
+
1
)
(
p
k
1
-
p
k
-
1
1
,
p
k
2
-
p
k
-
1
2
,
p
k
3
-
p
k
-
1
3
,
p
k
4
-
p
k
-
1
4
)
Set  x
*
=
P
[
min,max
]
(
b
-
λ
2
Λ
(
p
N
1
,
p
N
2
,
p
N
3
,
p
N
4
)
)
∑
12
12


## 3. Experimental Results

These experiments were conducted on images widely used in the computer vision literature. We selected two samples among this trial set, the “Cameraman” and the “Moon” pictures, to illustrate the effectiveness of the proposed method. These two images, by their different contents, are representative of the large spectrum of data sets that can be considered. A comparison of our methods (4-GP and 4-FGP) was performed with the GP, FGP [23], TGV [11], and NLTV [12] methods. The peak signal to noise ratio (PSNR), the convergence rate, and the robustness to noise were used for the evaluation of the denoised image quality. All algorithms have been implemented on a PC Intel Duo Core CPU E8400 3 GHz, RAM 8 GB with MATLAB R2011b.The300
×
300 gray-scale test image “Cameraman” and the 348
×
300 gray-scale test image “Moon” (Figures 2 and 3) were scaled in intensity to [
0
,
1
]. A normally distributed zero-mean Gaussian noise was then added, with standard deviations equal to 0.11 and 0.07 for the “Cameraman” image and the “Moon” image, respectively.Comparison of the different methods on the “Cameraman” image. (a) Original image; (b) noisy image, PSNR = 19.60 dB; (c) denoised image with GP, PSNR = 26.19 dB; (d) denoised image with FGP, PSNR = 26.19 dB; (e) denoised image with TGV, PSNR = 26.38 dB; (f) denoised image with NLTV, PSNR = 27.13 dB; (g) denoised image with 4-GP, PSNR = 27.41 dB; and (h) denoised image with 4-FGP, PSNR = 27.41 dB.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Results on the “Moon” image. (a) Original image; (b) noisy image, PSNR = 24.82 dB; (c) denoised image with GP, PSNR = 29.35 dB; (d) denoised image with FGP, PSNR = 29.35 dB; (e) denoised image with TGV, PSNR = 30.08 dB; (f) denoised image with NLTV, PSNR = 30.23 dB; (g) denoised image with 4-GP, PSNR = 30.35 dB; (h) denoised image with 4-FGP, PSNR = 30.36 dB.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)The parameters were set toλ
=
0.1 for the “Cameraman” image and to λ
=
0.05 for the “Moon” image in all experiments. The tolerance value for the convergence test was set to 0.0001 dB.For the TGV method [11], the parameters were set to α
0
=
1, α
1
=
2, σ
=
0.1, and τ
=
0.2 throughout this paper. For the NLTV method [12], we set the patch size as 5
×
5, the number of neighbors as k
=
20, and the searching window as 11
×
11.The MATLAB codes of the TGV method [11] and the NLTV method [12] are, respectively, available at http://www.imt.tugraz.at/research/mr-image-reconstruction and http://math.sjtu.edu.cn/faculty/xqzhang/html/code.html.The PSNR values obtained in the above cases for the different methods are indicated in Figures2 and 3 captions. They show that the convergence values when using 4-GP or 4-FGP are almost identical. They also show that the 4-GP and 4-FGP methods lead to a PSNR gain, ranging, for instance, from 0.3 dB to 1 dB for the “Cameraman” image.From Figures2 and 3, we can see that all methods have their own advantages and drawbacks. The GP and FGP methods remove the noise but still produce a staircasing effect in the flat and smooth regions. The TGV method clears up this effect while preserving the edges when the prior information is very close to the original image, but it leads to false image features when the prior information has been corrupted by a high level of noise. The NLTV method reduces the staircasing effect but some details are lost. The 4-GP and 4-FGP methods well preserve the edges and capture more details because the four different directional components are taken into consideration.As it was expected, the convergence is much faster when using 4-FGP instead of 4-GP. A detailed analysis of the convergence process makes it clear that the number of iterations is image dependent and much higher for GP than for FGP; for the “Cameraman” image, this number is equal to 124 for 4-GP and to 41 for 4-FGP; for the “Moon” image, they are, respectively, equal to 103 and 51. The time computation varies accordingly; it goes for the “Cameraman” image from about 55 seconds for 4-GP to approximately 18 seconds for 4-FGP and for the “Moon” image from about 53 seconds for 4-GP to approximately 27 seconds for 4-FGP.From Figure4, we can find that the convergence values of 4-GP and 4-FGP are almost the same. They are similar to those obtained by GP and FGP. So, we will just consider FGP and 4-FGP in the next experiments in order to avoid repetition.Comparison of the convergence behavior for the different methods with respect to PSNR. (a) The “Cameraman” image (standard deviation of zero-mean Gaussian white noise: 0.11) and (b) the “Moon” image (standard deviation of zero-mean Gaussian white noise: 0.07).
(a)
(b)In the examples described above, the values of the noise andλ were a priori set. Let us now examine the sensitivity of these parameters to explain the rationale behind our choices. We varied here the values of λ from 0 to 0.1 and from 0 to 0.05 for the “Cameraman” image and the “Moon” image, respectively.The results are provided in Figure5 and show that the 4-FGP method performs always better than the GP and NLTV methods in terms of PSNR for the two images, the only exception being for the TGV method when the value of λ is small.Evolution of PSNR withλ varying for the different methods. (a) The “Cameraman” image (the value of noise is 0.11) and (b) the “Moon” image (the value of noise is 0.07).
(a)
(b)By taking into account the previous results, the noise effect was analyzed. The selected parameters wereλ
=
0.1 for the “Cameraman” image and λ
=
0.05 for the “Moon” image. The noise level was varied from 0 to 0.3 in the first case and from 0 to 0.15 in the second case. Figure 6 depicts the evolution of PSNR.Evolution of PSNR when varying the noise level for the different methods. (a) The “Cameraman” image (λ
=
0.1) and (b) the “Moon” image (λ
=
0.05).
(a)
(b)The performances of the FGP, TGV, and NLTV methods are inferior to the performance obtained with 4-FGP method when the value of noise is large and superior to the performance of the 4-FGP method when the value of noise is rather low.


## 4. Discussion

These methods were compared on three other images and additional comments are provided in this section. The objective was to see if they were stable enough to be generalized to any type of images. Two images, “Lena” and “Woman,” were used by Chambolle in [21] and the third one, the “Louvre” image, was used in [20]. The experimental conditions were a Gaussian noise level 0.1 and λ
=
0.1 for all three images. Table 1 presents the PSNR obtained by the different methods. These results confirm the interest of 4-TV even if the benefits remain in a modest range.Table 1
PSNR (dB) for GP, FGP, TGV, NLTV, 4-GP, and 4-FGP on different images  (“Lena,” “Woman,” and “Louvre” images, standard deviation of zero-mean Gaussian white noise 0.1 andλ
=
0.1).
Image
Noisy image
GP
FGP
TGV
NLTV
4-GP
4-FGP
Lena
20.43
26.22
26.23
26.55
26.60
26.67
26.67
Woman
20.50
27.48
27.48
27.70
27.78
28.39
28.38
Louvre
20.35
27.02
27.02
26.91
27.45
27.50
27.49


## 5. Conclusion

In this paper, 4-GP and 4-FGP methods were proposed for image denoising. We added the diagonal components to the conventional TV model and we provided the complete mathematical proof of the relevance of the new model. Moreover, the 4-FGP algorithm makes for the first time use of the double information in both time and space domains. Experimental results show that the 4-GP and 4-FGP methods lead to better denoising results in most cases. In the future work, we will address the weighted 4-GP and 4-FGP frames.



---

*Source: 815132-2015-02-23.xml*