# Performance of Restarted Homotopy Perturbation Method for TV-Based Image Denoising Problem

Authors: Yu Du Han; Jae Heon Yun
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We first propose a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem which repeats HPM process by computing only the first few terms instead of computing infinite terms, and then we present an application of RHPM to TV- (Total Variation-) based image denoising problem. The main difficulty in applying RHPM to the nonlinear denoising problem is settled by using binomial series techniques. We also provide finite difference schemes for numerical implementation of RHPM. Lastly, numerical experiments for several test images are carried out to demonstrate the feasibility, efficiency, and reliability of RHPM by comparing the performance of RHPM with that of existing TM and recently proposed RHAM methods.


---

## Body

## 1. Introduction

Digital images generated by various digital devices or used in practical applications usually contain some noises, so denoising is necessary before analyzing the images [1]. Over the last decade, a number of image denoising techniques preserving important image information such as edges have been proposed to obtain the true image from the Noisy (known) image. These methods include filtering technique [2], wavelet techniques [3, 4], and variational methods [5–9].Denoising based on linear filters normally does not perform satisfactorily since both noise and edges contain high frequencies, while denoising based on the nonlinear models has been found to be successful. To preserve image edges, Rudin and Osher [8] and Rudin et al. [9] introduced a TV- (Total Variation-) based denoising model, which is based on a variational problem using the TV norm as a nonlinear nondifferentiable functional. Also Rudin et al. [9] used a nonlinear parabolic PDE (Partial Differential Equation) as a solution procedure for solving the Euler-Lagrange equation arising from the TV-based denoising model and proposed TM (Time Marching) method for solving the nonlinear parabolic PDE. This model is found to be a successful tool for image denoising and edge enhancement. High-order denoising models [10, 11] improving the TV denoising model have been proposed, but Euler-Lagrange equation associated with the TV functional for those models is a nonlinear PDE [12]. Due to the difficulties in obtaining analytic solutions of the nonlinear PDE, there has been many researches in developing implicit or explicit methods for the nonlinear PDE. Also, many good methods for image denoising problems such as high-order TV method [13, 14], total generalized variation (TGV) method, block-matching and 3D filtering (BM3D) method [15, 16], four-directional TV-based method [17], and an adaptive method based on Tikhonov and TV regularization [18] have been recently proposed.For the last 15 years, HPM (homotopy perturbation method), first introduced by He [19], has been used by many mathematicians and engineers to solve various functional equations. In this method, the solution is considered as the sum of an infinite series which converges rapidly to the exact solutions. Using homotopy technique in topology, a homotopy is constructed with an embedding parameter p
∈
[
0,1
] which is considered as a “small parameter.” Many researches have been conducted in applying the HPM to a class of linear and nonlinear equations. In addition, HPM was further applied to nonlinear oscillators with discontinuities [20], nonlinear wave equations [21], boundary value problem [22], and many other subjects [23–26]. Recently, Ghanbari et al. [27] proposed RHAM (restarted homotopy analysis method), based on HAM [28], for a nonlinear denoising model problem, which is a main motivation for this paper. The purpose of this paper is to propose a restarted homotopy perturbation method (RHPM), based on the HPM, for solving a nonlinear PDE problem and then to apply RHPM to the TV-based image denoising problem.This paper is organized as follows. In Section2, we review TM method for TV-based image denoising problem. In Section 3, we propose a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem. In Section 4, we present an application of RHPM to TV-based image denoising problem. The main difficulty in applying RHPM is settled by using binomial series techniques. In Section 5, numerical experiments for several test images are carried out to demonstrate the feasibility, efficiency, and reliability of RHPM by comparing the performance of RHPM with that of existing TM and recently proposed RHAM methods. Lastly, some conclusions are drawn.


## 2. Review of TM Method for TV-Based Image Denoising Problem

In this section, we only review TM method for solving a TV-based image denoising problem (see [27] for details of the recently proposed RHAM method). TM method has been proven to be a successful tool for image denoising and has been used extensively in many practical application.LetΩ be a bounded domain, let an intensity function z
(
x
,
y
) be the pixel values of an observed image, and let η
(
x
,
y
) be an unknown noise function. Then, we would like to construct the desired clean image u
(
x
,
y
) from the observed image z
(
x
,
y
) such that(1)
z
x
,
y
=
u
x
,
y
+
η
x
,
y
.Image denoising problem is to reconstruct u
(
x
,
y
) from z
(
x
,
y
) which contains the additive noise η
(
x
,
y
). It is well known that image denoising problem is an ill-conditioned problem, so regularization techniques should be used to approximate u
(
x
,
y
) from z
(
x
,
y
).Rudin et al. [9] proposed the use of the following TV-based minimization problem:(2)
m
i
n
u
∫
Ω
∇
u
x
,
y
d
x
d
y
+
λ
2
∫
Ω
u
x
,
y
-
z
x
,
y
2
d
x
d
y
,where λ
>
0 is a regularization parameter. The first term is the TV of u
(
x
,
y
), a regularizer, while the second term is a fidelity term ensuring that the denoised image u
(
x
,
y
) will be close to the given image z
(
x
,
y
). The regularization parameter λ is important for balancing denoising and smoothing. That is, it can be observed that the parameter λ may be small in the presence of high noise and λ may be large for little noise. The Euler-Lagrange equation corresponding to (2) leads to a nonlinear elliptic PDE with homogeneous Neumann boundary conditions as follows:(3)
g
u
x
,
y
=
∇
·
∇
u
x
,
y
∇
u
x
,
y
-
λ
u
x
,
y
-
z
x
,
y
=
0
,
x
,
y
∈
Ω
,
∇
u
x
,
y
·
n
→
=
0
,
x
,
y
∈
∂
Ω
,where n
→ is the unit normal vector exterior to the boundary ∂
Ω.To avoid division by zero in numerical implementation, we replace the nondifferentiable term∇
u
x
,
y by a smooth approximation as follows:(4)
∇
u
x
,
y
β
=
∇
u
x
,
y
2
+
β
2
for some small
β
>
0
.Hence, one obtains(5)
g
u
x
,
y
=
∇
·
∇
u
x
,
y
∇
u
x
,
y
β
-
λ
u
x
,
y
-
z
x
,
y
=
0
,
x
,
y
∈
Ω
,
∇
u
x
,
y
·
n
→
=
0
,
x
,
y
∈
∂
Ω
.In order to solve the Euler-Lagrange equation (5), Rudin et al. [9] proposed TM method for solving the following time-dependent nonlinear parabolic PDE(6)
∂
u
x
,
y
;
t
∂
t
=
∇
·
∇
u
x
,
y
;
t
∇
u
x
,
y
;
t
β
-
λ
u
x
,
y
;
t
-
z
x
,
y
,
t
>
0
,
x
,
y
∈
Ω
,
∇
u
·
n
→
=
0
,
x
,
y
∈
∂
Ω
,
where
u
x
,
y
,
0
is given
.For numerical implementation, let us assume that the domainΩ has been split into m
×
n cells, where the grid points are located at (
x
i
=
i
h
x
,
y
j
=
j
h
y
), 1
≤
i
≤
m, 1
≤
j
≤
n, t
k
=
k
Δ
t, where Δ
t and k
=
1,2
,
…, denote the time step and iteration number, respectively. We denote the values of u
(
x
,
y
;
t
) at the grid points (
x
i
,
y
j
;
t
k
) by u
i
j
k and u
i
j
0
=
z
(
x
i
,
y
j
).Without loss of generality, we can assume thath
x
=
h
y
=
1. Then, the nonlinear PDE (6) can be approximated by the following finite difference formula:(7)
u
i
j
k
+
1
=
u
i
j
k
+
Δ
t
·
κ
u
i
j
k
-
λ
u
i
j
k
-
z
i
j
,where(8)
κ
u
i
j
k
=
∇
·
∇
u
∇
u
β
i
j
k
=
∂
∂
x
u
x
∇
u
β
+
∂
∂
y
u
y
∇
u
β
i
j
k
=
Δ
-
x
Δ
+
x
u
i
j
k
Δ
+
x
u
i
j
k
2
+
Δ
+
y
u
i
j
k
2
+
β
2
+
Δ
-
y
Δ
+
y
u
i
j
k
Δ
+
x
u
i
j
k
2
+
Δ
+
y
u
i
j
k
2
+
β
2
i
jwith(9)
Δ
+
x
u
i
j
k
=
u
i
+
1
,
j
k
-
u
i
j
k
,
Δ
-
x
u
i
j
k
=
-
u
i
-
1
,
j
k
+
u
i
j
k
,
Δ
+
y
u
i
j
k
=
u
i
,
j
+
1
k
-
u
i
j
k
,
Δ
-
y
u
i
j
k
=
-
u
i
,
j
-
1
k
+
u
i
j
k
,and boundary conditions are(10)
u
0
j
k
=
u
1
j
k
,
u
m
+
1
,
j
k
=
u
m
j
k
,
u
i
0
k
=
u
i
1
k
,
u
i
,
n
+
1
k
=
u
i
,
n
k
,
for
1
≤
i
≤
m
,
1
≤
j
≤
n
.Hence, we obtain TM method, called Algorithm 1, for solving (6), where t
o
l is the tolerance for stopping criterion, m
a
x
i
t is the maximum number of iterations, z is the Noisy image, v is the true image, and PSNR is the peak signal-to-noise ratio between the true and restored images which is defined in Section 5.Algorithm 1:TM (Time-Marching) method.
(1)u
i
j
0
=
z
i
j for 1
≤
i
≤
m, 1
≤
j
≤
n
(2)for  k
=
0 to maxit  do
(3)Compute κ
u
i
j
k using  (8)  for 1
≤
i
≤
m, 1
≤
j
≤
n
(4)Compute u
i
j
k
+
1
=
u
i
j
k
+
Δ
t
⋅
κ
u
i
j
k
-
λ
u
i
j
k
-
z
i
j for 1
≤
i
≤
m, 1
≤
j
≤
n
(5)
if  ‖
u
k
-
u
k
+
1
‖
/
‖
u
k
‖
<
tol or PSNR
(
u
k
)
>
PSNR
(
u
k
+
1
)  then
(6) Stop
(7)
end if
(8)end for


