# A Cartoon-Texture Decomposition Based Multiplicative Noise Removal Method

Authors: Chenping Zhao; Xiangchu Feng; Weiwei Wang; Huazhu Chen
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We propose a new frame for multiplicative noise removal. To improve the multiplicative denoising performance, we add the regularization of texture component in the denoising model, designing a multiscale multiplicative noise removal model. The proposed model is jointly convex and can be easily solved by optimization algorithms. We introduce Douglas-Rachford splitting method to solve the proposed model. In the algorithm, we make full use of some important proximity operators, which have closed expression or can be executed in one time iteration. In particular, the proximity ofH
-
1 norm is deduced, which is just the Fourier domain filtering. In the process of simulation experiments, we first analyze and select the needed parameters and then test the experiments on several images using the designed algorithm and the given parameters. Finally, we compare the denoising performance of the proposed model with the existing models, in which the signal to noise ratio (SNR) and the peak signal to noise ratios (PSNRs) are applied to evaluate the noise suppressing effects. Experimental results demonstrate that the designed algorithms can solve the model perfectly and the recovery images of the proposed model have higher SNRs/PSNRs and better visual quality.


---

## Body

## 1. Introduction

Image denoising is a basic and important task in image processing. The relatively mature developed denoising model is the additive noise model in which case the noise is assumed to obey a Gaussian distribution, that is,(1)
g
=
f
+
n
,
n
~
N
0
,
σ
2
.However, the noises involved in many applications do not conform to the characteristic of the additive one; they may corrupt an image in other forms. In this paper, we are concerned with the denoising problem under the assumption that the original image has been corrupted by some multiplicative noise. The corresponding examples include the uneven phenomena during the magnetic resonance imaging (MRI) and the speckle noise in the ultrasonic and in synthetic aperture radar (SAR) image. The degradation model of the multiplicative noise can be expressed as (2)
g
=
f
·
n
,where · refers to the componentwise multiplication. To model the actual problems, we assume that the multiplicative noise obeys some random distribution; for example, the noise in SAR images is assumed to obey the Gamma distribution and the one in ultrasonic image obeys the Rayleigh distribution.


### 1.1. Multiplicative Noise Removal Models

There are many multiplicative noise removal models; one of the most classic models based on TV regularization is called RLO model proposed by Rudin et al. [1]:(3)
m
i
n
f
∫
Ω
∇
f
d
x
s.t.
∫
Ω
g
x
f
x
d
x
=
1
,
∫
Ω
g
x
f
x
-
1
2
d
x
=
σ
2
,where the mean of noise is assumed to be 1 and the variance is assumed to be σ
2. Since (3) is a nonconvex problem, it is very difficult to be solved.The second classic model is called AA model [2] which is given by Aubert and Aujol under the hypothesis of Gamma distribution: (4)
m
i
n
f
∫
Ω
log
⁡
f
+
g
x
f
x
d
x
+
λ
∇
f
1
.The objective function in (4) is also nonconvex. However, Aubert and Aujol showed the existence of minimizers of the objective function and employed a gradient method to solve (4) numerically. There are several improved works based on (4) developed in recent years.Recently, Shi and Osher [3] proposed considering a logarithm transformation on the noisy observation, log
⁡
g
=
l
o
g
+
log
⁡
n, and derived the TV minimization model for multiplicative noise removal problems. Huang et al. also proposed the log-domain denoising model using the transformation f
=
exp
⁡
(
z
) in (4). Consider (5)
m
i
n
z
∫
Ω
z
x
+
g
x
e
-
z
x
d
x
+
λ
z
T
V
,which is known as EXP model [4]. The objective function in (5) is strictly convex in z and the authors promoted an alternating minimization algorithm to solve the model and showed the convergence of the method simultaneously. Nevertheless, model (5) is convex in z rather than in f in the original image domain. At the same time, Durand et al. [5] proposed a method composed of several stages. They also used the log-image data and applied reasonable suboptimal hard thresholding on its curvelet transform; then they applied a variational method by minimizing a specialized hybrid criterion composed of an l
1 data-fidelity term of the thresholded curvelet coefficients and a TV regularization term in the log-image domain. The restored image can be obtained by using an exponential of the minimizer, weighted in such a way that the mean of the original image is preserved. Their restored images combine the advantages of shrinkage and variational methods. Besides the above approaches, dictionary learning methods and nonlocal mean methods have also been proposed and developed for multiplicative denoising [6–9].In [10], Zhao et al. developed a convex optimization model for multiplicative noise removal. The main idea is to rewrite a multiplicative noise equation such that both the image variable and the noise variable are decoupled. That is to say, rewrite problem (2) as (6)
g
=
N
f
,where N
=
diag
⁡
(
n
) is a diagonal matrix of which the main diagonal entries are given by [
n
]
i. According to (2), when there is no noise in the observed image, we obtain that n
=
e, a vector of all ones. When there is a multiplicative noise in the observed image, we expect that {
[
n
]
i
0
≠
0
,
∀
i
}, and moreover they are greater than zeros. So we can say that N is invertible, and (6) is equivalent to (7)
G
w
=
f
,where G
=
diag
⁡
(
g
) is the diagonalizable matrix of vector g and [
w
]
i
=
1
/
[
n
]
i. The convex model to solve (2) is proposed as (8)
m
i
n
w
,
f
1
2
w
-
μ
e
2
2
+
α
1
G
w
-
f
1
+
α
2
f
T
V
.In (8), the first term is to measure the variance of w, the second term is the data term, and the third term is the TV regularization term. If the first term is absent and w can be arbitrarily assigned, then minimizing (9)
m
i
n
w
,
f
α
1
G
w
-
f
1
+
α
2
f
T
Vmay lead to the trivial solution f
≡
C
,
w
≡
C
G
-
1, where C is a constant. So the existence of variance term (
1
/
2
)
w
-
μ
e
2
2 can promise a proper minimizer of (8) [10]. Furthermore, It is worth noting that the image variable and the noise variable are decoupled, and the model is jointly convex for (
f
,
w
), which make the model able to be easily solved by optimizing methods.


### 1.2. The Contribution

