# A Fast Alternating Minimization Algorithm for Nonlocal Vectorial Total Variational Multichannel Image Denoising

Authors: Rubing Xi; Zhengming Wang; Xia Zhao; Meihua Xie
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The variational models with nonlocal regularization offer superior image restoration quality over traditional method. But the processing speed remains a bottleneck due to the calculation quantity brought by the recent iterative algorithms. In this paper, a fast algorithm is proposed to restore the multichannel image in the presence of additive Gaussian noise by minimizing an energy function consisting of anl
2-norm fidelity term and a nonlocal vectorial total variational regularization term. This algorithm is based on the variable splitting and penalty techniques in optimization. Following our previous work on the proof of the existence and the uniqueness of the solution of the model, we establish and prove the convergence properties of this algorithm, which are the finite convergence for some variables and the q-linear convergence for the rest. Experiments show that this model has a fabulous texture-preserving property in restoring color images. Both the theoretical derivation of the computation complexity analysis and the experimental results show that the proposed algorithm performs favorably in comparison to the widely used fixed point algorithm.


---

## Body

## 1. Introduction

The restoration of multichannel image is one of the most important inverse problems in image processing. Among all the methods for restoration, the vectorial total variational (VTV) regularization is most popular for preserving discontinuities in restored images. VTV has been introduced by [1] into multichannel image restoration as a vectorial extension of total variation (TV) in gray scale image processing.The famous total variational (TV) methods have drawbacks due to its local property while calculating the gradient of the image. Thus the images that contain many small structures or textures cannot get satisfying restored results. As a result, TV method based on nonlocal (NL) operator is generally used in gray scale image processing because of its structure preserving property. This method is generalized from the Yaroslavsky denoising and the block method. The main idea is using the similar blocks in the whole image to estimate current pixel value. It was introduced by Efros and Leung [2] while studying texture synthesis. Then Buades et al. [3] generalized nonlocal method and proposed the well-known neighbor denoising filter, namely, nonlocal means (NLM). After that, Gilboa and Osher [4] united nonlocal method into the regularization framework by defining the nonlocal formed TV model including NL-TV-l
1 model, NL-TV-l
2 model, and NL-TV-G model and proposed the algorithms separately. NL-TV model can both preserve the edges and the other small structures in images and has excellent application effect.Our previews works extends NL-TV regularization to NL-VTV regularization for multichannel image restoration (see [5]). The corresponding algorithms are based on the variation calculus by solving the Euler-Lagrangian equations by fixed point iterations. This kind of methods is tested effectively but has a low convergence rate.Recently, the optimization algorithms based on variable-splitting and penalty techniques have attracted much attention in image and signal processing, such as the dual methods, the Bregman methods, and the alternating direction method of multipliers (ADMM) methods (see [6–8]). Correspondingly, these methods have been extended for multichannel images. Bresson and Chan [6] has extended Chambolle’s dual algorithm for TV problem to vectorial images by defining the standard dual definition of VTV. Wang et al. [9, 10] has applied the alternating minimization algorithm to VTV model for multichannel image deblurring.In this paper, we propose a fast alternating algorithm for NL-VTV model to restore color images from noisy observations. The following sections are organized as follows. Section2 is the brief introduction of the NL-VTV model and classical fixed point algorithm. In Section 3, we propose the alternating algorithm for NL-VTV model. In Section 4, we do the convergence analysis of this algorithm. In Section 5, we applied this algorithm to color image denoising to present the texture-preserving effect of the NL-VTV model and the high convergent rate of the proposed algorithm.


## 2. Nonlocal Vectorial Total Variation Model

In this section, we will briefly introduce the nonlocal vectorial total variational model for multichannel image processing. The usually used iterative algorithm will be introduced as well, which will be compared to our algorithm in Section5.


### 2.1. Vectorial Total Variation Model

Consider a vector valued functionu, such as the multichannel images, defined on a bounded open domain Ω
⊂
R
N:
(1)
u
:
Ω
⟶
R
M
x
⟼
u
(
x
)
=
(
u
1
(
x
)
,
u
2
(
x
)
,
…
,
u
M
(
x
)
)
,
where u
m
:
Ω
→
R
(
1
≤
m
≤
M
) denotes the mth channel of u.The regularized VTV model for multichannel image reconstruction could be written as(2)
inf
⁡
u
∈
R
M
E
(
u
)
≔
∫
Ω
|
∇
u
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
where u
0 is the observed noisy multichannel image and λ
>
0 is the regularized parameter. The vectorial gradient of u is defined as
(3)
∇
u
=
(
∇
u
1
,
…
,
∇
u
M
)
=
[
∂
u
1
∂
x
∂
u
2
∂
x
⋯
∂
u
M
∂
x
∂
u
1
∂
y
∂
u
2
∂
y
⋯
∂
u
M
∂
y
]
∈
R
2
×
M
,
|
∇
u
|
=
∥
∇
u
∥
F
=
∑
m
=
1
M
(
∂
u
m
∂
x
)
2
+
∑
m
=
1
M
(
∂
u
m
∂
y
)
2
=
∑
m
=
1
M
|
∇
u
m
|
2
.
We refer to [6] for more information about the well-posedness analysis of the VTV model.


### 2.2. Nonlocal Functional