## 3. Description of Restarted Homotopy Perturbation Method (RHPM)

In this section, we first describe the homotopy perturbation method (HPM) for solving nonlinear Partial Differential Equation(11)
A
u
-
f
r
=
0
,
r
∈
Ωwith the boundary conditions(12)
B
u
,
∂
u
∂
n
=
0
,
r
∈
Γ
,where A is a general differential operator, B is a boundary operator, f
(
r
) is a known analytic functions, and Γ is the boundary of the domain Ω. The operator A can be divided in two parts L and N, where L is linear and N is nonlinear. Therefore, (11) can be written as(13)
L
u
+
N
u
-
f
r
=
0
.By using homotopy technique, one can construct a homotopy v
(
r
,
p
)
:
Ω
×
[
0,1
]
→
R which satisfies(14)
H
v
,
p
=
1
-
p
L
v
-
L
u
0
+
p
A
v
-
f
r
=
0
,or, equivalently,(15)
H
v
,
p
=
L
v
-
L
u
0
+
p
L
u
0
+
p
N
v
-
f
r
=
0
,where p
∈
[
0,1
] is an embedding parameter and u
0 is the initial approximation of (11) which satisfies the boundary conditions. Clearly, we have(16)
H
v
,
0
=
L
v
-
L
u
0
=
0
,
H
v
,
1
=
A
v
-
f
r
=
0
.The changing process of p from zero to unity is just that of v
(
r
,
p
) changing from u
0
(
r
) to u
(
r
). This is called deformation and also L
(
v
)
-
L
(
u
0
) and A
(
v
)
-
f
(
r
) are called homotopic in topology. If the embedding parameter 0
≤
p
≤
1 is considered as a “small parameter” applying the classical perturbation technique, then the solution of (15) can be expressed as a power series in p; that is,(17)
v
=
v
r
,
p
=
v
0
+
p
v
1
+
p
2
v
2
+
⋯
.Setting p
→
1, the solution u of (11) can be expressed as the sum of an infinite series(18)
u
=
lim
p
→
1
⁡
v
r
,
p
=
v
0
+
v
1
+
v
2
+
⋯
.Further description and convergence for the HPM can be found in [19–22, 26].Before giving a general description of the restarted homotopy perturbation method (RHPM), we first provide a simple example of how RHPM can be applied to a PDE problem.Example 1.
Consider two-dimensional heat equation with variable coefficients of the form(19)
u
t
-
y
2
2
u
x
x
-
x
2
2
u
y
y
=
0with the initial condition u
(
x
,
y
,
0
)
=
y
2.Method 1 (HPM).
We first solve (19) using the traditional HPM. The operators L, N, and f in (13) can be viewed as(20)
L
u
=
u
t
,
N
u
=
-
y
2
2
u
x
x
-
x
2
2
u
y
y
,
f
r
=
0
.Let u
0
=
u
(
x
,
y
,
0
). From (15), we readily obtain a homotopy v
(
r
,
p
) which satisfies(21)
H
v
,
p
=
v
t
-
u
0
t
+
p
u
0
t
+
p
-
y
2
2
v
x
x
-
x
2
2
v
y
y
=
0
.Then, the solution v of (21) can be expressed as in (17). Substituting (17) into (21), one obtains(22)
v
0
t
+
p
v
1
t
+
p
2
v
2
t
+
⋯
=
u
0
t
-
p
u
0
t
+
p
y
2
2
v
0
x
x
+
p
v
1
x
x
+
⋯
+
x
2
2
v
0
y
y
+
p
v
1
y
y
+
⋯
=
u
0
t
+
p
-
u
0
t
+
y
2
2
v
0
x
x
+
x
2
2
v
0
y
y
+
p
2
y
2
2
v
1
x
x
+
x
2
2
v
1
y
y
+
⋯
.Comparing coefficients of terms with identical powers of p, we obtain that(23)
p
0
:
v
0
t
=
u
0
t
⟹
v
0
=
u
x
,
y
,
0
=
y
2
p
1
:
v
1
t
=
-
u
0
t
+
y
2
2
v
0
x
x
+
x
2
2
v
0
y
y
⟹
v
1
=
∫
0
t
-
v
0
t
+
1
2
y
2
v
0
x
x
+
x
2
v
0
y
y
d
t
p
2
n
:
v
2
n
t
=
y
2
2
v
2
n
-
1
x
x
+
x
2
2
v
2
n
-
1
y
y
n
≥
1
⟹
v
2
n
=
1
2
∫
0
t
y
2
v
2
n
-
1
x
x
+
x
2
v
2
n
-
1
y
y
d
t
p
2
n
+
1
:
v
2
n
+
1
t
=
y
2
2
v
2
n
x
x
+
x
2
2
v
2
n
y
y
n
≥
1
⟹
v
2
n
+
1
=
1
2
∫
0
t
y
2
v
2
n
x
x
+
x
2
v
2
n
y
y
d
t
.Hence, we have (24)
v
0
=
y
2
,
v
1
=
∫
0
t
x
2
d
t
=
t
x
2
v
2
n
=
∫
0
t
t
2
n
-
1
2
n
-
1
!
y
2
d
t
=
t
2
n
2
n
!
y
2
n
≥
1
v
2
n
+
1
=
∫
0
t
t
2
n
2
n
!
x
2
d
t
=
t
2
n
+
1
2
n
+
1
!
x
2
n
≥
1
.Therefore, we obtain the homotopy solution of (21) (25)
v
=
y
2
+
t
x
2
p
+
t
2
2
y
2
p
2
+
t
3
3
!
x
2
p
3
+
t
4
4
!
y
2
p
4
+
t
5
5
!
x
2
p
5
+
⋯
=
∑
n
=
0
∞
t
2
n
p
2
n
2
n
!
y
2
+
t
2
n
+
1
p
2
n
+
1
2
n
+
1
!
x
2
.Letting p
→
1, we obtain the exact solution of (19) which is given by (26)
u
=
lim
p
→
1
⁡
v
=
∑
n
=
0
∞
t
2
n
2
n
!
y
2
+
t
2
n
+
1
2
n
+
1
!
x
2
=
y
2
cosh
⁡
t
+
x
2
sinh
⁡
t
.Method 2 (RHPM).
We propose a new approach for solving (19), called RHPM, which is based on HPM. We choose u
0
=
u
(
x
,
y
,
0
)
=
y
2 as an initial approximation of (19). Let the operators L, N, and f be defined as in (20). Comparing coefficients of terms p
0 and p
1 in (22), one obtains(27)
p
0
:
v
0
t
=
u
0
t
⟹
v
0
=
u
0
=
u
x
,
y
,
0
=
y
2
p
1
:
v
1
t
=
-
u
0
t
+
y
2
2
v
0
x
x
+
x
2
2
v
0
y
y
⟹
v
1
=
∫
0
t
-
v
0
t
+
1
2
y
2
v
0
x
x
+
x
2
v
0
y
y
d
t
=
t
x
2
.Then, the first two term approximation of (18) is defined as follows: (28)
u
0
1
=
v
0
+
v
1
=
y
2
+
t
x
2
.From (27), using u
0
=
u
0
1 as an initial approximation, one easily obtains(29)
p
0
:
v
0
1
=
u
0
1
=
y
2
+
t
x
2
p
1
:
v
1
1
=
∫
0
t
-
v
0
1
t
+
1
2
y
2
v
0
1
x
x
+
x
2
v
0
1
y
y
d
t
=
t
2
2
y
2
.From (29), the new approximation u
0
2 is given by (30)
u
0
2
=
v
0
1
+
v
1
1
=
y
2
+
t
x
2
+
t
2
2
y
2
.From (27) with an initial approximation u
0
=
u
0
2, one obtains(31)
p
0
:
v
0
2
=
u
0
2
=
y
2
+
t
x
2
+
t
2
2
y
2
p
1
:
v
1
2
=
∫
0
t
-
v
0
2
t
+
1
2
y
2
v
0
2
x
x
+
x
2
v
0
2
y
y
d
t
=
t
3
3
!
x
2
.Hence, the new approximation u
0
3 is given by (32)
u
0
3
=
v
0
2
+
v
1
2
=
y
2
+
t
x
2
+
t
2
2
y
2
+
t
3
3
!
x
2
.Continuing this process, it can be seen that u
0
k converges to the exact solution u of (19), which is the same as the result of Method 1, as k
→
∞.Generally speaking, the computation ofv
k in the HPM becomes more complicated as k increases. To circumvent this problem, RHPM repeats the HPM process by computing only the first few terms instead of computing infinite terms of v
k. This is a big advantage of RHPM as compared with HPM. From the idea of RHPM introduced in Example 1, the general procedure of RHPM for solving a nonlinear equation is described below.Letu
0
0 be an initial approximation of a nonlinear equation. First, we compute the first kth term approximation u
0
1 in the following way: (33)
v
0
0
=
u
0
0
Compute
v
1
0
,
v
2
0
,
…
,
v
k
0
using  (
15
)  and  (
17
)
u
0
1
=
v
0
0
+
v
1
0
+
v
2
0
+
⋯
+
v
k
0
.Following the above computational process with u
0
1 as an initial approximation, a new approximation u
0
2 is computed as follows: (34)
v
0
1
=
u
0
1
Compute
v
1
1
,
v
2
1
,
…
,
v
k
1
using  (
15
)  and  (
17
)
u
0
2
=
v
0
1
+
v
1
1
+
v
2
1
+
⋯
+
v
k
1
.By repeating the above process, a new approximation u
0
l
+
1 can be computed from the initial approximation u
0
l, which was obtained at the lth step, in the following way:(35)
v
0
l
=
u
0
l
Compute
v
1
l
,
v
2
l
,
…
,
v
k
l
using  (
15
)  and  (
17
)
u
0
l
+
1
=
v
0
l
+
v
1
l
+
v
2
l
+
⋯
+
v
k
l
.Notice that v
0
l
,
v
1
l
,
…
,
v
k
l are coefficients of p
0
,
p
1
,
…
,
p
k, respectively. Then, the exact solution u can be approximated by u
0
l
+
1; that is, under some appropriate conditions, u
=
l
i
m
l
→
∞
u
0
l (see Example 1). In RHPM, which computes only k vectors v
1
l
,
v
2
l
,
…
,
v
k
l, every step is called R
H
P
M
(
k
) method. Clearly, RHPM introduced in Example 1 is RHPM(
1
).