We firstly analyze the performance of (8). It is observed that, if we fix w and set I
=
G
w, the minimizing problem (8) is just a TV-L1 problem, denoted as (10)
m
i
n
f
α
1
I
-
f
1
+
α
2
f
T
V
.In [11], it was pointed out that minimizing (10) can extract the large scale component and leave out the small scale contents of image I. In other words, we can extract cartoon component with different scales based on the selection of parameter λ
=
α
2
/
α
1. We verify the above conclusion through the following experiment. The example is shown in Figure 1(a), an image I with four squares with size of 3
×
3, 5
×
5, 20
×
20, and 80
×
80 pixels, respectively. The experiment results are the following: when λ
=
0.01, all the four squares can be extracted, as shown in Figure 1(b); when λ
=
0.05, three squares can be extracted while the one of 3
×
3 is lost, as shown in Figure 1(c); when λ
=
0.2, the two big squares appear, and another two are lost, as shown in Figure 1(d); when λ
=
2, only the biggest square 80
×
80 remains, while all other ones are lost, as shown in Figure 1(e). The above results are consistent with the theoretical analysis in [11]; the bigger λ is, the bigger the scale of extraction is. As a result, we have sound reasons to believe that the minimizer f of model (10) is almost the cartoon component of the whole restored image we expect.Figure 1
Numerical experiment results of (10) with different λ: (a) is the original image I; (b)–(e) are the minimizer f with λ
=
0.01
,
0.05
,
0.2
,
a
n
d
2, respectively.
(a)
I
(b)
λ
=
0.01
(c)
λ
=
0.05
(d)
λ
=
0.2
(e)
λ
=
2The above analysis shows that, if we want to recover more contents from an imageI, we should choose a small λ, which bring much more small scales and details. However, for a noisy image, the small scale component may contain much more noise, which may hinder the denoising quality. For this reason, we will improve the denoising model through adding the prior information of the texture component and design a cartoon-texture decomposition based denoising model. Then, we design the numerical algorithm based on the Douglas-Rachford splitting algorithm and three useful proximity operators. Before conducting the experiment, we analyze the selection of the corresponding parameters. Finally, we verify the performance of our proposed model through comparing the SNRs/PSNRs indexes with some existing models.The rest of the paper is organized as follows. We propose our new model in next section. In Section3, the numerical method is presented, in which the Douglas-Rachford (DR) splitting algorithm and some proximal operators used in this paper are firstly reviewed, and the algorithm is then designed and described in detail. In Section 4, the denoising experiments on two types of multiplicative noise are implemented, and the experiment results and performance analyses are unfolded. Finally, we conclude our work in Section 5.


## 1.1. Multiplicative Noise Removal Models

There are many multiplicative noise removal models; one of the most classic models based on TV regularization is called RLO model proposed by Rudin et al. [1]:(3)
m
i
n
f
∫
Ω
∇
f
d
x
s.t.
∫
Ω
g
x
f
x
d
x
=
1
,
∫
Ω
g
x
f
x
-
1
2
d
x
=
σ
2
,where the mean of noise is assumed to be 1 and the variance is assumed to be σ
2. Since (3) is a nonconvex problem, it is very difficult to be solved.The second classic model is called AA model [2] which is given by Aubert and Aujol under the hypothesis of Gamma distribution: (4)
m
i
n
f
∫
Ω
log
⁡
f
+
g
x
f
x
d
x
+
λ
∇
f
1
.The objective function in (4) is also nonconvex. However, Aubert and Aujol showed the existence of minimizers of the objective function and employed a gradient method to solve (4) numerically. There are several improved works based on (4) developed in recent years.Recently, Shi and Osher [3] proposed considering a logarithm transformation on the noisy observation, log
⁡
g
=
l
o
g
+
log
⁡
n, and derived the TV minimization model for multiplicative noise removal problems. Huang et al. also proposed the log-domain denoising model using the transformation f
=
exp
⁡
(
z
) in (4). Consider (5)
m
i
n
z
∫
Ω
z
x
+
g
x
e
-
z
x
d
x
+
λ
z
T
V
,which is known as EXP model [4]. The objective function in (5) is strictly convex in z and the authors promoted an alternating minimization algorithm to solve the model and showed the convergence of the method simultaneously. Nevertheless, model (5) is convex in z rather than in f in the original image domain. At the same time, Durand et al. [5] proposed a method composed of several stages. They also used the log-image data and applied reasonable suboptimal hard thresholding on its curvelet transform; then they applied a variational method by minimizing a specialized hybrid criterion composed of an l
1 data-fidelity term of the thresholded curvelet coefficients and a TV regularization term in the log-image domain. The restored image can be obtained by using an exponential of the minimizer, weighted in such a way that the mean of the original image is preserved. Their restored images combine the advantages of shrinkage and variational methods. Besides the above approaches, dictionary learning methods and nonlocal mean methods have also been proposed and developed for multiplicative denoising [6–9].In [10], Zhao et al. developed a convex optimization model for multiplicative noise removal. The main idea is to rewrite a multiplicative noise equation such that both the image variable and the noise variable are decoupled. That is to say, rewrite problem (2) as (6)
g
=
N
f
,where N
=
diag
⁡
(
n
) is a diagonal matrix of which the main diagonal entries are given by [
n
]
i. According to (2), when there is no noise in the observed image, we obtain that n
=
e, a vector of all ones. When there is a multiplicative noise in the observed image, we expect that {
[
n
]
i
0
≠
0
,
∀
i
}, and moreover they are greater than zeros. So we can say that N is invertible, and (6) is equivalent to (7)
G
w
=
f
,where G
=
diag
⁡
(
g
) is the diagonalizable matrix of vector g and [
w
]
i
=
1
/
[
n
]
i. The convex model to solve (2) is proposed as (8)
m
i
n
w
,
f
1
2
w
-
μ
e
2
2
+
α
1
G
w
-
f
1
+
α
2
f
T
V
.In (8), the first term is to measure the variance of w, the second term is the data term, and the third term is the TV regularization term. If the first term is absent and w can be arbitrarily assigned, then minimizing (9)
m
i
n
w
,
f
α
1
G
w
-
f
1
+
α
2
f
T
Vmay lead to the trivial solution f
≡
C
,
w
≡
C
G
-
1, where C is a constant. So the existence of variance term (
1
/
2
)
w
-
μ
e
2
2 can promise a proper minimizer of (8) [10]. Furthermore, It is worth noting that the image variable and the noise variable are decoupled, and the model is jointly convex for (
f
,
w
), which make the model able to be easily solved by optimizing methods.


## 1.2. The Contribution

We firstly analyze the performance of (8). It is observed that, if we fix w and set I
=
G
w, the minimizing problem (8) is just a TV-L1 problem, denoted as (10)
m
i
n
f
α
1
I
-
f
1
+
α
2
f
T
V
.In [11], it was pointed out that minimizing (10) can extract the large scale component and leave out the small scale contents of image I. In other words, we can extract cartoon component with different scales based on the selection of parameter λ
=
α
2
/
α
1. We verify the above conclusion through the following experiment. The example is shown in Figure 1(a), an image I with four squares with size of 3
×
3, 5
×
5, 20
×
20, and 80
×
80 pixels, respectively. The experiment results are the following: when λ
=
0.01, all the four squares can be extracted, as shown in Figure 1(b); when λ
=
0.05, three squares can be extracted while the one of 3
×
3 is lost, as shown in Figure 1(c); when λ
=
0.2, the two big squares appear, and another two are lost, as shown in Figure 1(d); when λ
=
2, only the biggest square 80
×
80 remains, while all other ones are lost, as shown in Figure 1(e). The above results are consistent with the theoretical analysis in [11]; the bigger λ is, the bigger the scale of extraction is. As a result, we have sound reasons to believe that the minimizer f of model (10) is almost the cartoon component of the whole restored image we expect.Figure 1
Numerical experiment results of (10) with different λ: (a) is the original image I; (b)–(e) are the minimizer f with λ
=
0.01
,
0.05
,
0.2
,
a
n
d
2, respectively.
(a)
I
(b)
λ
=
0.01
(c)
λ
=
0.05
(d)
λ
=
0.2
(e)
λ
=
2The above analysis shows that, if we want to recover more contents from an imageI, we should choose a small λ, which bring much more small scales and details. However, for a noisy image, the small scale component may contain much more noise, which may hinder the denoising quality. For this reason, we will improve the denoising model through adding the prior information of the texture component and design a cartoon-texture decomposition based denoising model. Then, we design the numerical algorithm based on the Douglas-Rachford splitting algorithm and three useful proximity operators. Before conducting the experiment, we analyze the selection of the corresponding parameters. Finally, we verify the performance of our proposed model through comparing the SNRs/PSNRs indexes with some existing models.The rest of the paper is organized as follows. We propose our new model in next section. In Section3, the numerical method is presented, in which the Douglas-Rachford (DR) splitting algorithm and some proximal operators used in this paper are firstly reviewed, and the algorithm is then designed and described in detail. In Section 4, the denoising experiments on two types of multiplicative noise are implemented, and the experiment results and performance analyses are unfolded. Finally, we conclude our work in Section 5.