The neighborhood filter NLM proposed by Buades et al. developed to a famous image denoising method [3]. It uses the pixels similar to the current pixel to estimate the value of the current pixel. Let the reference image be ν. The similarity of the two pixels in the image x, y could be measured as the following weight function w
ν
(
x
,
y
):
(4)
w
ν
(
x
,
y
)
=
exp
⁡
(
-
(
G
a
*
|
ν
(
x
+
·
)
-
ν
(
y
+
·
)
|
2
)
(
0
)
h
2
)
,
where G
a is the Gaussian kernel with the standard deviation of a:
(5)
(
G
a
*
|
ν
(
x
+
·
)
-
ν
(
y
+
·
)
|
2
)
(
0
)
=
∫
R
2
G
a
(
t
)
|
ν
(
x
+
t
)
-
ν
(
y
+
t
)
|
2
d
t
,
where h is the filter parameter. Generally, h corresponds to the noise level, which is usually set as the standard deviation of the noise. It is shown that the weight function above could be calculated by the image blocks centered around x, y. The value of the weight function reflects the similarity of the two blocks. The estimation of the value at pixel x by the nonlocal means is
(6)
N
L
ν
(
x
)
≔
1
C
(
x
)
∫
Ω
w
ν
(
x
,
y
)
ν
(
y
)
d
y
,
where the normalized factor C
(
x
)
=
∫
Ω
w
ν
(
x
,
y
)
d
y. It can be seen that this method is a kind of special filter algorithm, and the estimation of the value at pixel x is weighted by the values of all the pixels y
(
y
∈
Ω
) in the image.Let the 2-dimensional image domain beΩ
⊂
R
2, pixels x
∈
Ω, and u
(
x
) is a real function on Ω
→
R. The derivation under the nonlocal frame is
(7)
∂
y
u
(
x
)
=
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
where w
:
Ω
×
Ω
→
R is the nonnegative symmetrical weight function defined as (4), which is calculated by the reference image.The nonlocal gradient∇
w
u
(
x
,
y
) is defined as a vector composed by all the partial derivations at position x: ∇
w
u
(
x
,
·
), where ∇
w
u
:
Ω
→
Ω
×
Ω:
(8)
∇
w
u
(
x
,
y
)
≔
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
∀
y
∈
Ω
.Define the nonlocal divergence as(9)
di
v
w
v
→
(
x
)
≔
∫
Ω
(
v
(
x
,
y
)
-
v
(
y
,
x
)
)
w
(
x
,
y
)
d
y
,
which satisfies the conjugation with the nonlocal gradient; that is,
(10)
〈
∇
w
u
,
v
→
〉
=
-
〈
u
,
di
v
w
v
→
〉
,
∀
u
:
Ω
⟶
R
,
∀
v
→
:
Ω
×
Ω
⟶
R
.Thus the nonlocal Laplacian is defined as(11)
Δ
w
u
(
x
)
:
=
1
2
di
v
w
(
∇
w
u
(
x
)
)
=
∫
Ω
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
d
y
.
The NL-TV norm of u
∈
L
1
(
Ω
) based on nonlocal operators, that is, the weighted gradient ∇
w
u, is defined by
(12)
J
NL-TV
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
=
∫
Ω
∫
Ω
(
u
(
x
)
-
u
(
y
)
)
2
w
(
x
,
y
)
d
y
d
x
,
where function w
(
x
,
y
) is defined as (4). The NL-TV model can be written as
(13)
∥
u
∥
L
1
(
Ω
)
+
J
NL-TV
(
u
)
.
The well-posedness of NL-TV model is guaranteed in [11].


### 2.3. Nonlocal Vectorial Total Variation Model and the Fixed Point Iterative Algorithm

In this section, the nonlocal gradient will be introduced to reconstruct the variation norm in the vectorial total variation model to deduce the NL-VTV model.The vectorial nonlocal gradient could be derived by generalizing the nonlocal gradient defined previously:(14)
∇
w
u
=
(
∇
w
u
1
,
∇
w
u
2
,
…
,
∇
w
u
M
)
,
where ∇
w
u
m
=
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
), m
=
1
,
…
,
M, represents the gradient of each channel separately, which is calculated by the image of its own channel.Use the vectorial nonlocal gradient∇
w
u to replace the gradient ∇
u in VTV norm; we get the NL-VTV model:
(15)
inf
⁡
u
∈
R
M
E
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
where ∫
Ω
|
∇
w
u
|
d
x is the VTV regularization based on the nonlocal gradient, which could be denoted by
(16)
J
(
u
)
=
NL-VTV
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
≔
∫
Ω
∑
m
=
1
M
∫
Ω
(
u
m
(
x
)
-
u
m
(
y
)
)
2
w
m
(
x
,
y
)
d
y
d
x
.The existence and the uniqueness of the solution to (15) have been proven in our previews work, in which we proposed a fixed point iterative algorithm to solve (15) as well. Here we will give out a brief introduction as follows.The Euler-Lagrangian equation of the functional (15) is
(17)
-
div
⁡
w
(
∇
w
u
|
∇
w
u
|
)
+
λ
(
u
-
u
0
)
=
0
.Calculate the nonlocal curvature channel by channel; that is,(18)
κ
w
m
≔
div
⁡
w
(
∇
w
u
m
|
∇
w
u
|
)
=
∫
Ω
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
)
×
(
1
|
∇
w
u
|
(
x
)
+
1
|
∇
w
u
|
(
y
)
)
d
y
,
where |
∇
w
u
|
(
q
)
≔
∑
m
=
1
M
∫
Ω
(
u
m
(
z
)
-
u
m
(
q
)
)
2
w
m
(
q
,
z
)
d
z, m
=
1
,
…
,
M. Then the Euler-Lagrange equation becomes
(19)
-
∫
Ω
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
)
[
1
|
∇
w
u
|
(
x
)
+
1
|
∇
w
u
|
(
y
)
]
d
y
+
λ
(
u
-
u
0
)
=
0
,
m
=
1
,
…
,
M
.According to (19), when the weight function is fixed, the Euler-Lagrange equation is linear, which could be solved by the gradient descent kind of algorithm as follows.First of all, we discretize (19). Denote u
i
m to be the pixel value at position i of the mth channel image, where 1
≤
m
≤
M, 1
≤
i
≤
N (if the size of the image is a
×
b, N
=
a
×
b). w
i
,
j
m is the discrete value of the weight function w
(
x
,
y
), that is, the weight function value with respect to pixel i and pixel j. The pixel j traverses the image channel where the pixel i is in. Denote the neighbor set to be N
i
:
=
{
j
:
w
i
,
j
>
τ
}, where τ denotes a threshold value as a judgment of the similarity between the two blocks centered at pixel i and pixel j, respectively. Let ∇
w
d be the discrete ∇
w:
(20)
∇
w
d
(
u
i
m
)
=
(
u
j
m
-
u
i
m
)
w
i
,
j
m
,
j
∈
N
i
.
And let di
v
w
d be the discrete di
v
w: di
v
w
d
(
p
i
,
j
)
≔
∑
j
∈
N
i
(
p
i
,
j
-
p
j
,
i
)
w
i
,
j
m.As for functions, the discrete inner product is〈
u
,
v
〉
≔
∑
i
(
u
i
v
i
), while, as for vectors, the discrete dot product is (
p
,
q
)
i
≔
∑
j
(
p
i
,
j
q
i
,
j
), the discrete inner product is 〈
p
,
q
〉
≔
∑
i
∑
j
(
p
i
,
j
q
i
,
j
). The norm of vector is |
p
i
|
=
∑
j
(
p
i
,
j
)
2.Then for channelm, where 1
≤
m
≤
M, the discrete nonlocal curvature is
(21)
κ
w
m
≔
div
⁡
w
m
(
∇
w
u
m
|
∇
w
u
|
)
=
∑
j
(
u
j
m
-
u
i
m
)
w
i
,
j
m
[
1
|
∇
w
u
|
(
i
)
+
1
|
∇
w
u
|
(
j
)
]
,
where |
∇
w
u
|
(
i
)
=
∑
m
=
1
M
∑
j
(
u
j
m
-
u
i
m
)
2
w
i
,
j
m and |
∇
w
u
|
(
j
)
=
∑
m
=
1
M
∑
i
(
u
j
m
-
u
i
m
)
2
w
i
,
j
m. Then the Euler-Lagrange equation becomes
(22)
-
∑
j
(
u
j
m
-
u
i
m
)
w
i
,
j
m
[
1
|
∇
w
u
|
(
i
)
+
1
|
∇
w
u
|
(
j
)
]
+
λ
(
u
i
m
-
u
0
i
m
)
=
0
,
m
=
1
,
…
,
M
.Denotew
i
,
j
m
[
(
1
/
|
∇
w
u
|
(
i
)
)
+
(
1
/
|
∇
w
u
|
(
j
)
)
]
≜
W
i
,
j
m. Then the Euler-Lagrange equation could be shortly written as
(23)
-
∑
j
W
i
,
j
m
(
u
j
m
-
u
i
m
)
+
λ
(
u
i
m
-
u
0
i
m
)
=
0
,
m
=
1
,
…
,
M
.
Furthermore, it could be written as
(24)
u
i
m
=
∑
j
W
i
,
j
m
∑
j
W
i
,
j
m
+
λ
u
j
m
+
λ
∑
j
W
i
,
j
m
+
λ
u
0
i
m
,
m
=
1
,
…
,
M
.According to the fixed point iteration, we use the multichannel imageu
(
k
) in step k to update image u
(
k
+
1
) in step k
+
1. The weights are calculated by the original observation image. Then we get the iteration formula as
(25)
u
(
0
)
=
u
0
,
u
i
(
k
+
1
)
=
∑
j
W
i
,
j
m
(
k
)
∑
j
W
i
,
j
m
(
k
)
+
λ
u
j
(
k
)
+
λ
∑
j
W
i
,
j
m
(
k
)
+
λ
u
0
i
,
where i
∈
Ω, j
∈
N
i.The convergence property of the fixed point algorithm has been proved in our in-press work. However, when|
∇
w
u
|
=
0, the NLTV functional is nondifferentiable, which is similar to the classical TV model. Thus a regularization |
∇
w
u
|
2
+
ϵ is involved in the calculation of W
i
,
j
m to escape the denominator to be zero. Consequently, this algorithm takes many iteration steps to obtain the best restoration result. Thus we consider designing a fast alternating minimization algorithm instead.