## 4. Application of RHPM to TV-Based Image Denoising Problem

In this section, we describe an application of RHPM to a TV-based image denoising problem. Before doing this, we first describe the binomial series which is used for an application of RHPM. The difficulty in applying RHPM to the TV-based image denoising problem comes from (45) (i.e., (
v
x
2
+
v
y
2
)
-
1
/
2 and (
v
x
2
+
v
y
2
)
-
3
/
2). In order to handle this difficulty, the following binomial series is used. For x
<
1,(36)
1
+
x
n
=
∑
k
=
0
∞
n
k
x
k
=
1
+
n
x
+
n
n
-
1
2
!
x
2
+
n
n
-
1
n
-
2
3
!
x
3
+
⋯
.If n
=
-
1
/
2 in (36), then(37)
1
+
x
-
1
/
2
=
1
1
+
x
=
1
-
1
2
x
+
3
8
x
2
-
5
16
x
3
+
⋯
.Letting x
=
τ
-
1 in (37), for 0
<
τ
<
2,(38)
1
τ
=
1
-
1
2
τ
-
1
+
3
8
τ
-
1
2
-
5
16
τ
-
1
3
+
⋯
=
35
16
-
35
16
τ
+
21
16
τ
2
-
5
16
τ
3
+
⋯
.If τ
>
0 is large, then, by introducing a small parameter ɛ
>
0 and using (38), 1
/
τ can be approximated by(39)
1
τ
=
ɛ
ɛ
τ
≈
ɛ
35
16
-
35
ɛ
16
τ
+
21
ɛ
2
16
τ
2
-
5
ɛ
3
16
τ
3
.Similarly, if n
=
-
3
/
2 and x
=
τ
-
1 in (36), then, for x
<
1 or 0
<
τ
<
2,(40)
1
1
+
x
3
=
1
+
x
-
3
/
2
=
1
-
3
2
x
+
15
8
x
2
-
35
16
x
3
+
⋯
,
1
τ
3
=
95
16
-
179
16
τ
+
135
16
τ
2
-
35
16
τ
3
+
⋯
.If τ
>
0 is large, then, by introducing a small parameter ɛ
>
0 and using (40), 1
/
τ
3 can be approximated by(41)
1
τ
3
=
ɛ
3
ɛ
τ
3
≈
ɛ
3
95
16
-
179
ɛ
16
τ
+
135
ɛ
2
16
τ
2
-
35
ɛ
3
16
τ
3
.Note that a small parameter ɛ
>
0 is used to guarantee the convergence of the binomial series when τ
>
0 is large.Now, we describe how to apply RHPM to the TV-based image denoising problem (6). For simplicity of exposition, we describe computational process for RHPM(
2
) (i.e., k
=
2). Since we only use four terms of the binomial series for numerical computation of (45), division by zero does not occur and thus β
=
0 is used for RHPM (see (39) and (41)). Let the operators L, N, and f be defined as follows:(42)
L
v
=
∂
v
x
,
y
;
t
∂
t
,
N
v
=
-
∇
·
∇
v
x
,
y
;
t
∇
v
x
,
y
;
t
,
f
r
=
λ
v
x
,
y
;
t
-
z
x
,
y
,where (43)
∇
v
x
,
y
;
t
=
∂
v
∂
x
,
∂
v
∂
y
T
,
∇
v
x
,
y
;
t
=
v
x
2
+
v
y
2
.From (15), one obtains(44)
L
v
=
L
u
0
-
p
L
u
0
-
p
N
v
-
f
r
=
∂
u
0
∂
t
-
p
∂
u
0
∂
t
-
p
-
∇
·
∇
v
x
,
y
;
t
∇
v
x
,
y
;
t
-
λ
v
x
,
y
;
t
-
z
x
,
y
=
∂
u
0
∂
t
-
p
∂
u
0
∂
t
+
p
∇
·
∇
v
x
,
y
;
t
∇
v
x
,
y
;
t
+
p
λ
v
x
,
y
;
t
-
z
x
,
y
,where (45)
∇
·
∇
v
∇
v
=
v
x
x
+
v
y
y
v
x
2
+
v
y
2
-
v
x
2
v
x
x
+
v
y
2
v
y
y
v
x
2
+
v
y
2
3
.Since the solution v of (44) is of the form (17), by simple computation (46)
v
x
=
v
0
x
+
p
v
1
x
+
p
2
v
2
x
+
⋯
,
v
x
x
=
v
0
x
x
+
p
v
1
x
x
+
p
2
v
2
x
x
+
⋯
,
v
y
=
v
0
y
+
p
v
1
y
+
p
2
v
2
y
+
⋯
,
v
y
y
=
v
0
y
y
+
p
v
1
y
y
+
p
2
v
2
y
y
+
⋯
,
v
x
2
=
v
0
x
2
+
2
v
0
x
v
1
x
p
+
⋯
,
v
y
2
=
v
0
y
2
+
2
v
0
y
v
1
y
p
+
⋯
.From these equations, one easily obtains (47)
v
x
2
+
v
y
2
=
v
0
x
2
+
v
0
y
2
+
2
v
0
x
v
1
x
+
v
0
y
v
1
y
p
+
⋯
=
D
1
+
2
D
2
p
+
⋯
v
x
2
+
v
y
2
2
=
D
1
2
+
4
D
1
D
2
p
+
⋯
v
x
2
+
v
y
2
3
=
D
1
3
+
6
D
1
2
D
2
p
+
⋯
v
x
x
+
v
y
y
=
v
0
x
x
+
v
0
y
y
+
v
1
x
x
+
v
1
y
y
p
+
⋯
=
D
3
+
D
4
p
+
⋯
v
x
2
v
x
x
+
v
y
2
v
y
y
=
D
5
+
D
6
p
+
⋯
,where (48)
D
1
=
v
0
x
2
+
v
0
y
2
,
D
2
=
v
0
x
v
1
x
+
v
0
y
v
1
y
,
D
3
=
v
0
x
x
+
v
0
y
y
,
D
4
=
v
1
x
x
+
v
1
y
y
,
D
5
=
v
0
x
2
v
0
x
x
+
v
0
y
2
v
0
y
y
,
D
6
=
v
0
x
2
v
1
x
x
+
v
0
y
2
v
1
y
y
+
2
v
0
x
v
0
x
x
v
1
x
+
2
v
0
y
v
0
y
y
v
1
y
.The L
(
v
) and p
λ
v
-
z in (44) can be expressed as(49)
L
v
=
∂
v
∂
t
=
∂
v
0
∂
t
+
p
∂
v
1
∂
t
+
p
2
∂
v
2
∂
t
+
⋯
,
p
λ
v
-
z
=
λ
p
v
0
-
z
+
λ
p
2
v
1
+
λ
p
3
v
2
+
⋯
.Substituting τ
=
(
v
x
)
2
+
(
v
y
)
2 into (39), one obtains(50)
1
∇
v
=
ɛ
ɛ
τ
=
ɛ
ɛ
v
x
2
+
v
y
2
=
ɛ
35
16
-
35
ɛ
16
v
x
2
+
v
y
2
+
21
ɛ
2
16
v
x
2
+
v
y
2
2
-
5
ɛ
3
16
v
x
2
+
v
y
2
3
+
⋯
=
35
ɛ
16
-
35
ɛ
ɛ
16
v
0
x
2
+
v
0
y
2
+
2
v
0
x
v
1
x
+
v
0
y
v
1
y
p
+
⋯
+
21
ɛ
2
ɛ
16
v
0
x
2
+
v
0
y
2
2
+
4
v
0
x
2
+
v
0
y
2
v
0
x
v
1
x
+
v
0
y
v
1
y
p
+
⋯
-
5
ɛ
3
ɛ
16
v
0
x
2
+
v
0
y
2
3
+
6
v
0
x
2
+
v
0
y
2
2
v
0
x
v
1
x
+
v
0
y
v
1
y
p
+
⋯
+
⋯
=
35
ɛ
16
-
35
ɛ
ɛ
16
D
1
+
2
D
2
p
+
⋯
+
21
ɛ
2
ɛ
16
D
1
2
+
4
D
1
D
2
p
+
⋯
-
5
ɛ
3
ɛ
16
D
1
3
+
6
D
1
2
D
2
p
+
⋯
+
⋯
=
A
1
+
B
1
D
2
p
+
⋯
,where (51)
A
1
=
ɛ
35
16
-
35
ɛ
16
D
1
+
21
ɛ
2
16
D
1
2
-
5
ɛ
3
16
D
1
3
B
1
=
ɛ
ɛ
-
35
8
+
21
ɛ
4
D
1
-
15
ɛ
2
8
D
1
2
.Hence, the first term of the right hand side in (45) can be written as(52)
v
x
x
+
v
y
y
∇
v
=
A
1
D
3
+
A
1
D
4
p
+
B
1
D
2
D
3
p
+
⋯
.Substituting τ
=
(
v
x
)
2
+
(
v
y
)
2 into (41), one obtains (53)
1
∇
v
3
=
ɛ
3
ɛ
τ
3
=
ɛ
3
ɛ
3
v
x
2
+
v
y
2
3
=
ɛ
3
95
16
-
179
ɛ
16
v
x
2
+
v
y
2
+
135
ɛ
2
16
v
x
2
+
v
y
2
2
-
35
ɛ
3
16
v
x
2
+
v
y
2
3
+
⋯
=
95
ɛ
3
16
-
179
ɛ
ɛ
3
16
D
1
+
2
D
2
p
+
⋯
+
135
ɛ
2
ɛ
3
16
D
1
2
+
4
D
1
D
2
p
+
⋯
-
35
ɛ
3
ɛ
3
16
D
1
3
+
6
D
1
2
D
2
p
+
⋯
+
⋯
=
A
2
+
B
2
D
2
p
+
⋯
,where (54)
A
2
=
ɛ
3
95
16
-
179
ɛ
16
D
1
+
135
ɛ
2
16
D
1
2
-
35
ɛ
3
16
D
1
3
,
B
2
=
ɛ
ɛ
3
-
179
8
+
135
ɛ
4
D
1
-
105
ɛ
2
8
D
1
2
.Hence, the second term of the right hand side in (45) can be written as(55)
v
x
2
v
x
x
+
v
y
2
v
y
y
∇
v
3
=
A
2
D
5
+
B
2
D
2
D
5
p
+
A
2
D
6
p
+
⋯
.From (45), (52), and (55), (44) can be expressed as(56)
∂
v
0
∂
t
+
p
∂
v
1
∂
t
+
p
2
∂
v
2
∂
t
+
⋯
=
∂
u
0
∂
t
-
p
∂
u
0
∂
t
+
p
v
x
x
+
v
y
y
∇
v
-
v
x
2
v
x
x
+
v
x
2
v
y
y
∇
v
3
+
p
λ
v
0
-
z
+
p
2
λ
v
1
+
⋯
=
∂
u
0
∂
t
-
p
∂
u
0
∂
t
-
λ
v
0
-
z
+
p
2
λ
v
1
+
⋯
+
A
1
D
3
p
-
A
2
D
5
p
+
A
1
D
4
p
2
+
B
1
D
2
D
3
p
2
-
B
2
D
2
D
5
p
2
-
A
2
D
6
p
2
+
⋯
.Comparing the coefficients of p
0, p
1, and p
2 in (56),(57)
p
0
:
∂
v
0
∂
t
=
∂
u
0
∂
t
⟹
v
0
=
u
0
=
z
:  given an initial approximation
(58)
p
1
:
∂
v
1
∂
t
=
-
∂
v
0
∂
t
+
λ
v
0
-
z
+
A
1
D
3
-
A
2
D
5
(59)
p
2
:
∂
v
2
∂
t
=
λ
v
1
+
A
1
D
4
+
B
1
D
3
-
B
2
D
5
D
2
-
A
2
D
6
.Applying the inverse operator of L to (58), one obtains(60)
v
1
=
-
v
0
-
z
+
λ
v
0
-
z
t
+
A
1
D
3
t
-
A
2
D
5
t
=
-
F
+
λ
F
+
G
t
,where F
=
v
0
-
z and G
=
(
A
1
)
(
D
3
)
-
(
A
2
)
(
D
5
). From (60), the partial derivatives of v
1 are(61)
v
1
x
=
-
F
x
+
t
λ
F
+
G
x
=
-
F
x
+
λ
F
x
+
G
x
t
=
-
F
x
+
H
1
t
,
v
1
y
=
-
F
x
+
t
λ
F
+
G
y
=
-
F
x
+
λ
F
y
+
G
y
t
=
-
F
y
+
H
2
t
,
v
1
x
x
=
-
F
x
x
+
t
λ
F
+
G
x
x
=
-
F
x
x
+
λ
F
x
x
+
G
x
x
t
=
-
F
x
x
+
H
3
t
,
v
1
y
y
=
-
F
y
y
+
t
λ
F
+
G
y
y
=
-
F
y
y
+
λ
F
y
y
+
G
y
y
t
=
-
F
y
y
+
H
4
t
,where (
H
1
)
=
λ
F
x
+
G
x
,
(
H
2
)
=
λ
F
y
+
G
y, (
H
3
)
=
λ
F
x
x
+
G
x
x, and (
H
4
)
=
λ
F
y
y
+
G
y
y. From (61), (
D
2
), (
D
4
), and (
D
6
) can be expressed as(62)
D
2
=
v
0
x
-
F
x
+
H
1
t
+
v
0
y
-
F
y
+
H
2
t
,
D
4
=
-
F
x
x
+
F
y
y
+
H
3
+
H
4
t
,
D
6
=
v
0
x
2
-
F
x
x
+
H
3
t
+
v
0
y
2
-
F
y
y
+
H
4
t
+
2
v
0
x
v
0
x
x
-
F
x
+
H
1
t
+
2
v
0
y
v
0
y
y
-
F
y
+
H
2
t
=
-
v
0
x
2
F
x
x
+
v
0
y
2
F
y
y
+
2
v
0
x
v
0
x
x
F
x
+
v
0
y
v
0
y
y
F
y
+
v
0
x
2
H
3
+
v
0
y
2
H
4
+
2
v
0
x
v
0
x
x
H
1
+
v
0
y
v
0
y
y
H
2
t
.Substituting (60) and (62) into (59) and applying the inverse operator of L to (59), one obtains(63)
v
2
=
-
λ
F
+
A
1
F
x
x
+
F
y
y
+
B
1
D
3
-
B
2
D
5
v
0
x
F
x
+
v
0
y
F
y
t
+
A
2
v
0
x
2
F
x
x
+
v
0
y
2
F
y
y
+
2
v
0
x
v
0
x
x
F
x
+
v
0
y
v
0
y
y
F
y
t
+
λ
λ
F
+
G
1
2
t
2
+
A
1
H
3
+
H
4
1
2
t
2
+
B
1
D
3
-
B
2
D
5
v
0
x
H
1
+
v
0
y
H
2
1
2
t
2
-
A
2
v
0
x
2
H
3
+
v
0
y
2
H
4
+
2
v
0
x
v
0
x
x
H
1
+
2
v
0
y
v
0
y
y
H
2
1
2
t
2
.Now, update u
0 whose new value is set to v
0
+
v
1
+
v
2. By repeating the aforementioned process with the new updated u
0, we obtain RHPM(
2
) method.For numerical implementation of RHPM(
2
) method, let us assume that the domain Ω has been split into m
×
n cells, where the grid points are located at (
x
i
=
i
h
x, y
j
=
j
h
y
), 1
≤
i
≤
m, 1
≤
j
≤
n. Let t
k
=
k
Δ
t, where Δ
t and k
=
1,2
,
…, refer to the time step and iteration number, respectively. We denote the values of u
(
x
,
y
;
t
) at the grid points (
x
i
,
y
j
;
t
k
) by (
u
k
)
i
j. Without loss of generality, we can assume that h
=
h
x
=
h
y
=
1. For simplicity of exposition, we define the following notation: (64)
v
0
k
=
v
0
k
i
j
=
v
0
x
i
,
y
j
;
t
k
∣
1
≤
i
≤
m
,
1
≤
j
≤
n
,
k
=
1,2
,
…
,
v
1
k
=
v
1
k
i
j
=
v
1
x
i
,
y
j
;
t
k
∣
1
≤
i
≤
m
,
1
≤
j
≤
n
,
k
=
1,2
,
…
,
v
2
k
=
v
2
k
i
j
=
v
2
x
i
,
y
j
;
t
k
∣
1
≤
i
≤
m
,
1
≤
j
≤
n
,
k
=
1,2
,
…
.Boundary conditions are defined the same as those in TM method. Assume that (
u
0
0
)
i
j
=
z
(
x
i
,
y
j
)
=
z
i
j and (
u
0
k
)
i
j has been computed for all i and j. Then, (
v
0
k
)
i
j
=
(
u
0
k
)
i
j and discretization of (60) is given by(65)
v
1
k
i
j
=
-
F
i
j
+
λ
F
i
j
+
G
i
j
t
k
.Also, discretization of (63) is given by (66)
v
2
k
i
j
=
-
λ
F
i
j
+
A
1
i
j
F
x
x
+
F
y
y
+
B
1
i
j
D
3
i
j
-
B
2
i
j
D
5
i
j
v
0
x
F
x
+
v
0
y
F
y
t
k
+
A
2
i
j
v
0
x
2
F
x
x
+
v
0
y
2
F
y
y
+
2
v
0
x
v
0
x
x
F
x
+
v
0
y
v
0
y
y
F
y
t
k
+
λ
2
F
i
j
+
λ
G
i
j
+
A
1
i
j
H
3
i
j
+
H
4
i
j
1
2
t
k
2
+
B
1
i
j
D
3
i
j
-
B
2
i
j
D
5
i
j
v
0
x
H
1
i
j
+
v
0
y
H
2
i
j
1
2
t
k
2
-
A
2
i
j
v
0
x
2
H
3
i
j
+
2
v
0
x
v
0
x
x
H
1
i
j
+
v
0
y
2
H
4
i
j
+
2
v
0
y
v
0
y
y
H
2
i
j
1
2
t
k
2
,where (67)
v
0
i
j
=
v
0
k
i
j
=
v
0
k
-
1
i
j
+
v
1
k
-
1
i
j
+
v
2
k
-
1
i
j
,
v
0
x
=
v
0
x
i
j
=
v
0
i
+
1
j
-
v
0
i
j
,
v
0
x
x
=
v
0
x
x
i
j
=
v
0
i
-
1
j
-
2
v
0
i
j
+
v
0
i
+
1
j
,
v
0
y
=
v
0
y
i
j
=
v
0
i
j
+
1
-
v
0
i
j
,
v
0
y
y
=
v
0
y
y
i
j
=
v
0
i
j
-
1
-
2
v
0
i
j
+
v
0
i
j
+
1
,
F
i
j
=
v
0
i
j
-
z
i
,
j
,
G
i
j
=
A
1
i
j
D
3
i
j
-
A
2
i
j
D
5
i
j
,
F
x
=
F
x
i
j
=
F
i
+
1
j
-
F
i
j
,
F
x
x
=
F
x
x
i
j
=
F
i
-
1
j
-
2
F
i
j
+
F
i
+
1
j
,
F
y
=
F
y
i
j
=
F
i
j
+
1
-
F
i
j
,
F
y
y
=
F
y
y
i
j
=
F
i
j
-
1
-
2
F
i
j
+
F
i
j
+
1
,
D
1
i
j
=
v
0
i
+
1
j
-
v
0
i
j
2
+
v
0
i
j
+
1
-
v
0
i
j
2
,
D
3
i
j
=
v
0
i
-
1
j
-
2
v
0
i
j
+
v
0
i
+
1
j
+
v
0
i
j
-
1
-
2
v
0
i
j
+
v
0
i
j
+
1
,
D
5
i
j
=
v
0
i
+
1
j
-
v
0
i
j
2
v
0
i
-
1
j
-
2
v
0
i
j
+
v
0
i
+
1
j
+
v
0
i
j
+
1
-
v
0
i
j
2
v
0
i
j
-
1
-
2
v
0
i
j
+
v
0
i
j
+
1
,
(68)
A
1
i
j
=
ɛ
35
16
-
35
ɛ
16
D
1
i
j
+
21
ɛ
2
16
D
1
i
j
2
-
5
ɛ
3
16
D
1
i
j
3
,
A
2
i
j
=
ɛ
3
95
16
-
179
ɛ
16
D
1
i
j
+
135
ɛ
2
16
D
1
i
j
2
-
35
ɛ
3
16
D
1
i
j
3
,
B
1
i
j
=
ɛ
3
-
35
8
+
21
ɛ
4
D
1
i
j
-
15
ɛ
2
8
D
1
i
j
2
,
B
2
i
j
=
ɛ
5
-
179
8
+
135
ɛ
4
D
1
i
j
-
105
ɛ
2
8
D
1
i
j
2
,
(69)
H
1
i
j
=
λ
F
i
+
1
j
-
F
i
j
+
G
i
+
1
j
-
G
i
j
,
H
2
i
j
=
λ
F
i
j
+
1
-
F
i
j
+
G
i
j
+
1
-
G
i
j
,
H
3
i
j
=
λ
F
i
-
1
j
-
2
F
i
j
+
F
i
+
1
j
+
G
i
-
1
j
-
2
G
i
j
+
G
i
+
1
j
,
H
4
i
j
=
λ
F
i
j
-
1
-
2
F
i
j
+
F
i
j
+
1
+
G
i
j
-
1
-
2
G
i
j
+
G
i
j
+
1
.We now compute u
0
k
+
1
=
v
0
k
+
v
1
k
+
v
2
k. Repeating the above process with u
0
k
+
1 as an initial approximation, we can obtain RHPM(
2
) method, called Algorithm 2, for solving (6).Algorithm 2:RHPM(2) method.
(1)(
u
0
0
)
i
j
=
z
i
j for all i and j
(2)for  k
=
0 to maxit do
(3)Set v
0
k
=
u
0
k
(4)Compute v
1
k using (65) for all i and j
(5)Compute v
2
k using (66) for all i and j
(6)
u
0
k
+
1
=
v
0
k
+
v
1
k
+
v
2
k
(7)
if  ‖
u
0
k
-
u
0
k
+
1
‖
/
‖
u
0
k
‖
<
tol or PSNR
(
u
0
k
)
>
PSNR
(
u
0
k
+
1
)  then
(8) Stop
(9)
end if
(10)end forSince the computation ofv
2
k in RHPM(
2
) method is so complicated, we also provide RHPM(
1
) method, called Algorithm 3, whose computational cost is much cheaper than Algorithm 2. In RHPM(
2
) and RHPM(
1
) methods, t
o
l, m
a
x
i
t, z, and u are the tolerance for stopping criterion, the maximum number of iterations, the Noisy image, and the true image, respectively. In (68), the coefficients of (
D
1
) may have large values which cause amplification of some incorrect results. Hence, we have used an intensity range [
0,1
] for all test images used in RHPM methods to guarantee numerical stability.Algorithm 3:RHPM(1) method.
(1)(
u
0
0
)
i
j
=
z
i
j for all i and j
(2)for  k
=
0 to maxit do
(3)Set v
0
k
=
u
0
k
(4)Compute v
1
k using (65) for all i and j
(5)
u
0
k
+
1
=
v
0
k
+
v
1
k
(6)
if  ‖
u
0
k
-
u
0
k
+
1
‖
/
‖
u
0
k
‖
<
tol or PSNR
(
u
0
k
)
>
PSNR
(
u
0
k
+
1
)  then
(7) Stop
(8)
end if
(9)end for