## 2. The Proposed Model

Based on the analysis in Section1.2, the recovery image obtained by solving (8) mainly contains the cartoon component of what we expect, while most of the significant texture components are lost. To optimize the denoising performance, we improve model (8) by adding the texture information and alter the regularization term f
T
V to be u
T
V, where u denotes the cartoon component of the restored image f. Taking into consideration the prior information of texture v
=
f
-
u, we select v
H
-
1
2 as the regularization term about v, for that minimizing ·
H
-
1
2 can extract the texture component very well. In a word, the new model is (11)
w
∗
,
u
∗
,
v
∗
=
a
r
g
m
i
n
w
,
u
,
v
1
2
w
-
μ
e
2
2
+
α
1
G
w
-
u
+
v
1
+
α
2
u
T
V
+
ρ
v
H
-
1
2
,and the restored image is f
∗
=
u
∗
+
v
∗. In (11), α
1, α
2, and ρ are the positive regularization parameters to control the balance among the terms in the objective function.The proposed model possesses the following superiorities. Firstly, we add the termρ
v
H
-
1
2 in the regularization part, which can help to recover more image information and improve the quality of the restored image. Secondly, α
2
/
α
1 and ρ can be adjusted according to the noise level, which makes the model have multiscale property. Thirdly, when ρ
→
+
∞, then v
→
0, and (11) degrades to model (8); that is, (11) is the modified and improved version of (8). Finally, the proposed model keeps the structure of (8), is jointly convex for (
w
,
u
,
v
), and can be easily solved by optimal methods. In the following section, we use alternative iteration method and Douglas-Rachford method to deal with it, which combine the application of several proximity operators.


## 3. The Numerical Method

### 3.1. Douglas-Rachford Algorithm and Some Proximity Operators

Before proposing the numerical scheme of (11), we will review some basic algorithms and operators. The first one is the Douglas-Rachford splitting algorithm.Definition 1.
LetΦ
(
x
) and Ψ
(
x
) be the proper, l.s.c., and convex functions such that (12)
r
i
d
o
m
Ψ
∩
r
i
d
o
m
Φ
≠
∅
,and Ψ
(
x
)
+
Φ
(
x
)
→
+
∞ as x
→
+
∞; then the problem (13)
m
i
n
x
∈
R
N
Ψ
x
+
Φ
xadmits at least one solution and, for any γ
∈
(
0
,
+
∞
), its solutions are characterized by Algorithm 2.Algorithm 2 (Douglas-Rachford algorithm).
Consider the following:Fixingϵ
∈
(
0,1
)
,
γ
>
0
,
y
0
∈
R
N,
Forn
=
0,1
,
…
,
d
o
x
n
=
p
r
o
x
γ
Φ
y
n
λ
n
∈
[
ϵ
,
2
-
ϵ
]
y
n
+
1
=
y
n
+
λ
n
(
p
r
o
x
γ
Ψ
(
2
x
n
-
y
)
-
x
n
),
end forwhere the operator p
r
o
x
Φ will be defined in the following.Secondly, we review the definition of proximity operator and present several special examples that will be used in the numerical scheme.Definition 3.
LetΦ be a proper, l.s.c., and convex function; for every x
∈
R
N, the minimization problem (14)
m
i
n
u
Φ
u
+
1
2
u
-
x
2
2admits a unique solution, which is denoted by p
r
o
x
Φ
(
x
), and the operator p
r
o
x
Φ
(
x
)
:
R
N
→
R
N thus defined is called the proximity operator of Φ.In the following, we list three examples of proximity operator which will be used in this paper:(
1
)
Φ
(
x
)
=
λ
T
x
1
,
λ
=
(
λ
1
,
λ
2
,
…
,
λ
N
)
T
∈
R
N. Set y
=
p
r
o
x
γ
Φ
(
x
); then y can be componentwise expressed in the close form: (15)
y
i
=
sign
⁡
x
i
x
i
-
λ
i
γ
,
x
i
≥
λ
i
γ
0
,
e
l
s
e
i
=
1,2
,
…
,
N
,which is known as the soft thresholding operator, and we denote it in short by y
=
S
T
λ
γ
(
x
).(
2
)
C
o
n
s
i
d
e
r
Φ
(
x
)
=
x
T
V. Computing p
r
o
x
γ
Φ
(
x
) is equivalent to solving the TV-L2 problem (16)
m
i
n
u
1
2
u
-
x
2
2
+
γ
u
T
V
.There are many strategies to solve it. One of the most classical methods is the semi-implicit gradient descent algorithm proposed by Chambolle in [12]; another is the Forward-Backward method developed in [5]. Both of the above methods need inner loop in the processing of numerical scheme. In this paper, we adopt the variable splitting method and express the problem in discrete form: 16
′
m
i
n
u
1
2
u
-
x
2
2
+
γ
∇
u
1
,where ∇ is the gradient operator given by (17)
∇
=
∇
1
,
∇
2
,
∇
u
i
=
∇
1
u
i
,
∇
2
u
i
,
i
=
1,2
,
…
,
N
,
∇
u
1
=
∑
i
=
1
N
∇
1
u
i
2
+
∇
2
u
i
2
.We introduce an auxiliary variable p
∈
R
N
×
2 to transfer ∇
u: (18)
p
i
=
p
i
,
1
,
p
i
,
2
=
∇
u
i
=
∇
1
u
i
,
∇
2
u
i
,
i
=
1,2
,
…
,
N
;then we get the following approximation model to 16
′: (19)
m
i
n
u
,
p
1
2
u
-
x
2
2
+
γ
p
1
+
β
2
p
-
∇
u
2
2with a sufficiently large penalty parameter β. We solve (19) by an alternating minimization scheme given in Algorithm 4.Algorithm 4 (proximity of TV).
Consider the following:
(
1
) Fixing u, compute p componentwisemin
p
i
⁡
γ
p
i
2
+
(
β
/
2
)
p
i
-
∇
u
i
2
⇒
p
i
=
max
⁡
∇
u
i
2
-
γ
/
β
,
0
(
∇
u
i
/
∇
u
i
2
),
i
=
1
,
…
,
N
(
2
) Fixing p, compute u
min
u
⁡
(
1
/
2
)
u
-
x
2
2
+
(
1
/
β
)
p
-
∇
u
2
2
⇒
∇
T
∇
+
(
1
/
β
)
I
u
=
∇
T
p
+
(
1
/
β
)
x
⇒
u
=
∇
T
∇
+
(
1
/
β
)
I
-
1
∇
T
p
+
(
1
/
β
)
x
.It is worth mentioning that the above alternating minimization scheme is embedded into the whole algorithm with just one time iteration, which makes the algorithm have no inner iteration.(
3
)
C
o
n
s
i
d
e
r
Φ
(
x
)
=
x
H
-
1
2. We will deduce p
r
o
x
γ
Φ
(
x
) in the continuous situation, with the purpose of describing the derivation process clearly; this does not affect the processing of the numerical implementation. The proximity of Φ
(
x
) is expressed as (20)
p
r
o
x
γ
Φ
x
=
a
r
g
m
i
n
v
1
2
v
-
x
2
2
+
γ
v
H
-
1
2
.Denote the objective function in (20) by (21)
H
v
=
1
2
v
-
x
2
2
+
γ
v
H
-
1
2
, while minimizing H
(
v
) is equivalent to minimizing H
^
(
v
^
) in Fourier domain: (22)
H
^
v
^
=
1
2
v
^
-
x
^
2
2
+
γ
∫
v
^
ξ
2
2
π
ξ
2
d
ξ
; minimizing H
^
(
v
^
) in v
^ yields the unique solution v
^
=
L
^
x
^, where (23)
L
^
ξ
=
2
π
ξ
2
2
π
ξ
2
+
γ
. Taking inverse Fourier transform, we have (24)
v
∗
=
p
r
o
x
γ
1
Φ
1
z
=
L
∗
z
.It is indicated from (24) that the proximity of H
-
1 norm is just the Fourier domain filtering.