## 2.1. Vectorial Total Variation Model

Consider a vector valued functionu, such as the multichannel images, defined on a bounded open domain Ω
⊂
R
N:
(1)
u
:
Ω
⟶
R
M
x
⟼
u
(
x
)
=
(
u
1
(
x
)
,
u
2
(
x
)
,
…
,
u
M
(
x
)
)
,
where u
m
:
Ω
→
R
(
1
≤
m
≤
M
) denotes the mth channel of u.The regularized VTV model for multichannel image reconstruction could be written as(2)
inf
⁡
u
∈
R
M
E
(
u
)
≔
∫
Ω
|
∇
u
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
where u
0 is the observed noisy multichannel image and λ
>
0 is the regularized parameter. The vectorial gradient of u is defined as
(3)
∇
u
=
(
∇
u
1
,
…
,
∇
u
M
)
=
[
∂
u
1
∂
x
∂
u
2
∂
x
⋯
∂
u
M
∂
x
∂
u
1
∂
y
∂
u
2
∂
y
⋯
∂
u
M
∂
y
]
∈
R
2
×
M
,
|
∇
u
|
=
∥
∇
u
∥
F
=
∑
m
=
1
M
(
∂
u
m
∂
x
)
2
+
∑
m
=
1
M
(
∂
u
m
∂
y
)
2
=
∑
m
=
1
M
|
∇
u
m
|
2
.
We refer to [6] for more information about the well-posedness analysis of the VTV model.


## 2.2. Nonlocal Functional

The neighborhood filter NLM proposed by Buades et al. developed to a famous image denoising method [3]. It uses the pixels similar to the current pixel to estimate the value of the current pixel. Let the reference image be ν. The similarity of the two pixels in the image x, y could be measured as the following weight function w
ν
(
x
,
y
):
(4)
w
ν
(
x
,
y
)
=
exp
⁡
(
-
(
G
a
*
|
ν
(
x
+
·
)
-
ν
(
y
+
·
)
|
2
)
(
0
)
h
2
)
,
where G
a is the Gaussian kernel with the standard deviation of a:
(5)
(
G
a
*
|
ν
(
x
+
·
)
-
ν
(
y
+
·
)
|
2
)
(
0
)
=
∫
R
2
G
a
(
t
)
|
ν
(
x
+
t
)
-
ν
(
y
+
t
)
|
2
d
t
,
where h is the filter parameter. Generally, h corresponds to the noise level, which is usually set as the standard deviation of the noise. It is shown that the weight function above could be calculated by the image blocks centered around x, y. The value of the weight function reflects the similarity of the two blocks. The estimation of the value at pixel x by the nonlocal means is
(6)
N
L
ν
(
x
)
≔
1
C
(
x
)
∫
Ω
w
ν
(
x
,
y
)
ν
(
y
)
d
y
,
where the normalized factor C
(
x
)
=
∫
Ω
w
ν
(
x
,
y
)
d
y. It can be seen that this method is a kind of special filter algorithm, and the estimation of the value at pixel x is weighted by the values of all the pixels y
(
y
∈
Ω
) in the image.Let the 2-dimensional image domain beΩ
⊂
R
2, pixels x
∈
Ω, and u
(
x
) is a real function on Ω
→
R. The derivation under the nonlocal frame is
(7)
∂
y
u
(
x
)
=
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
where w
:
Ω
×
Ω
→
R is the nonnegative symmetrical weight function defined as (4), which is calculated by the reference image.The nonlocal gradient∇
w
u
(
x
,
y
) is defined as a vector composed by all the partial derivations at position x: ∇
w
u
(
x
,
·
), where ∇
w
u
:
Ω
→
Ω
×
Ω:
(8)
∇
w
u
(
x
,
y
)
≔
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
∀
y
∈
Ω
.Define the nonlocal divergence as(9)
di
v
w
v
→
(
x
)
≔
∫
Ω
(
v
(
x
,
y
)
-
v
(
y
,
x
)
)
w
(
x
,
y
)
d
y
,
which satisfies the conjugation with the nonlocal gradient; that is,
(10)
〈
∇
w
u
,
v
→
〉
=
-
〈
u
,
di
v
w
v
→
〉
,
∀
u
:
Ω
⟶
R
,
∀
v
→
:
Ω
×
Ω
⟶
R
.Thus the nonlocal Laplacian is defined as(11)
Δ
w
u
(
x
)
:
=
1
2
di
v
w
(
∇
w
u
(
x
)
)
=
∫
Ω
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
d
y
.
The NL-TV norm of u
∈
L
1
(
Ω
) based on nonlocal operators, that is, the weighted gradient ∇
w
u, is defined by
(12)
J
NL-TV
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
=
∫
Ω
∫
Ω
(
u
(
x
)
-
u
(
y
)
)
2
w
(
x
,
y
)
d
y
d
x
,
where function w
(
x
,
y
) is defined as (4). The NL-TV model can be written as
(13)
∥
u
∥
L
1
(
Ω
)
+
J
NL-TV
(
u
)
.
The well-posedness of NL-TV model is guaranteed in [11].