## 5. Numerical Experiments

In this section, numerical experiments for several test images are carried out to compare the performance of RHPM(
2
) and RHPM(
1
) with those of existing TM and recently proposed RHAM2 in [27]. Numerical tests are performed using MATLAB R2014a on a personal computer with a processer at 3.6 GHz and 8 GB RAM. The Noisy images are generated by adding Gaussian white noise to the clean images using the built-in MATLAB function r
a
n
d
n. For example, adding 15
% of Gaussian white noise to the true image B can be done as follows: (70)
E
=
r
a
n
d
n
m
,
n
;
E
=
E
E
F
;
B
~
=
B
+
0.15
×
B
F
×
E
,where ·
F refers to the Frobenius norm.In order to illustrate efficiency and reliability of the proposed denoising algorithms (i.e., Algorithms2 and 3), we provide experimental results for several images such as artificial, synthetic, and real images. To evaluate the quality of the restored images, we use the peak signal-to-noise ratio (PSNR) between the restored image and Original image which is defined by (71)
P
S
N
R
=
10
log
10
⁡
m
·
n
·
max
⁡
i
,
j
⁡
u
i
j
2
u
-
u
k
F
2
,where u and u
k are the Original and restored images, respectively. Also u
i
j stands for the value of Original image u at the grid point (
i
,
j
) and m
·
n is the total number of pixels. It is clear that the larger PSNR is, the better the quality denoising image is.We have tested the denoising algorithms for several images with an intensity range[
0,255
] for TM and RHAM2 methods and [
0,1
] for RHPM(
2
) and RHPM(
1
) methods, and Noisy images are generated by adding 15
% or 20
% of Gaussian white noise to true images. For all of the 4 algorithms, parameters t
o
l
=
1
0
-
4 and m
a
x
i
t
=
5000 have been used. In TM and RHAM2 methods, parameters λ
=
1
0
-
6 and β
=
1
0
-
3 are used. In RHPM(
2
) and RHPM(
1
) methods, parameters λ
=
0.35 and ϵ
=
1
0
-
4 are used. Note that the parameters are not optimal and are chosen as the best one by tries.Tables1–4 provide numerical results for TM, RHAM2, RHPM(
2
), and RHPM(
1
) methods corresponding to Δ
t
=
0.15, 0.10, 0.05, and 0.01, respectively. In Tables 1–4, column PSNR contains the PSNR values, column Iter contains the number of iterations, and column CPU contains the elapsed CPU time in seconds. Figures 1 to 7 show the images denoised by TM, RHAM2, RHPM(
2
), and RHPM(
1
) methods with Δ
t
=
0.05 for Noisy images with 20
% of Gaussian white noise. Note that the denoised images for other values of Δ
t and 15
% of Gaussian white noise are almost the same. The first row images are Original image and Noisy image with 20
% of Gaussian white noise and the image denoised by TM, and the second row images are the images denoised by RHAM2, RHPM(
2
), and RHPM(
1
) methods. In Figures 1 to 4, Cameraman, Jet-plane, Lake, and Mandrill images with 512
×
512 size are used. In Figure 5, Nodules image with 389
×
366 size is used. In Figures 6 and 7, Rice and Caribou images with 256
×
256 size are used.Table 1
Numerical results for image denoising problems using TM, RHAM2, RHPM(2), and RHPM(1) methods withΔ
t
=
0.15.
Image
Noise
TM
RHAM2
RHPM(2)
RHPM(1)
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
Cameraman
15%
30.24
134
6.58
30.27
88
17.03
30.14
27
3.58
30.33
22
1.36
20%
28.60
188
8.80
28.64
185
35.42
28.69
27
3.54
28.78
22
1.33
Jet-plane
15%
26.79
203
9.85
26.80
219
42.65
27.08
27
3.62
27.16
22
1.36
20%
25.22
287
13.42
25.31
652
123.81
25.54
27
3.55
25.52
22
1.32
Lake
15%
27.35
134
6.50
27.35
88
16.47
27.63
26
3.47
27.65
21
1.31
20%
25.89
191
8.95
25.94
192
35.71
26.33
27
3.52
26.43
22
1.31
Mandrill
15%
26.18
104
5.08
26.17
56
10.74
26.70
25
3.34
26.70
21
1.31
20%
24.62
153
7.31
24.61
113
20.97
25.41
26
3.47
25.38
21
1.26
Nodules
15%
28.20
195
5.16
28.22
199
19.81
28.69
27
1.97
28.67
22
0.76
20%
26.60
271
7.03
26.66
549
54.14
26.79
28
2.01
26.73
22
0.74
Rice
15%
30.62
104
1.34
30.57
57
2.86
31.64
27
0.95
31.85
22
0.41
20%
29.01
147
1.86
28.96
105
5.06
30.16
27
0.94
30.14
22
0.39
Caribou
15%
27.48
125
1.61
27.48
78
3.83
28.42
26
0.94
28.35
22
0.41
20%
26.22
183
2.27
26.23
172
8.02
27.06
27
0.94
27.22
22
0.40Table 2
Numerical results for image denoising problems using TM, RHAM2, RHPM(2), and RHPM(1) methods withΔ
t
=
0.10.
Image
Noise
TM
RHAM2
RHPM(2)
RHPM(1)
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
Cameraman
15%
30.27
200
9.54
30.27
67
13.14
30.19
39
5.14
30.25
31
1.82
20%
28.65
282
13.24
28.56
117
22.44
28.76
40
5.18
28.89
32
1.86
Jet-plane
15%
26.85
304
14.50
26.85
135
26.34
27.10
39
5.14
27.14
32
1.84
20%
25.32
431
20.09
25.31
286
54.55
25.63
40
5.17
25.81
32
1.85
Lake
15%
27.35
201
9.70
27.31
67
12.64
27.60
38
5.02
27.61
30
1.75
20%
25.93
287
13.48
25.90
121
22.69
26.36
39
5.07
26.39
31
1.80
Mandrill
15%
26.15
156
7.58
26.16
48
9.40
26.64
36
4.73
26.56
29
1.70
20%
24.61
228
10.75
24.60
82
15.36
25.38
38
4.90
25.35
31
1.78
Nodules
15%
28.10
289
7.65
28.20
124
13.59
28.80
40
2.89
28.96
32
1.03
20%
26.60
407
10.41
26.62
251
25.04
27.07
40
2.84
27.22
32
1.02
Rice
15%
30.44
154
1.99
30.44
48
2.44
31.72
39
1.55
31.67
31
0.52
20%
28.87
218
2.67
28.83
76
3.67
30.20
40
1.35
30.28
32
0.52
Caribou
15%
27.53
189
2.45
27.56
62
3.08
28.47
38
1.33
28.43
31
0.52
20%
26.17
277
3.33
26.23
113
5.39
27.15
39
1.32
27.17
32
0.52Table 3
Numerical results for image denoising problems using TM, RHAM2, RHPM(2), and RHPM(1) methods withΔ
t
=
0.05.
Image
Noise
TM
RHAM2
RHPM(2)
RHPM(1)
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
Cameraman
15%
30.29
402
19.04
30.16
58
11.73
30.15
76
9.82
30.10
59
3.26
20%
28.57
566
26.26
28.63
85
16.35
28.84
77
9.84
28.88
60
3.32
Jet-plane
15%
26.78
609
28.64
26.82
93
18.16
27.12
77
9.77
27.18
60
3.34
20%
25.38
861
40.20
25.29
150
28.78
25.72
78
9.85
25.87
61
3.35
Lake
15%
27.30
398
18.68
27.35
58
11.06
27.59
73
9.30
27.51
57
3.17
20%
25.95
576
26.73
25.93
86
17.44
26.38
76
9.64
26.45
59
3.24
Mandrill
15%
26.16
311
14.54
26.17
46
8.92
26.61
71
9.07
26.46
56
3.15
20%
24.65
457
21.19
24.59
66
12.49
25.33
74
9.39
25.24
58
3.22
Nodules
15%
28.21
581
14.90
28.19
88
9.04
28.92
78
5.41
29.15
61
1.84
20%
26.68
815
20.58
26.70
139
14.05
27.26
79
5.55
27.67
61
1.82
Rice
15%
30.61
312
3.76
30.53
46
2.32
31.59
76
2.47
31.55
59
0.89
20%
28.92
438
5.25
28.91
63
3.20
30.24
77
2.50
30.27
60
0.88
Caribou
15%
27.49
374
4.53
27.52
55
2.76
28.44
75
2.54
28.42
58
0.87
20%
26.21
552
6.56
26.21
83
4.04
27.20
77
2.53
27.25
60
0.89Table 4
Numerical results for image denoising problems using TM, RHAM2, RHPM(2), and RHPM(1) methods withΔ
t
=
0.01.
Image
Noise
TM
RHAM2
RHPM(2)
RHPM(1)
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
PSNR
Iter
CPU
Cameraman
15%
30.18
1994
92.34
30.22
88
17.01
30.00
367
46.31
29.69
275
14.49
20%
28.68
2823
129.78
28.61
107
20.70
28.73
376
47.41
28.47
281
14.83
Jet-plane
15%
26.85
3033
139.50
26.83
113
23.16
27.08
373
46.50
26.83
280
14.76
20%
25.35
4292
197.00
25.29
142
27.31
25.77
380
47.58
25.75
284
14.88
Lake
15%
27.39
2003
92.20
27.33
88
17.78
27.53
358
44.91
27.30
270
14.25
20%
25.97
2860
131.34
25.92
108
20.30
26.40
369
46.39
26.21
277
14.47
Mandrill
15%
26.19
1554
72.43
26.18
77
14.52
26.53
346
43.47
26.21
261
13.94
20%
24.62
2270
104.40
24.63
95
17.82
25.32
360
45.15
24.93
271
14.23
Nodules
15%
28.18
2900
72.85
28.10
109
11.05
28.94
379
25.60
28.85
284
8.03
20%
26.56
4067
102.17
26.63
136
13.71
27.44
384
25.99
27.77
288
8.14
Rice
15%
30.62
1550
17.99
30.55
77
3.77
31.46
368
11.53
31.06
275
3.74
20%
28.88
2181
25.48
28.98
92
4.44
30.07
376
11.69
29.76
281
3.74
Caribou
15%
27.52
1865
21.67
27.50
85
4.11
28.45
363
11.42
28.10
273
3.65
20%
26.12
2725
31.68
26.22
106
5.10
27.24
373
11.56
27.01
280
3.74Figure 1
Example of the Cameraman image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 2
Example of the Jet-plane image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Example of the Lake image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
Example of the Mandrill image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Example of the Nodules image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Example of the Rice image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Example of the Caribou image: Original image (a), Noisy image (b), TM denoising (c) RHAM2 denoising (d), RHPM(
2
) denoising (e), and RHPM(
1
) denoising (f).
(a)
(b)
(c)
(d)
(e)
(f)As can be seen in Tables1–4, Δ
t
=
0.15 yields the best performance of 4 different cases of Δ
t for TM, RHPM(
2
), and RHPM(
1
) methods, but Δ
t
=
0.05 yields the best performance for RHAM2 method. It can be also seen that RHPM methods perform much faster than TM and RHAM2 methods. On the other hand, PSNR values are close to one another. Overall, RHPM methods give a little bit higher PSNR values than TM and RHAM2 methods. In particular, RHPM gives about 1 higher PSNR value than TM and RHAM2 methods for Rice and Caribou images which have small sizes compared to the other images. According to numerical experiments, this phenomenon occurs because of the shape of image, not the size of image. That is, PSNR values are image-dependent, not size-dependent.Observe that RHPM(
1
) computing 2 terms v
0, v
1 every iteration performs more efficiently than RHPM(
2
) computing 3 terms of v
0, v
1, v
2 every iteration. The reason why RHPM(
2
) performs worse than RHPM(
1
) is as follows: the computational cost for computing v
2 is expensive as compared with v
1, and its computation involves several approximation formulas, which use only a few terms of the infinite series, due to the difficulty in handling highly nonlinear terms. To this end, convergence rate of RHPM(
2
) is slower than RHPM(
1
).


## 6. Conclusion

In this paper, we proposed a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem, and we presented an application of RHPM to TV-based image denoising problem. We also provided numerical methods RHPM(
2
) and RHPM(
1
) which use binomial series techniques to settle the main difficulty in handling nonlinear terms. Numerical experiments show that RHPM(
2
) and RHPM(
1
) are effective, feasible, and reliable for TV-based image denoising problem as compared with TM and RHAM2 methods. Notice that RHPM(
1
) is very efficient compared to other methods since it has smaller iteration number and less CPU time than others. Also RHPM(
1
) is much simpler to implement than RHPM(
2
). So, we recommend the use of RHPM(
1
) with λ
=
0.35 and Δ
t
=
0.15 over RHPM(
2
) for TV-based image denoising problem. This work only studied the effectiveness, feasibility, and reliability of RHPM. So, future work will include comparison results between the proposed RHPM methods and other good methods in the literature.



---

*Source: 207541-2015-09-13.xml*