### 3.2. Solving the Proposed Model Numerically

We solve (11) by alternatively updating w
,
u
,
v in turn. The framework is given as in Algorithm 1.Algorithm 1
Initialization:u
1
=
g, v
1
=
0, f
1
=
u
1
+
v
1,
(1) for k
=
1,2
,
…,  do
(2)  w
k
+
1
=
arg
⁡
min
w
⁡
(
1
/
2
)
w
-
μ
e
2
2
+
α
1
G
w
-
f
k
1
(3)  u
k
+
1
=
arg
⁡
min
u
⁡
α
1
G
w
k
+
1
-
(
u
+
v
k
)
1
+
α
2
u
T
V
(4)  v
k
+
1
=
arg
⁡
min
v
⁡
α
1
G
w
k
+
1
-
(
u
k
+
1
+
v
)
1
+
α
2
ρ
v
H
-
1
2.
(5)  f
k
+
1
=
u
k
+
1
+
v
k
+
1.
(6) if stopping criterion is satisfied, then
(7)  stop the loop and output f
k
+
1.
(8) end forIn Algorithm1, updating w is to compute the proximity of ·
1, which is just the soft thresholding operator and can be realized by (15). Updating u is to solve a TV-L1 problem, and there are many methods to deal with it; to improve the algorithm efficiency, we realize the updating using Douglas-Rachford splitting algorithm, combining the computing of p
r
o
x
∥
·
∥
1 and p
r
o
x
∥
·
∥
T
V. Updating v is to solve an L1-H
-
1 problem, which is also dealt with by Douglas-Rachford splitting algorithm, combining the soft thresholding and Fourier domain filtering. A detailed description of the whole process is shown as in Algorithm 2.Algorithm 2
Initialization:u
(
1
)
=
g, v
(
1
)
=
0, f
=
u
+
v, y
=
u, z
=
0
(1) For k
=
1,2
,
…, do
(2) update w using soft-thresholding method
w
k
+
1
=
G
-
1
f
+
S
T
α
1
G
(
μ
e
-
G
-
1
f
k
)
(3) update u
y
k
+
1
=
y
k
+
λ
1
(
p
r
o
x
γ
1
α
1
⋅
-
(
G
w
k
+
1
-
v
k
)
1
(
2
u
k
-
y
k
)
-
u
k
)
u
k
+
1
=
p
r
o
x
γ
1
α
2
⋅
T
V
(
y
k
+
1
)
(4) update v
z
k
+
1
=
z
k
+
μ
1
(
prox
γ
2
α
1
⋅
-
(
G
w
k
+
1
-
u
k
+
1
)
1
(
2
v
k
-
z
k
+
1
)
-
v
k
)
v
k
+
1
=
p
r
o
x
γ
2
α
2
ρ
⋅
H
-
1
2
(
z
k
+
1
)
(5)    f
k
+
1
=
u
k
+
1
+
v
k
+
1
(6) end For, if stopping criterion is satisfied
(7) output f
=
u
+
v.There are some issues that need to be illustrated.Remark 5.
For step(
3
) in Algorithm 1, when we update u as u
k
+
1
=
p
r
o
x
γ
1
α
2
∥
·
∥
T
V
(
y
k
+
1
), we realize it by Algorithm 4 in one time iteration, which can avoid the inner loop and can improve the efficiency of the algorithm.Remark 6.
In this paper, the stopping criterion is used as follows:(25)
f
k
+
1
-
f
k
2
f
k
2
≤
1
e
-
4
.


## 3.1. Douglas-Rachford Algorithm and Some Proximity Operators