## 2.3. Nonlocal Vectorial Total Variation Model and the Fixed Point Iterative Algorithm

In this section, the nonlocal gradient will be introduced to reconstruct the variation norm in the vectorial total variation model to deduce the NL-VTV model.The vectorial nonlocal gradient could be derived by generalizing the nonlocal gradient defined previously:(14)
∇
w
u
=
(
∇
w
u
1
,
∇
w
u
2
,
…
,
∇
w
u
M
)
,
where ∇
w
u
m
=
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
), m
=
1
,
…
,
M, represents the gradient of each channel separately, which is calculated by the image of its own channel.Use the vectorial nonlocal gradient∇
w
u to replace the gradient ∇
u in VTV norm; we get the NL-VTV model:
(15)
inf
⁡
u
∈
R
M
E
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
where ∫
Ω
|
∇
w
u
|
d
x is the VTV regularization based on the nonlocal gradient, which could be denoted by
(16)
J
(
u
)
=
NL-VTV
(
u
)
=
∫
Ω
|
∇
w
u
|
d
x
≔
∫
Ω
∑
m
=
1
M
∫
Ω
(
u
m
(
x
)
-
u
m
(
y
)
)
2
w
m
(
x
,
y
)
d
y
d
x
.The existence and the uniqueness of the solution to (15) have been proven in our previews work, in which we proposed a fixed point iterative algorithm to solve (15) as well. Here we will give out a brief introduction as follows.The Euler-Lagrangian equation of the functional (15) is
(17)
-
div
⁡
w
(
∇
w
u
|
∇
w
u
|
)
+
λ
(
u
-
u
0
)
=
0
.Calculate the nonlocal curvature channel by channel; that is,(18)
κ
w
m
≔
div
⁡
w
(
∇
w
u
m
|
∇
w
u
|
)
=
∫
Ω
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
)
×
(
1
|
∇
w
u
|
(
x
)
+
1
|
∇
w
u
|
(
y
)
)
d
y
,
where |
∇
w
u
|
(
q
)
≔
∑
m
=
1
M
∫
Ω
(
u
m
(
z
)
-
u
m
(
q
)
)
2
w
m
(
q
,
z
)
d
z, m
=
1
,
…
,
M. Then the Euler-Lagrange equation becomes
(19)
-
∫
Ω
(
u
m
(
y
)
-
u
m
(
x
)
)
w
m
(
x
,
y
)
[
1
|
∇
w
u
|
(
x
)
+
1
|
∇
w
u
|
(
y
)
]
d
y
+
λ
(
u
-
u
0
)
=
0
,
m
=
1
,
…
,
M
.According to (19), when the weight function is fixed, the Euler-Lagrange equation is linear, which could be solved by the gradient descent kind of algorithm as follows.First of all, we discretize (19). Denote u
i
m to be the pixel value at position i of the mth channel image, where 1
≤
m
≤
M, 1
≤
i
≤
N (if the size of the image is a
×
b, N
=
a
×
b). w
i
,
j
m is the discrete value of the weight function w
(
x
,
y
), that is, the weight function value with respect to pixel i and pixel j. The pixel j traverses the image channel where the pixel i is in. Denote the neighbor set to be N
i
:
=
{
j
:
w
i
,
j
>
τ
}, where τ denotes a threshold value as a judgment of the similarity between the two blocks centered at pixel i and pixel j, respectively. Let ∇
w
d be the discrete ∇
w:
(20)
∇
w
d
(
u
i
m
)
=
(
u
j
m
-
u
i
m
)
w
i
,
j
m
,
j
∈
N
i
.
And let di
v
w
d be the discrete di
v
w: di
v
w
d
(
p
i
,
j
)
≔
∑
j
∈
N
i
(
p
i
,
j
-
p
j
,
i
)
w
i
,
j
m.As for functions, the discrete inner product is〈
u
,
v
〉
≔
∑
i
(
u
i
v
i
), while, as for vectors, the discrete dot product is (
p
,
q
)
i
≔
∑
j
(
p
i
,
j
q
i
,
j
), the discrete inner product is 〈
p
,
q
〉
≔
∑
i
∑
j
(
p
i
,
j
q
i
,
j
). The norm of vector is |
p
i
|
=
∑
j
(
p
i
,
j
)
2.Then for channelm, where 1
≤
m
≤
M, the discrete nonlocal curvature is
(21)
κ
w
m
≔
div
⁡
w
m
(
∇
w
u
m
|
∇
w
u
|
)
=
∑
j
(
u
j
m
-
u
i
m
)
w
i
,
j
m
[
1
|
∇
w
u
|
(
i
)
+
1
|
∇
w
u
|
(
j
)
]
,
where |
∇
w
u
|
(
i
)
=
∑
m
=
1
M
∑
j
(
u
j
m
-
u
i
m
)
2
w
i
,
j
m and |
∇
w
u
|
(
j
)
=
∑
m
=
1
M
∑
i
(
u
j
m
-
u
i
m
)
2
w
i
,
j
m. Then the Euler-Lagrange equation becomes
(22)
-
∑
j
(
u
j
m
-
u
i
m
)
w
i
,
j
m
[
1
|
∇
w
u
|
(
i
)
+
1
|
∇
w
u
|
(
j
)
]
+
λ
(
u
i
m
-
u
0
i
m
)
=
0
,
m
=
1
,
…
,
M
.Denotew
i
,
j
m
[
(
1
/
|
∇
w
u
|
(
i
)
)
+
(
1
/
|
∇
w
u
|
(
j
)
)
]
≜
W
i
,
j
m. Then the Euler-Lagrange equation could be shortly written as
(23)
-
∑
j
W
i
,
j
m
(
u
j
m
-
u
i
m
)
+
λ
(
u
i
m
-
u
0
i
m
)
=
0
,
m
=
1
,
…
,
M
.
Furthermore, it could be written as
(24)
u
i
m
=
∑
j
W
i
,
j
m
∑
j
W
i
,
j
m
+
λ
u
j
m
+
λ
∑
j
W
i
,
j
m
+
λ
u
0
i
m
,
m
=
1
,
…
,
M
.According to the fixed point iteration, we use the multichannel imageu
(
k
) in step k to update image u
(
k
+
1
) in step k
+
1. The weights are calculated by the original observation image. Then we get the iteration formula as
(25)
u
(
0
)
=
u
0
,
u
i
(
k
+
1
)
=
∑
j
W
i
,
j
m
(
k
)
∑
j
W
i
,
j
m
(
k
)
+
λ
u
j
(
k
)
+
λ
∑
j
W
i
,
j
m
(
k
)
+
λ
u
0
i
,
where i
∈
Ω, j
∈
N
i.The convergence property of the fixed point algorithm has been proved in our in-press work. However, when|
∇
w
u
|
=
0, the NLTV functional is nondifferentiable, which is similar to the classical TV model. Thus a regularization |
∇
w
u
|
2
+
ϵ is involved in the calculation of W
i
,
j
m to escape the denominator to be zero. Consequently, this algorithm takes many iteration steps to obtain the best restoration result. Thus we consider designing a fast alternating minimization algorithm instead.


## 3. The Alternating Minimization Algorithm

For the sake of conquering the nondifferentiability of the VTV functional, we design an alternating minimization algorithm in this section to solve the problem (15). The main idea of the alternating minimization is introducing an auxiliary variable. Then the original problem could be split into two subproblems of which one is differential with respect to the original variable, while the other is convex with respect to the auxiliary variable, respectively.Yang et al. [10] have already studied the alternating minimization algorithm for a general model:
(26)
min
⁡
u
∑
i
α
i
∥
G
i
u
∥
+
μ
2
∥
K
u
-
f
∥
2
.
Model (26) reduces to (15) when G
i
=
I
m
⊗
∇
w, α
i
=
1, and K
=
I, where I
m is the identity matrix of order m and ⊗ represents the Kronecker product. Yang et al. pointed out in [10] that if G
i is the finite difference operator D
i, then, under some boundary conditions for u and some other assumptions for the blurring kernel K, the problem (26) can be fast transformed by fast Fourier transform, which diagonalizes the block circulant matrices: D
i, K, and their transpose, and then can be solved by the fast algorithm such as the Gaussian elimination method.However, the nonlocal gradient operator∇
w in (15) is not block circulant. We will figure out another method to speed up the algorithm. We briefly propose the algorithm in this section. The convergence of this algorithm will be analysed in Section 4, as well as others properties of this algorithm.


### 3.1. Variable Splitting (Half-Quadratic Formulation of (15))

We introduce an auxiliary variabled to approximate the nonlocal gradient of u, that is, ∇
w
u. Then the minimization problem (15) is transformed into an equivalent constrained problem:
(27)
min
⁡
u
∈
R
M
E
(
u
)
=
∫
Ω
|
d
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
subject
to
d
=
∇
w
u
.In order to solve problem (27), we use the penalty method to reformulate it as
(28)
min
⁡
u
,
d
E
(
u
,
d
)
=
∫
Ω
|
d
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
,
where β is the penalty parameter. It could be easily proved that, as β
→
∞, the solution of problem (28) is convergent to the solution of problem (15).


### 3.2. Alternating Minimization Algorithm

The problem of minimizing (28) with respect to d for fixed u becomes
(29)
min
⁡
d
∫
Ω
|
d
|
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
.
The functional in (29) is convex with respect to d. According to Lemma  3.3 in [10], the problem (29) has the unique solutions given by
(30)
d
k
+
1
=
shrink
(
∇
w
u
k
+
1
,
1
β
)
,
where the multidimensional shrinkage operator shrink
(
x
,
λ
)
=
(
x
/
|
x
|
)
·
max
⁡
{
|
x
|
-
λ
,
0
}.On the other hand, minimizing (28) with respect to u for fixed d becomes
(31)
min
⁡
u
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
.
The functional in (31) is differential with respect to u. Thus the solution to (31) satisfies the first order optimized condition; that is, the Euler-Lagrange equation of the functional (31) is
(32)
λ
(
u
-
u
0
)
+
β
div
⁡
w
(
d
-
∇
w
u
)
=
0
.
After some steps of simple calculation, we get
(33)
u
-
β
λ
div
⁡
w
(
∇
w
u
)
=
u
0
-
β
λ
div
⁡
w
d
.
Then according to (11), we get
(34)
(
I
-
2
β
λ
Δ
w
)
u
=
u
0
-
β
λ
div
⁡
w
d
.
Finally, the multichannel image u is calculated as follows:
(35)
(
I
-
2
β
λ
Δ
w
)
u
k
+
1
=
u
0
-
β
λ
div
⁡
w
d
k
.Since the graph LaplacianΔ
w is negative semidefinite, the operator I
-
(
2
β
/
λ
)
Δ
w is diagonally dominant. Therefore we can solve u
k
+
1 by a Gauss-Seidel algorithm. However, the subproblem (35) consists in solving the linear equations with the scale of N
×
N for each image channel. As a result, the calculation complexity of the Gauss-Seidel algorithm is heavy. We are inspired by Chen and Cheng [12] to simplify, and, without loss of generality as well, we chose λ
=
β
2. Furthermore, using Taylor expansion to approximate, there is
(36)
(
I
-
2
β
λ
Δ
w
)
-
1
=
(
I
-
2
β
Δ
w
)
-
1
=
I
+
2
β
Δ
w
+
O
(
β
-
2
)
≈
I
+
2
β
Δ
w
=
I
+
2
β
(
W
-
I
)
=
(
1
-
2
β
)
I
+
2
β
W
,
where W
=
(
w
i
j
)
i
,
j
=
1
n
2 is the weight matrix of similar blocks of each channel, which is normalized to a symmetrical matrix with the sum of every row equal to 1. Thus the solution to (35) is
(37)
u
k
+
1
=
(
1
-
2
β
)
I
(
u
0
-
1
β
div
⁡
w
d
k
)
+
2
β
W
(
u
0
-
1
β
div
⁡
w
d
k
)
.
According to (9), the divergence at pixel i is div
⁡
w
d
i
m
=
∑
j
(
d
i
j
m
-
d
j
i
m
)
w
i
j, m
=
1
,
…
,
M.To conclude this section, we rewrite the alternating minimization algorithm as in Table1.Table 1
Comparison of the highest SNR obtained.
SNR
Noisy image
VTV
NLVTV
Barbara image
8.9013
13.1212
13.2873
Baboon image
9.2373
13.1618
12.5936


## 3.1. Variable Splitting (Half-Quadratic Formulation of (15))