Before proposing the numerical scheme of (11), we will review some basic algorithms and operators. The first one is the Douglas-Rachford splitting algorithm.Definition 1.
LetΦ
(
x
) and Ψ
(
x
) be the proper, l.s.c., and convex functions such that (12)
r
i
d
o
m
Ψ
∩
r
i
d
o
m
Φ
≠
∅
,and Ψ
(
x
)
+
Φ
(
x
)
→
+
∞ as x
→
+
∞; then the problem (13)
m
i
n
x
∈
R
N
Ψ
x
+
Φ
xadmits at least one solution and, for any γ
∈
(
0
,
+
∞
), its solutions are characterized by Algorithm 2.Algorithm 2 (Douglas-Rachford algorithm).
Consider the following:Fixingϵ
∈
(
0,1
)
,
γ
>
0
,
y
0
∈
R
N,
Forn
=
0,1
,
…
,
d
o
x
n
=
p
r
o
x
γ
Φ
y
n
λ
n
∈
[
ϵ
,
2
-
ϵ
]
y
n
+
1
=
y
n
+
λ
n
(
p
r
o
x
γ
Ψ
(
2
x
n
-
y
)
-
x
n
),
end forwhere the operator p
r
o
x
Φ will be defined in the following.Secondly, we review the definition of proximity operator and present several special examples that will be used in the numerical scheme.Definition 3.
LetΦ be a proper, l.s.c., and convex function; for every x
∈
R
N, the minimization problem (14)
m
i
n
u
Φ
u
+
1
2
u
-
x
2
2admits a unique solution, which is denoted by p
r
o
x
Φ
(
x
), and the operator p
r
o
x
Φ
(
x
)
:
R
N
→
R
N thus defined is called the proximity operator of Φ.In the following, we list three examples of proximity operator which will be used in this paper:(
1
)
Φ
(
x
)
=
λ
T
x
1
,
λ
=
(
λ
1
,
λ
2
,
…
,
λ
N
)
T
∈
R
N. Set y
=
p
r
o
x
γ
Φ
(
x
); then y can be componentwise expressed in the close form: (15)
y
i
=
sign
⁡
x
i
x
i
-
λ
i
γ
,
x
i
≥
λ
i
γ
0
,
e
l
s
e
i
=
1,2
,
…
,
N
,which is known as the soft thresholding operator, and we denote it in short by y
=
S
T
λ
γ
(
x
).(
2
)
C
o
n
s
i
d
e
r
Φ
(
x
)
=
x
T
V. Computing p
r
o
x
γ
Φ
(
x
) is equivalent to solving the TV-L2 problem (16)
m
i
n
u
1
2
u
-
x
2
2
+
γ
u
T
V
.There are many strategies to solve it. One of the most classical methods is the semi-implicit gradient descent algorithm proposed by Chambolle in [12]; another is the Forward-Backward method developed in [5]. Both of the above methods need inner loop in the processing of numerical scheme. In this paper, we adopt the variable splitting method and express the problem in discrete form: 16
′
m
i
n
u
1
2
u
-
x
2
2
+
γ
∇
u
1
,where ∇ is the gradient operator given by (17)
∇
=
∇
1
,
∇
2
,
∇
u
i
=
∇
1
u
i
,
∇
2
u
i
,
i
=
1,2
,
…
,
N
,
∇
u
1
=
∑
i
=
1
N
∇
1
u
i
2
+
∇
2
u
i
2
.We introduce an auxiliary variable p
∈
R
N
×
2 to transfer ∇
u: (18)
p
i
=
p
i
,
1
,
p
i
,
2
=
∇
u
i
=
∇
1
u
i
,
∇
2
u
i
,
i
=
1,2
,
…
,
N
;then we get the following approximation model to 16
′: (19)
m
i
n
u
,
p
1
2
u
-
x
2
2
+
γ
p
1
+
β
2
p
-
∇
u
2
2with a sufficiently large penalty parameter β. We solve (19) by an alternating minimization scheme given in Algorithm 4.Algorithm 4 (proximity of TV).
Consider the following:
(
1
) Fixing u, compute p componentwisemin
p
i
⁡
γ
p
i
2
+
(
β
/
2
)
p
i
-
∇
u
i
2
⇒
p
i
=
max
⁡
∇
u
i
2
-
γ
/
β
,
0
(
∇
u
i
/
∇
u
i
2
),
i
=
1
,
…
,
N
(
2
) Fixing p, compute u
min
u
⁡
(
1
/
2
)
u
-
x
2
2
+
(
1
/
β
)
p
-
∇
u
2
2
⇒
∇
T
∇
+
(
1
/
β
)
I
u
=
∇
T
p
+
(
1
/
β
)
x
⇒
u
=
∇
T
∇
+
(
1
/
β
)
I
-
1
∇
T
p
+
(
1
/
β
)
x
.It is worth mentioning that the above alternating minimization scheme is embedded into the whole algorithm with just one time iteration, which makes the algorithm have no inner iteration.(
3
)
C
o
n
s
i
d
e
r
Φ
(
x
)
=
x
H
-
1
2. We will deduce p
r
o
x
γ
Φ
(
x
) in the continuous situation, with the purpose of describing the derivation process clearly; this does not affect the processing of the numerical implementation. The proximity of Φ
(
x
) is expressed as (20)
p
r
o
x
γ
Φ
x
=
a
r
g
m
i
n
v
1
2
v
-
x
2
2
+
γ
v
H
-
1
2
.Denote the objective function in (20) by (21)
H
v
=
1
2
v
-
x
2
2
+
γ
v
H
-
1
2
, while minimizing H
(
v
) is equivalent to minimizing H
^
(
v
^
) in Fourier domain: (22)
H
^
v
^
=
1
2
v
^
-
x
^
2
2
+
γ
∫
v
^
ξ
2
2
π
ξ
2
d
ξ
; minimizing H
^
(
v
^
) in v
^ yields the unique solution v
^
=
L
^
x
^, where (23)
L
^
ξ
=
2
π
ξ
2
2
π
ξ
2
+
γ
. Taking inverse Fourier transform, we have (24)
v
∗
=
p
r
o
x
γ
1
Φ
1
z
=
L
∗
z
.It is indicated from (24) that the proximity of H
-
1 norm is just the Fourier domain filtering.


## 3.2. Solving the Proposed Model Numerically

We solve (11) by alternatively updating w
,
u
,
v in turn. The framework is given as in Algorithm 1.Algorithm 1
Initialization:u
1
=
g, v
1
=
0, f
1
=
u
1
+
v
1,
(1) for k
=
1,2
,
…,  do
(2)  w
k
+
1
=
arg
⁡
min
w
⁡
(
1
/
2
)
w
-
μ
e
2
2
+
α
1
G
w
-
f
k
1
(3)  u
k
+
1
=
arg
⁡
min
u
⁡
α
1
G
w
k
+
1
-
(
u
+
v
k
)
1
+
α
2
u
T
V
(4)  v
k
+
1
=
arg
⁡
min
v
⁡
α
1
G
w
k
+
1
-
(
u
k
+
1
+
v
)
1
+
α
2
ρ
v
H
-
1
2.
(5)  f
k
+
1
=
u
k
+
1
+
v
k
+
1.
(6) if stopping criterion is satisfied, then
(7)  stop the loop and output f
k
+
1.
(8) end forIn Algorithm1, updating w is to compute the proximity of ·
1, which is just the soft thresholding operator and can be realized by (15). Updating u is to solve a TV-L1 problem, and there are many methods to deal with it; to improve the algorithm efficiency, we realize the updating using Douglas-Rachford splitting algorithm, combining the computing of p
r
o
x
∥
·
∥
1 and p
r
o
x
∥
·
∥
T
V. Updating v is to solve an L1-H
-
1 problem, which is also dealt with by Douglas-Rachford splitting algorithm, combining the soft thresholding and Fourier domain filtering. A detailed description of the whole process is shown as in Algorithm 2.Algorithm 2
Initialization:u
(
1
)
=
g, v
(
1
)
=
0, f
=
u
+
v, y
=
u, z
=
0
(1) For k
=
1,2
,
…, do
(2) update w using soft-thresholding method
w
k
+
1
=
G
-
1
f
+
S
T
α
1
G
(
μ
e
-
G
-
1
f
k
)
(3) update u
y
k
+
1
=
y
k
+
λ
1
(
p
r
o
x
γ
1
α
1
⋅
-
(
G
w
k
+
1
-
v
k
)
1
(
2
u
k
-
y
k
)
-
u
k
)
u
k
+
1
=
p
r
o
x
γ
1
α
2
⋅
T
V
(
y
k
+
1
)
(4) update v
z
k
+
1
=
z
k
+
μ
1
(
prox
γ
2
α
1
⋅
-
(
G
w
k
+
1
-
u
k
+
1
)
1
(
2
v
k
-
z
k
+
1
)
-
v
k
)
v
k
+
1
=
p
r
o
x
γ
2
α
2
ρ
⋅
H
-
1
2
(
z
k
+
1
)
(5)    f
k
+
1
=
u
k
+
1
+
v
k
+
1
(6) end For, if stopping criterion is satisfied
(7) output f
=
u
+
v.There are some issues that need to be illustrated.Remark 5.
For step(
3
) in Algorithm 1, when we update u as u
k
+
1
=
p
r
o
x
γ
1
α
2
∥
·
∥
T
V
(
y
k
+
1
), we realize it by Algorithm 4 in one time iteration, which can avoid the inner loop and can improve the efficiency of the algorithm.Remark 6.
In this paper, the stopping criterion is used as follows:(25)
f
k
+
1
-
f
k
2
f
k
2
≤
1
e
-
4
.


## 4. Implementation Details and Experimental Results