We introduce an auxiliary variabled to approximate the nonlocal gradient of u, that is, ∇
w
u. Then the minimization problem (15) is transformed into an equivalent constrained problem:
(27)
min
⁡
u
∈
R
M
E
(
u
)
=
∫
Ω
|
d
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
,
subject
to
d
=
∇
w
u
.In order to solve problem (27), we use the penalty method to reformulate it as
(28)
min
⁡
u
,
d
E
(
u
,
d
)
=
∫
Ω
|
d
|
d
x
+
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
,
where β is the penalty parameter. It could be easily proved that, as β
→
∞, the solution of problem (28) is convergent to the solution of problem (15).


## 3.2. Alternating Minimization Algorithm

The problem of minimizing (28) with respect to d for fixed u becomes
(29)
min
⁡
d
∫
Ω
|
d
|
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
.
The functional in (29) is convex with respect to d. According to Lemma  3.3 in [10], the problem (29) has the unique solutions given by
(30)
d
k
+
1
=
shrink
(
∇
w
u
k
+
1
,
1
β
)
,
where the multidimensional shrinkage operator shrink
(
x
,
λ
)
=
(
x
/
|
x
|
)
·
max
⁡
{
|
x
|
-
λ
,
0
}.On the other hand, minimizing (28) with respect to u for fixed d becomes
(31)
min
⁡
u
λ
2
∫
Ω
|
u
-
u
0
|
2
d
x
+
β
2
∫
Ω
|
d
-
∇
w
u
|
2
d
x
.
The functional in (31) is differential with respect to u. Thus the solution to (31) satisfies the first order optimized condition; that is, the Euler-Lagrange equation of the functional (31) is
(32)
λ
(
u
-
u
0
)
+
β
div
⁡
w
(
d
-
∇
w
u
)
=
0
.
After some steps of simple calculation, we get
(33)
u
-
β
λ
div
⁡
w
(
∇
w
u
)
=
u
0
-
β
λ
div
⁡
w
d
.
Then according to (11), we get
(34)
(
I
-
2
β
λ
Δ
w
)
u
=
u
0
-
β
λ
div
⁡
w
d
.
Finally, the multichannel image u is calculated as follows:
(35)
(
I
-
2
β
λ
Δ
w
)
u
k
+
1
=
u
0
-
β
λ
div
⁡
w
d
k
.Since the graph LaplacianΔ
w is negative semidefinite, the operator I
-
(
2
β
/
λ
)
Δ
w is diagonally dominant. Therefore we can solve u
k
+
1 by a Gauss-Seidel algorithm. However, the subproblem (35) consists in solving the linear equations with the scale of N
×
N for each image channel. As a result, the calculation complexity of the Gauss-Seidel algorithm is heavy. We are inspired by Chen and Cheng [12] to simplify, and, without loss of generality as well, we chose λ
=
β
2. Furthermore, using Taylor expansion to approximate, there is
(36)
(
I
-
2
β
λ
Δ
w
)
-
1
=
(
I
-
2
β
Δ
w
)
-
1
=
I
+
2
β
Δ
w
+
O
(
β
-
2
)
≈
I
+
2
β
Δ
w
=
I
+
2
β
(
W
-
I
)
=
(
1
-
2
β
)
I
+
2
β
W
,
where W
=
(
w
i
j
)
i
,
j
=
1
n
2 is the weight matrix of similar blocks of each channel, which is normalized to a symmetrical matrix with the sum of every row equal to 1. Thus the solution to (35) is
(37)
u
k
+
1
=
(
1
-
2
β
)
I
(
u
0
-
1
β
div
⁡
w
d
k
)
+
2
β
W
(
u
0
-
1
β
div
⁡
w
d
k
)
.
According to (9), the divergence at pixel i is div
⁡
w
d
i
m
=
∑
j
(
d
i
j
m
-
d
j
i
m
)
w
i
j, m
=
1
,
…
,
M.To conclude this section, we rewrite the alternating minimization algorithm as in Table1.Table 1
Comparison of the highest SNR obtained.
SNR
Noisy image
VTV
NLVTV
Barbara image
8.9013
13.1212
13.2873
Baboon image
9.2373
13.1618
12.5936


## 4. Convergence Analysis

Since it has been proved in our in-press work that the model (15) has the unique solution, in this section, we prove that the solution obtained by Algorithm 1 is the unique one. First, the convergence of the quadratic penalty method as the penalty parameter goes to infinity is well known [13]. That is, as β
→
∞, the solution of (31) converges to that of (15). We analyze the convergence properties of Algorithm 1 for a fixed β
>
0. In this section, we show that, for any fixed β, the algorithm of minimizing u and d alternately has finite convergence for some variables and q-linear convergence for the rest.Algorithm 1:Alternating minimization algorithm for nonlocal VTV denoising.
Initialization:
fork
=
0 to K, do
Solved
k
+
1
=
shrink
(
∇
w
u
k
+
1
,
1
/
β
).
Solveu
k
+
1
=
(
1
-
(
2
/
β
)
)
I
(
u
0
-
(
1
/
β
)
di
v
w
d
k
)
+
(
2
/
β
)
W
(
u
0
-
(
1
/
β
)
di
v
w
d
k
).
end forTheorem 1.
The iterative sequence{
(
d
k
,
u
k
)
} generated by (30) and (37) from any starting point {
(
d
0
,
u
0
)
} converges to the solution of the minimization problem (31).Proof.
According to (30) and (37), the functional E
(
u
,
d
) in (28) satisfies
(38)
E
(
u
k
+
1
,
d
k
+
1
)
≤
E
(
u
k
,
d
k
+
1
)
≤
E
(
u
k
,
d
k
)
.
Therefore, the sequence{
E
(
u
k
,
d
k
)
} is nonincreasing. Thus it is convergent in R
N because of its nonnegativity as well. We denote the limit of the sequence by m. We prove that m
=
E
(
u
*
,
d
*
) followed, where u
* and d
* are the solution of (31).
BecauseE is coercive and the sequence {
E
(
u
k
,
d
k
)
} is convergent, the sequence {
(
u
k
,
d
k
)
} is bounded. Thus there exists a subsequence {
(
u
k
,
d
k
)
}, which converges to an accumulation point (
u
^
,
d
^
) when n
k
→
∞. According to (30) and (37), for all n
k
∈
N, d,
(39)
E
(
u
n
k
,
d
n
k
+
1
)
≤
E
(
u
n
k
,
d
)
,
and, for all n
k
∈
N, u,
(40)
E
(
u
n
k
+
1
,
d
n
k
+
1
)
≤
E
(
u
,
d
n
k
+
1
)
.
We denoted
~ as the accumulation point of the sequence {
d
n
k
+
1
}, because of that the functional E is continuous, according to (38),
(41)
m
=
E
(
u
^
,
d
~
)
=
E
(
u
^
,
d
^
)
.
Because of that the threshold function, shrink, in (30), is continuous, we compute the limits of both sides of the equation and get
(42)
d
~
=
shrink
(
∇
w
u
~
,
1
β
)
.
Thusd
~
=
arg
min
⁡
d
E
(
u
^
,
d
). According to (41), d
^
=
arg
min
⁡
d
E
(
u
^
,
d
). On the other hand, E is strictly convex with respect to d, then we have d
~
=
d
^ because of the uniqueness of the minimizer of min
⁡
d
E
(
u
^
,
d
). Similarly, we prove that u
~
=
u
^.
We compute the limits of (39) and (40),
(43)
inf
⁡
d
E
(
u
^
,
d
)
=
E
(
u
^
,
d
^
)
≤
E
(
u
^
,
d
)
inf
⁡
u
E
(
u
,
d
^
)
=
E
(
u
^
,
d
^
)
≤
E
(
u
,
d
^
)
.
According to (43), (
0
,
0
)
T
∈
∂
E
(
u
^
,
d
^
). Thus the limit points of all of the subsequence of {
(
u
k
,
d
k
)
} are the solutions to (28). As a result, {
(
u
k
,
d
k
)
} converges to (
u
*
,
d
*
).LetS
1
/
β
(
∇
w
u
)
=
shrink
(
∇
w
u
,
1
/
β
) and let P
=
(
I
-
(
2
β
/
λ
)
Δ
w
)
-
1. We define set H
=
{
(
y
,
y
′
)
∈
R
m
×
n
2
⊗
R
m
×
n
2
∣
∥
W
·
(
y
-
y
′
)
∥
2
≤
ν
∥
y
-
y
′
∥
2
}, where ν
<
1 is a positive constant. Then we have the following conclusion.Theorem 2 (q-linear convergence).
Assume that the sequence{
u
k
}
⊆
H; then there is a constant β
0
>
0, such that for all β
>
β
0, the sequence {
u
k
} generated by (30) and (37) isq-linear convergence.Proof.
Combine (30) and (37); we have
(44)
u
k
+
1
=
P
(
u
0
-
1
β
div
⁡
w
(
S
1
/
β
(
∇
w
u
k
)
)
)
.
We know that u
* is a fixed point of (44). Let Q
=
P
-
1
=
I
-
(
2
β
/
λ
)
Δ
w
=
(
1
+
(
2
/
β
)
)
I
-
(
2
/
β
)
W. Then
(45)
Q
(
u
k
+
1
-
u
*
)
=
-
1
β
div
⁡
w
(
S
1
/
β
(
∇
w
u
k
)
-
S
1
/
β
(
∇
w
u
*
)
)
.
Since the operator S
1
/
β is nonexpansive,
(46)
∥
-
1
β
div
⁡
w
(
S
1
/
β
(
∇
w
u
k
)
-
S
1
/
β
(
∇
w
u
*
)
)
∥
2
≤
∥
div
⁡
w
∥
2
∥
∇
w
∥
2
β
∥
u
k
-
u
*
∥
2
.
Because {
u
k
}
⊆
H,
(47)
∥
Q
(
u
k
+
1
-
u
*
)
∥
2
=
∥
(
(
1
+
2
β
)
I
-
2
β
W
)
(
u
k
+
1
-
u
*
)
∥
2
≥
(
1
+
2
β
)
∥
u
k
+
1
-
u
*
∥
2
-
2
β
∥
W
(
u
k
+
1
-
u
*
)
∥
2
≥
(
1
+
2
β
-
2
β
ν
)
∥
u
k
+
1
-
u
*
∥
2
.
We writeτ
=
1
/
(
1
+
(
2
/
β
)
-
(
2
/
β
)
ν
)
<
1. Then by uniting (46) and (47), we have
(48)
∥
u
k
+
1
-
u
*
∥
2
≤
τ
∥
div
⁡
w
∥
2
∥
∇
w
∥
2
β
∥
u
k
-
u
*
∥
2
.
Letβ
0
=
τ
∥
div
⁡
w
∥
2
∥
∇
w
∥
2
>
0. Then for all β
>
β
0, the sequence {
u
k
} generated by (30) and (37) is q-linear convergent to u
*.


## 5. Experiments

In this section, we apply Algorithm1 to color image denoising to test its efficiency. We compare this method to that in [10] to show the detail-preserving property of the nonlocal VTV model. Then we compare the convergent rate of this algorithm with the fixed point algorithm presented in Section 2.


### 5.1. Texture-Preserving Property of the NLVTV Model

We choose the color Barbara picture size300
×
300
×
3 and the color Baboon picture size 256
×
256
×
3 for testing, both of which contain many textures, such as the woman’s scarf and the baboon’s beard. We add the Gaussian white noise with standard deviation 20 to the two original pictures. Then we use the VTV model and the NLVTV model to restore them, separately. In this experiment, the signal to noise rate (SNR) is used to estimate the restoration effect of the methods, which is defined as follows:
(49)
SNR
=
10
·
lo
g
10
(
var
(
u
)
var
(
u
k
-
u
)
)
,
where var
(
x
) represents the variance of all the elements of matrix x.As is shown in Figure1 and Table 1, the NLVTV model obtains similar SNR to VTV model. In the experiment of restoring the color Barbara picture, we choose the restored image by NLVTV model under the SNR (13.0678) lower than but close to that of the VTV model (SNR = 13.1212) instead of the restored image under the best SNR (13.2873). Figure 2 shows the results. Figure 3 shows the zoom part of the restored image, that is, the woman’s scarf. We can see that the NVTV model preserves the details better than the VTV model. Figures 4 and 5 show the results of the experiment of restoring the color Baboon picture. Although the NLVTV model gets lower SNR than the VTV model, the textures in this picture are better preserved by the NLVTV model rather than the VTV model.Figure 1
Comparison of SNR.Comparison between the NLVTV and the VTV in the Barbara image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV (SNR = 13.1212)
(d)
Denoised by NLVTV (SNR = 13.0678)Zoom in the Barbara image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV
(d)
Denoised by NLVTVComparison between the NLVTV and the VTV in the Baboon image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV (SNR=13.1618)
(d)
Denoised by NVTV (SNR=12.5936)Zoom in the Baboon image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV
(d)
Denoised by NVTV


### 5.2. Computation Complexity Analysis and the Convergent Rate Comparison