This section is mainly devoted to numerical simulation of image restoration in the presence of multiplicative noise. We test the performance of our model under the corruption of two types of multiplicative noise: Gamma and Rayleigh.Gamma Noise. The probability density function of Gamma noise n is given by (26)
p
n
=
n
k
-
1
Γ
k
θ
k
e
-
n
/
θ
,where θ
>
0
,
k
>
1. The mean and the variance of n are (27)
E
n
=
k
θ
,
V
n
=
k
θ
2
.As multiplicative noise of an image, the mean of n is set to be 1, so θ
=
1
/
k, and the variance of n is k
θ
2
=
1
/
k.In model (10), w
=
1
/
n
(
i
)
i
=
1
N, and the value of its mean is estimated as [10] (28)
μ
=
E
1
n
=
1
θ
k
-
1
=
k
k
-
1
.Rayleigh Noise. The probability density function of Rayleigh noise n is given by (29)
p
n
=
n
σ
2
exp
⁡
-
n
2
2
σ
2
,where σ is a positive parameter. The mean value of n is equal to π
/
2
σ, the variance of n is equal to (
2
-
π
/
2
)
/
σ
2, and the mean value of the inverse of Rayleigh noise can be estimated as [10] (30)
μ
=
E
1
n
=
π
2
σ
2
.We test experiments on three images: Cameraman (size of256
×
256), Lena (size of 512
×
512), and Aerial image of some city (size of 512
×
512); see Figure 2.Figure 2
The original images: (a) Cameraman, (b) Lena, and (c) Aerial image of some city.
(a)
(b)
(c)In this paper, we use the signal to noise ratio (SNR) and the peak signal to noise ratio (PSNR) as the evaluation indicators. These indicators are measured between clean and restored images. LetI
n
e
w represent the image restored from the noisy image I
n
o
i
s
y, and let σ
I
0
2 represent the average variance of the clean image I
0. Then, we define the SNR indicator by the formula as follows: (31)
S
N
R
=
10
log
10
⁡
d
σ
I
0
2
I
n
e
w
-
I
0
2
2
, and the PSNR is defined as (32)
P
S
N
R
=
10
log
10
⁡
d
I
0
∞
2
I
n
e
w
-
I
0
2
2
, where d denotes the length of the images.


### 4.1. Parameters Selection

Some necessary parametersα
1, α
2, ρ, λ
1, μ
1, γ
1, and γ
2 are needed to be given to start up the new algorithm; α
1, α
2, and ρ are regularization parameters to keep balance among the terms in model (11). In particular, λ
=
α
2
/
α
1 is an important indicator to trade off between the data term and the regularization term. λ
1 and μ
1 are two relaxation parameters of the Douglas-Rachford splitting algorithms in Algorithm 2 and are empirically set to be 0.5. γ
1 and γ
2 are the parameters of the proximity operator of ·
T
V and ·
H
-
1
2, respectively, and are empirically set to be 10.According to the experimental analysis in [10], the ratio λ
=
α
2
/
α
1 in (8) is almost a constant depending on the type of noise and was empirically estimated to be 5
/
6, 5
/
8 for Gamma noise and Rayleigh noise, respectively. From a number of experiments, we find that the best value of α
2
/
α
1 can be chosen from the range [
0.6,0.85
], and the exact selection depends on the different noise levels and images. We give the values of α
2
/
α
1 in Table 1.Table 1
PSNR and SNR results for the Gamma noise removal.
Image
k
λ
ρ
Method
PSNR
SNR
Cameraman
k
=
15
—
—
Noisy image
17.3027
5.1705
0.18
—
AA
24.9257
12.7825
0.8
—
Convex
25.5732
13.4701
0.8
20
Ours
25.7435
13.5874
Lena
k
=
15
—
—
Noisy image
17.0867
2.9056
0.18
—
AA
27.5628
11.5464
0.6
—
Convex
27.4350
11.0819
0.6
20
Ours
27.7890
11.6649
Aerial
k
=
15
—
—
Noisy image
24.2567
4.5260
0.18
—
AA
26.8085
6.0724
0.85
—
Convex
28.1713
6.2617
0.85
20
Ours
28.5635
7.5391
Cameraman
k
=
10
—
—
Noisy image
15.5091
3.3768
0.18
—
AA
24.1401
11.9376
0.8
—
Convex
24.3942
11.5807
0.8
50
Ours
24.7390
11.7593
Lena
k
=
10
—
—
Noisy image
15.2946
1.1407
0.18
—
AA
26.6863
12.3616
0.6
—
Convex
26.7026
12.4127
0.6
50
Ours
27.0777
12.3698
Aerial
k
=
10
—
—
Noisy image
22.5160
2.7469
0.18
—
AA
25.4918
5.7500
0.85
—
Convex
27.3225
6.1550
0.85
50
Ours
27.7026
7.2175
Cameraman
k
=
4
—
—
Noisy image
11.4822
−0.6501
0.18
—
AA
22.2638
9.8301
0.6
—
Convex
22.4161
10.6513
0.6
80
Ours
22.6186
10.4654
Lena
k
=
4
—
—
Noisy image
11.3537
−2.8313
0.18
—
AA
24.4908
9.7435
0.6
—
Convex
24.4912
9.0932
0.6
80
Ours
24.5983
10.2984
Aerial
k
=
4
—
—
Noisy image
18.5211
−1.2087
0.18
—
AA
24.4714
4.6684
0.85
—
Convex
24.5657
4.7386
0.85
80
Ours
25.0209
5.0016Forρ, it is a key parameter to promote the performance of model (11). Since a proper value of ρ is critical to the recovery of a satisfying image, it is necessary for us to emphatically discuss its impact on the denoising results. Here, we show its influences on image decomposition by using the image Cameraman. Fixing other parameters, we decompose the image into cartoon component and texture component with different ρ and put the experimental results in Figure 3. It is clear from Figure 3 that the different selections of ρ can bring completely different results. If ρ is selected too large (ρ
≥
100), there are very few contents in the texture v, and many details are found in u. With the decrease of ρ, more and more features get involved in the texture component v. However, when ρ is selected too little (e.g., ρ
≤
0.1), the block effect in cartoon component u is very serious, and some major features are included in the texture component v. Moreover, in the process of denoising, the too little value of ρ may bring much more details or small scale parts appearing in v, including some noise we do not expect. For this reason, we experimentally select ρ in the interval [
20,80
] according to the different noise level, and the exact selection can be found in Table 1.Figure 3
The influence ofρ on decomposition. (a) is image f. (b–e) are the cartoon u’s with different selection of ρ, which correspond to ρ
=
0.1,1
,
10, and 100, respectively. (f–i) are the texture v’s, which correspond to the u’s on (b–e).
(a)
Original image
(b)
Cartoon:ρ
=
100
(c)
Cartoon:ρ
=
10
(d)
Cartoon:ρ
=
1
(e)
Cartoon:ρ
=
0.1
(f)
Texture:ρ
=
100
(g)
Texture:ρ
=
10
(h)
Texture:ρ
=
1
(i)
Texture:ρ
=
0.1


### 4.2. Experimental Results in the Assumption of Gamma Noise