In this section, we compare the computation complexity and the convergent rate of this algorithm to the fixed point algorithm.First, we analyze the computation complexities of the two algorithms, respectively. Assume the size of the multichannel image isn
×
n
×
M. The size of the nonlocal searching window is m
×
m. We do not consider the computation complexity of the calculation of the nonlocal weight, because that is executed exactly the same in both algorithms.At each iterative step of Algorithm1, namely, (30) and (37), the dominant calculation of (30) contains the calculation of the nonlocal gradient and the threshold computation of M
(
n
m
)
2 times. Thus the computation complexity is O
(
M
n
4
)
+
O
(
M
(
n
m
)
2
). The dominant calculation of (37) contains a step of nonlocal divergence computation with the computation complexity of O
(
M
n
2
) and a matrix multiplication with the computation complexity of O
(
M
n
4
). As a result, the total computation complexity for each iteration is O
(
M
n
4
)
+
O
(
M
(
n
m
)
2
)
+
O
(
M
n
2
)
+
O
(
M
n
4
).On the other hand, each iterative step of the fixed point iteration contains the nonlocal curvature with the computation complexity ofO
(
M
n
4
) and three matrix multiplications with the computation complexity of O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
). As a result, the total computation complexity for each iteration is O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
).Consequently, the computation complexity of the proposed algorithm is much less than the fixed point algorithm.The convergence rate can be represented by the relationship between the iterative step and the standard deviation of the restored image, intuitively. Figure6 shows the comparison of the convergence rate between the fixed point iteration and the proposed alternating minimization algorithm. The horizontal axis represents the iterative step. And the vertical axis represents the standard deviation of the restored image in each step. The standard deviation is calculated as
(50)
std
(
|
u
k
-
u
k
-
1
|
)
,
where u
k and u
k
-
1 represent the restored images of step k and step k
-
1, respectively. When calculating std
(
x
), we reshape matrix x into a vector x whose elements are taken column-wise from matrix x. Then
(51)
std
(
x
)
=
(
1
n
-
1
∑
i
=
1
n
(
x
i
-
x
-
)
2
)
1
/
2
.
Figure 6 shows that the alternating minimization algorithm converges faster.Figure 6
Comparison of the convergence rate.


## 5.1. Texture-Preserving Property of the NLVTV Model

We choose the color Barbara picture size300
×
300
×
3 and the color Baboon picture size 256
×
256
×
3 for testing, both of which contain many textures, such as the woman’s scarf and the baboon’s beard. We add the Gaussian white noise with standard deviation 20 to the two original pictures. Then we use the VTV model and the NLVTV model to restore them, separately. In this experiment, the signal to noise rate (SNR) is used to estimate the restoration effect of the methods, which is defined as follows:
(49)
SNR
=
10
·
lo
g
10
(
var
(
u
)
var
(
u
k
-
u
)
)
,
where var
(
x
) represents the variance of all the elements of matrix x.As is shown in Figure1 and Table 1, the NLVTV model obtains similar SNR to VTV model. In the experiment of restoring the color Barbara picture, we choose the restored image by NLVTV model under the SNR (13.0678) lower than but close to that of the VTV model (SNR = 13.1212) instead of the restored image under the best SNR (13.2873). Figure 2 shows the results. Figure 3 shows the zoom part of the restored image, that is, the woman’s scarf. We can see that the NVTV model preserves the details better than the VTV model. Figures 4 and 5 show the results of the experiment of restoring the color Baboon picture. Although the NLVTV model gets lower SNR than the VTV model, the textures in this picture are better preserved by the NLVTV model rather than the VTV model.Figure 1
Comparison of SNR.Comparison between the NLVTV and the VTV in the Barbara image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV (SNR = 13.1212)
(d)
Denoised by NLVTV (SNR = 13.0678)Zoom in the Barbara image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV
(d)
Denoised by NLVTVComparison between the NLVTV and the VTV in the Baboon image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV (SNR=13.1618)
(d)
Denoised by NVTV (SNR=12.5936)Zoom in the Baboon image.
(a)
Original image
(b)
Noisy image
(c)
Denoised by VTV
(d)
Denoised by NVTV


## 5.2. Computation Complexity Analysis and the Convergent Rate Comparison

In this section, we compare the computation complexity and the convergent rate of this algorithm to the fixed point algorithm.First, we analyze the computation complexities of the two algorithms, respectively. Assume the size of the multichannel image isn
×
n
×
M. The size of the nonlocal searching window is m
×
m. We do not consider the computation complexity of the calculation of the nonlocal weight, because that is executed exactly the same in both algorithms.At each iterative step of Algorithm1, namely, (30) and (37), the dominant calculation of (30) contains the calculation of the nonlocal gradient and the threshold computation of M
(
n
m
)
2 times. Thus the computation complexity is O
(
M
n
4
)
+
O
(
M
(
n
m
)
2
). The dominant calculation of (37) contains a step of nonlocal divergence computation with the computation complexity of O
(
M
n
2
) and a matrix multiplication with the computation complexity of O
(
M
n
4
). As a result, the total computation complexity for each iteration is O
(
M
n
4
)
+
O
(
M
(
n
m
)
2
)
+
O
(
M
n
2
)
+
O
(
M
n
4
).On the other hand, each iterative step of the fixed point iteration contains the nonlocal curvature with the computation complexity ofO
(
M
n
4
) and three matrix multiplications with the computation complexity of O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
). As a result, the total computation complexity for each iteration is O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
)
+
O
(
M
n
4
).Consequently, the computation complexity of the proposed algorithm is much less than the fixed point algorithm.The convergence rate can be represented by the relationship between the iterative step and the standard deviation of the restored image, intuitively. Figure6 shows the comparison of the convergence rate between the fixed point iteration and the proposed alternating minimization algorithm. The horizontal axis represents the iterative step. And the vertical axis represents the standard deviation of the restored image in each step. The standard deviation is calculated as
(50)
std
(
|
u
k
-
u
k
-
1
|
)
,
where u
k and u
k
-
1 represent the restored images of step k and step k
-
1, respectively. When calculating std
(
x
), we reshape matrix x into a vector x whose elements are taken column-wise from matrix x. Then
(51)
std
(
x
)
=
(
1
n
-
1
∑
i
=
1
n
(
x
i
-
x
-
)
2
)
1
/
2
.
Figure 6 shows that the alternating minimization algorithm converges faster.Figure 6
Comparison of the convergence rate.


## 6. Conclusion

In this paper, we propose a fast alternating minimization algorithm for multichannel image processing, based on the new nonlocal vectorial total variational model introduced recently. We use the variable splitting and penalty techniques to split the energy minimization problem into two optimization subproblems. One is a simplel
1 minimization problem, which is solved by the threshold method. The other is a simple l
2 minimization problem, which is solved by the variation calculus. As a result, this algorithm is faster than the classical algorithm because of avoiding dealing with the nondifferentiability of the nonlocal TV norm.The main contribution of this paper is introducing the alternating minimization algorithm for the nonlocal vectorial total variational model, which is proved fast and can be applied to many multichannel image processing problems. In this paper, we apply this algorithm to the RGB color image denoising as an example. It can be easily extended to other vectorial problems, such as color image deblurring, and to other color image models, such as hue-saturation-value (HSV) model, as well.



---

*Source: 731272-2014-08-26.xml*