In the assumption of Gamma noise corruption, we test denoising performances of AA model, convex model, and the proposed model on the images of Figure2 and compare the SNRs/PSNRs of three methods. The AA model is formulated as (33)
m
i
n
f
∫
Ω
log
⁡
f
+
g
x
f
x
d
x
+
λ
∇
f
1
,which is solved by the gradient descent method (34)
f
i
+
1
-
f
i
Δ
t
=
div
⁡
Δ
f
i
Δ
f
i
+
λ
g
-
f
i
f
i
2
,
i
=
1,2
,
…
,where Δ
t is step size and is set as 0.01 for all experiments and λ is the regularization parameter and is adjusted according to the images and noise level; the exact value can be found in Table 1. The convex model is solved by the ADMM method; see [10].In the test, the original images are corrupted by Gamma noise withk
=
4, k
=
10, and k
=
15, respectively. In Table 1, we list values of related parameters and obtained SNRs/PSNRs by taking average of ten noise cases under the same experimental setting. It is clear that the proposed model performs quite better in terms of PSNR values and SNR values than AA model and convex model.We further show the denoised results obtained by three methods in Figures4–6. It is observed that when k
=
15 and k
=
10, as shown in Figures 4 and 5, the restored images from the new model possess more clear contents and much better visual effects than the AA model and the convex model. However, when k
=
4, that is, the noise level is higher, we can observe from Figure 6 that the new model possesses little advantage over the convex model in the visual effect. We also observe from Table 1 that, when k
=
4, the SNRs/PSNRs of our model and the convex model are neck to neck. The main reason for this phenomenon is that the larger the variance of the noise is, the more the texture components are seriously destroyed and regarded as noise part, and the bigger value ρ should be adjusted. As described in the analysis (c) of the new model, there will be a little content in the v parts, and denoising results of the new model will approximate that of the convex model.Figure 4
Experiment results under the corruption of Gamma noise,K
=
15.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
OursFigure 5
Experiment results under the corruption of Gamma noise,K
=
10.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
OursFigure 6
Experiment results under the corruption of Gamma noise,K
=
4.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
Ours


### 4.3. Experimental Results in the Assumption of Rayleigh Noise

The Rayleigh noise is generated byσ
-
2
log
⁡
(
1
-
b
), where b is a uniformly distributed random variable generated by “rand” in MATLAB. The mean of Rayleigh noise is set to be 1, and the variance is 0.2732.Since AA model is deduced under the assumption of Gamma noise, we just test the denoising performance of the convex model and the proposed model. Table2 shows the SNRs/PSNRs results of the two models. we notice that the proposed model gets higher SNRs and PSNRs than the convex model. The experiment results are shown as in Figure 7. It can be observed that the restored images from the new model possess more clear contents and much better visual effects.Table 2
PSNR and SNR results for the Rayleigh noise removal.
Image
Method
PSNR
SNR
Cameraman
Noisy image
11.1192
−1.013
Convex
22.4597
10.3275
Ours
22.6852
10.5529
Lena
Noisy image
10.9307
−3.2543
Convex
24.6532
10.4683
Ours
25.2379
11.0529
Aerial
Noisy image
18.1219
−1.6079
Convex
25.0408
5.3110
Ours
25.0688
5.3390Figure 7
Experiment results under the corruption of Rayleigh noise with the mean1 and the variance 0.2732.
(a)
Original
(b)
Noisy
(c)
Convex
(d)
Ours
(e)
Original
(f)
Noisy
(g)
Convex
(h)
Ours
(i)
Original
(j)
Noisy
(k)
Convex
(l)
Ours


### 4.4. Further Analysis of the Restored Images

In this section, we further analyze the cartoon partu and the texture part v of the restored image f of model (11). Take the case under Gamma noise as an example; we exhibit the experiment results when k
=
10.It can be observed from Figure8 that there are rich contents in the texture parts v, which enhance the image quality in both the SNRs/PSNRs index and the visual effects. So, the proposed model can also be regarded as a cartoon-texture decomposition method under the multiplicative noise corruption, on which we need further research.Figure 8
Experiment results of the proposed model. From (a–l), (a, e, i) are the noisy images corrupted by Gamma noise (K
=
10), (b, f, j) are the restored cartoon parts, (c, g, k) are the restored texture parts, and (d, h, l) are the whole recovery images.
(a)
Noisy
(b)
Cartoon
(c)
Texture
(d)
Ours
(e)
Noisy
(f)
Cartoon
(g)
Texture
(h)
Ours
(i)
Noisy
(j)
Cartoon
(k)
Texture
(l)
Ours


## 4.1. Parameters Selection

Some necessary parametersα
1, α
2, ρ, λ
1, μ
1, γ
1, and γ
2 are needed to be given to start up the new algorithm; α
1, α
2, and ρ are regularization parameters to keep balance among the terms in model (11). In particular, λ
=
α
2
/
α
1 is an important indicator to trade off between the data term and the regularization term. λ
1 and μ
1 are two relaxation parameters of the Douglas-Rachford splitting algorithms in Algorithm 2 and are empirically set to be 0.5. γ
1 and γ
2 are the parameters of the proximity operator of ·
T
V and ·
H
-
1
2, respectively, and are empirically set to be 10.According to the experimental analysis in [10], the ratio λ
=
α
2
/
α
1 in (8) is almost a constant depending on the type of noise and was empirically estimated to be 5
/
6, 5
/
8 for Gamma noise and Rayleigh noise, respectively. From a number of experiments, we find that the best value of α
2
/
α
1 can be chosen from the range [
0.6,0.85
], and the exact selection depends on the different noise levels and images. We give the values of α
2
/
α
1 in Table 1.Table 1
PSNR and SNR results for the Gamma noise removal.
Image
k
λ
ρ
Method
PSNR
SNR
Cameraman
k
=
15
—
—
Noisy image
17.3027
5.1705
0.18
—
AA
24.9257
12.7825
0.8
—
Convex
25.5732
13.4701
0.8
20
Ours
25.7435
13.5874
Lena
k
=
15
—
—
Noisy image
17.0867
2.9056
0.18
—
AA
27.5628
11.5464
0.6
—
Convex
27.4350
11.0819
0.6
20
Ours
27.7890
11.6649
Aerial
k
=
15
—
—
Noisy image
24.2567
4.5260
0.18
—
AA
26.8085
6.0724
0.85
—
Convex
28.1713
6.2617
0.85
20
Ours
28.5635
7.5391
Cameraman
k
=
10
—
—
Noisy image
15.5091
3.3768
0.18
—
AA
24.1401
11.9376
0.8
—
Convex
24.3942
11.5807
0.8
50
Ours
24.7390
11.7593
Lena
k
=
10
—
—
Noisy image
15.2946
1.1407
0.18
—
AA
26.6863
12.3616
0.6
—
Convex
26.7026
12.4127
0.6
50
Ours
27.0777
12.3698
Aerial
k
=
10
—
—
Noisy image
22.5160
2.7469
0.18
—
AA
25.4918
5.7500
0.85
—
Convex
27.3225
6.1550
0.85
50
Ours
27.7026
7.2175
Cameraman
k
=
4
—
—
Noisy image
11.4822
−0.6501
0.18
—
AA
22.2638
9.8301
0.6
—
Convex
22.4161
10.6513
0.6
80
Ours
22.6186
10.4654
Lena
k
=
4
—
—
Noisy image
11.3537
−2.8313
0.18
—
AA
24.4908
9.7435
0.6
—
Convex
24.4912
9.0932
0.6
80
Ours
24.5983
10.2984
Aerial
k
=
4
—
—
Noisy image
18.5211
−1.2087
0.18
—
AA
24.4714
4.6684
0.85
—
Convex
24.5657
4.7386
0.85
80
Ours
25.0209
5.0016Forρ, it is a key parameter to promote the performance of model (11). Since a proper value of ρ is critical to the recovery of a satisfying image, it is necessary for us to emphatically discuss its impact on the denoising results. Here, we show its influences on image decomposition by using the image Cameraman. Fixing other parameters, we decompose the image into cartoon component and texture component with different ρ and put the experimental results in Figure 3. It is clear from Figure 3 that the different selections of ρ can bring completely different results. If ρ is selected too large (ρ
≥
100), there are very few contents in the texture v, and many details are found in u. With the decrease of ρ, more and more features get involved in the texture component v. However, when ρ is selected too little (e.g., ρ
≤
0.1), the block effect in cartoon component u is very serious, and some major features are included in the texture component v. Moreover, in the process of denoising, the too little value of ρ may bring much more details or small scale parts appearing in v, including some noise we do not expect. For this reason, we experimentally select ρ in the interval [
20,80
] according to the different noise level, and the exact selection can be found in Table 1.Figure 3
The influence ofρ on decomposition. (a) is image f. (b–e) are the cartoon u’s with different selection of ρ, which correspond to ρ
=
0.1,1
,
10, and 100, respectively. (f–i) are the texture v’s, which correspond to the u’s on (b–e).
(a)
Original image
(b)
Cartoon:ρ
=
100
(c)
Cartoon:ρ
=
10
(d)
Cartoon:ρ
=
1
(e)
Cartoon:ρ
=
0.1
(f)
Texture:ρ
=
100
(g)
Texture:ρ
=
10
(h)
Texture:ρ
=
1
(i)
Texture:ρ
=
0.1


## 4.2. Experimental Results in the Assumption of Gamma Noise

In the assumption of Gamma noise corruption, we test denoising performances of AA model, convex model, and the proposed model on the images of Figure2 and compare the SNRs/PSNRs of three methods. The AA model is formulated as (33)
m
i
n
f
∫
Ω
log
⁡
f
+
g
x
f
x
d
x
+
λ
∇
f
1
,which is solved by the gradient descent method (34)
f
i
+
1
-
f
i
Δ
t
=
div
⁡
Δ
f
i
Δ
f
i
+
λ
g
-
f
i
f
i
2
,
i
=
1,2
,
…
,where Δ
t is step size and is set as 0.01 for all experiments and λ is the regularization parameter and is adjusted according to the images and noise level; the exact value can be found in Table 1. The convex model is solved by the ADMM method; see [10].In the test, the original images are corrupted by Gamma noise withk
=
4, k
=
10, and k
=
15, respectively. In Table 1, we list values of related parameters and obtained SNRs/PSNRs by taking average of ten noise cases under the same experimental setting. It is clear that the proposed model performs quite better in terms of PSNR values and SNR values than AA model and convex model.We further show the denoised results obtained by three methods in Figures4–6. It is observed that when k
=
15 and k
=
10, as shown in Figures 4 and 5, the restored images from the new model possess more clear contents and much better visual effects than the AA model and the convex model. However, when k
=
4, that is, the noise level is higher, we can observe from Figure 6 that the new model possesses little advantage over the convex model in the visual effect. We also observe from Table 1 that, when k
=
4, the SNRs/PSNRs of our model and the convex model are neck to neck. The main reason for this phenomenon is that the larger the variance of the noise is, the more the texture components are seriously destroyed and regarded as noise part, and the bigger value ρ should be adjusted. As described in the analysis (c) of the new model, there will be a little content in the v parts, and denoising results of the new model will approximate that of the convex model.Figure 4
Experiment results under the corruption of Gamma noise,K
=
15.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
OursFigure 5
Experiment results under the corruption of Gamma noise,K
=
10.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
OursFigure 6
Experiment results under the corruption of Gamma noise,K
=
4.
(a)
Noisy
(b)
AA
(c)
Convex
(d)
Ours
(e)
Noisy
(f)
AA
(g)
Convex
(h)
Ours
(i)
Noisy
(j)
AA
(k)
Convex
(l)
Ours


## 4.3. Experimental Results in the Assumption of Rayleigh Noise

The Rayleigh noise is generated byσ
-
2
log
⁡
(
1
-
b
), where b is a uniformly distributed random variable generated by “rand” in MATLAB. The mean of Rayleigh noise is set to be 1, and the variance is 0.2732.Since AA model is deduced under the assumption of Gamma noise, we just test the denoising performance of the convex model and the proposed model. Table2 shows the SNRs/PSNRs results of the two models. we notice that the proposed model gets higher SNRs and PSNRs than the convex model. The experiment results are shown as in Figure 7. It can be observed that the restored images from the new model possess more clear contents and much better visual effects.Table 2
PSNR and SNR results for the Rayleigh noise removal.
Image
Method
PSNR
SNR
Cameraman
Noisy image
11.1192
−1.013
Convex
22.4597
10.3275
Ours
22.6852
10.5529
Lena
Noisy image
10.9307
−3.2543
Convex
24.6532
10.4683
Ours
25.2379
11.0529
Aerial
Noisy image
18.1219
−1.6079
Convex
25.0408
5.3110
Ours
25.0688
5.3390Figure 7
Experiment results under the corruption of Rayleigh noise with the mean1 and the variance 0.2732.
(a)
Original
(b)
Noisy
(c)
Convex
(d)
Ours
(e)
Original
(f)
Noisy
(g)
Convex
(h)
Ours
(i)
Original
(j)
Noisy
(k)
Convex
(l)
Ours


## 4.4. Further Analysis of the Restored Images

In this section, we further analyze the cartoon partu and the texture part v of the restored image f of model (11). Take the case under Gamma noise as an example; we exhibit the experiment results when k
=
10.It can be observed from Figure8 that there are rich contents in the texture parts v, which enhance the image quality in both the SNRs/PSNRs index and the visual effects. So, the proposed model can also be regarded as a cartoon-texture decomposition method under the multiplicative noise corruption, on which we need further research.Figure 8
Experiment results of the proposed model. From (a–l), (a, e, i) are the noisy images corrupted by Gamma noise (K
=
10), (b, f, j) are the restored cartoon parts, (c, g, k) are the restored texture parts, and (d, h, l) are the whole recovery images.
(a)
Noisy
(b)
Cartoon
(c)
Texture
(d)
Ours
(e)
Noisy
(f)
Cartoon
(g)
Texture
(h)
Ours
(i)
Noisy
(j)
Cartoon
(k)
Texture
(l)
Ours


## 5. Conclusions

A novel cartoon-texture decomposition based multiplicative noise removal model has been proposed in this paper. The main advantage of our work is embodied in the following two aspects: the one in which we take most of the a priori information of texture component into the denoising model, which makes the model possess more flexibility and efficiency. Another is that we dexterously solve the new convex model using optimal splitting algorithms and Fourier domain filtering method. The experiment results demonstrate that the proposed model can better remove the noise than the two other models, and the recovery images have higher SNRs/PSNRs and better visual effects.



---

*Source: 5130346-2016-08-29.xml*