A NALYSIS AND S IMULATION OF A C OUPLED D IFFUSION BASED
I MAGE D ENOISING M ODEL
arXiv:1907.04526v2 [math.NA] 7 Aug 2019
A P REPRINT
Subit K. Jain
School of Basic Sciences
Indian Institute of Technology Mandi
PIN 175005, INDIA
jain.subit@gmail.com
Sudeb Majee
School of Basic Sciences
Indian Institute of Technology Mandi
PIN 175005, INDIA
sudebmajee@gmail.com
Rajendra K. Ray
School of Basic Sciences
Indian Institute of Technology Mandi
PIN 175005, INDIA
rajendra@iitmandi.ac.in
Ananta K. Majee
Department of Mathematics
Indian Institute of Technology Delhi
PIN 110016, INDIA
majee@maths.iitd.ac.in
August 8, 2019
A BSTRACT
In this study, a new coupled Partial Differential Equation (CPDE) based image denoising model incorporating space-time regularization into non-linear diffusion is proposed. This proposed model is
fitted with additive Gaussian noise which performs efficient image smoothing along with the preservation of edges and fine structures. For this purpose, we propose a new functional minimization
framework to remove the image noise, which results in solving a system of three partial differential
equations (PDEs). Our proposed model is dissimilar from the existing CPDE models as it includes
two additional evolution equations to handle edge strength function and data fidelity term. These
two evolution equations control the smoothing process and force the resultant denoised solution
to be close to the initial solution. To the best of our knowledge, the proposed model is the only
work, which deciphers the combined effect of both the terms using separate PDEs. Furthermore,
we establish the existence and uniqueness of a weak solution of the proposed system using the time
discretization method with H 1 initial data. Finally, we used a generalized weighted average finite
difference scheme to efficiently solve the coupled system and experiment results show the effectiveness of the proposed CPDE model.
Keywords Image denoising · Space-Time regularization · Non-linear diffusion · Semi discretization · Weak solution
1 Introduction
Reduction of generated noise in the restored images without loss of spatial resolution is highly desirable in many
image processing applications. The behavior of noise reduction methods is entirely different for large areas with welldefined borders as well as for poorly defined borders. It has been observed that the goal of noise reduction with feature
preservation and applicability to the multiple acquisition systems is difficult to achieve through a single method [1, 2].
Therefore, the accurate reconstruction of key features in an image from the noisy measurement is the main objective
of our study.
The noise caused by image acquisition devices is often modeled by Gaussian random distribution. Generally, the
degraded image can be modeled as
J = I + η,
(1)
A
PREPRINT
- AUGUST 8, 2019
where J is the noisy image, degraded by Gaussian noise with zero mean, defined on an image space Ω, and I depicts
the clean image, and η is the additive noise function. Our main task is to obtain the reconstructed image from the
observed noisy image.
In recent years, PDE based methods have become widely applicable for noise removal and signal reconstruction, due
to their well-established mathematical properties and their general wellposedness [1, 3, 4, 5, 6, 7, 8]. Obtaining the
steady-state solution is the basic objective of PDE based image restoration problem. This goal can be achieved through
two approaches. In the first approach, we directly find a steady state solution of the diffusion model. Whereas in the
second approach, the steady-state solution is obtained from the evolution equation of the Euler-Lagrange equation
associated with the energy functional. Several techniques are employed to handle the challenges posed by restoration
problems. The availability of vast literature does not allow us to review all of them. However, we briefly report some
of the robust techniques such that anisotropic diffusion [3, 6, 7, 9, 10, 11], fourth order diffusion [12, 13, 14], complex
diffusion [15, 16] and variational based methods [17, 18, 19, 20, 21, 22]. In all of the above-mentioned methods, spatial
regularization is used for diffusion function, which is not able to inject the past information into the iterative process
of diffusion. Therefore, the scheme may not be effective for solving image denoising problems which contain huge
textures in their image domain. To handle this, applications of time-delay regularization with spatial regularization
are also available in the literature [23, 24, 25, 26, 27, 28]. The space-time regularization is more suitable approach
than Gaussian filtering as it enables to incorporate the information obtained along with the scales into the diffusion
process. This type of regularization is used to recover fine structures, especially for images degraded with higher noise
levels. Hence apart from a single PDE based approaches, several studies are utilizing coupled PDE to perform image
restoration [29]. In this regard, Nitzberg and Shiota [27] introduce a delay in time to calculate diffusivity term as,
1
It = ∇
∇I ,
in Ω × (0, T ),
(2)
1 + u2
where u in diffusion function is calculated using the relation: ut = ϕ(Gξ ∗ |∇I|2 − u) and T is a fixed time. Here,
ξ is a positive smoothing parameter required for the analytical well-posedness of the model. So, the term u updates
the edge information from the gradient of the updated image and past information of itself. Hence, instead of spatial
regularization, u represents a temporal regularization in the formation of the diffusion term. This coupled PDE system
is quite similar to the Perona-Malik model [10] as there is no spatial smoothing term. In [27], the authors show that
the coupled PDE system admits a unique classical solution (I, u) in any dimension and satisfy the maximum principle.
They also explain that the system does not produce spurious features. A similar idea of coupled PDE system using
time-delay regularization is proposed by Luo et al. [26], for image denoising. They propose the estimation of a better
edge map by substituting the isotropic diffusion by a nonlinear diffusion equation, to produce a regularized version of
the diffusion function for image smoothing. Subsequently, Guo et al. [25] propose and analyze the well-posedness of
a reaction-diffusion system for image denoising. They deploy the H −1 norm instead of L1/L2 norm with the total
variation model [30], to preserve the oscillatory and texture patterns, as follows,
Z
(3)
|∇I| dΩ + λ||J − I||2H −1 (Ω) .
EH −1 =
Ω
R
Here,
= Ω |∇∆ (.)| dΩ. This model is efficient for denoising and decomposing the image into cartoon
plus texture. Now evolution equation of the Euler-Lagrange equation of the above functional, in which two PDEs are
interacting with each other to calculate noise fidelity term, can be given as,
∇I
It = ∇
− 2λv,
in Ω × (0, T ),
|∇I|
(4)
vt = ∆v − (J − I),
in Ω × (0, T ).
||.||2H −1 (Ω)
−1
2
In the above model, another PDE is used to calculate the fidelity (v) term efficiently. The main difficulty with this
1
model is that the term |∇I|
is singular when |∇I| = 0. However, the numerical simulations involve the regularization
of |∇I| with |∇I| + ε. A good survey on the existing CPDE approaches can be found in [29].
The objective of the present work is to systematically develop a new framework of coupled nonlinear diffusion model.
The present approach differs from previous existing CPDE based works as it utilizes separate evolution equations to
calculate edge strength function as well as the data fidelity term. The originality of the proposed coupled non-linear
PDEs algorithm lies in the combined effect of:
i. diffusion coefficient with space-time regularization instead of space regularization only;
ii. H −1 norm based time-varying fidelity term.
2
A
PREPRINT
- AUGUST 8, 2019
Moreover, due to the interest in theoretical research, we establish the well-posedness of the proposed CPDE model
(10)-(13). To prove the well-posedness, we use a time discretization method along with essential a-priori estimates
and classical results of compact inclusion in Sobolev spaces [31]; by utilizing some important inequalities, lemmas
and theorems [32] and then tends the time discretization parameter to zero. This imparts a major role of this study as
it is very important for the numerical computation. Further, to simulate the image denoising, we apply a higher order
accurate generalized weighted average finite difference scheme with an advanced iterative solver (Hybrid Bi-Conjugate
Gradient Stabilized method) to solve the algebraic system of equations generated from the discretization process [33].
The image quality of denoised images using the proposed CPDE model is compared against several existing nonlinear diffusion models, variational based models and CPDE based models. In practice, our calculations indicate that
the proposed approach allows effective image smoothing on fine scales even when the images are degraded with a
higher level of noise.
The rest of the paper is arranged in the following sequence. Section 2 describes the mathematical formulation of the
proposed coupled nonlinear diffusion model. In section 3 we establish the existence and uniqueness of a weak solution
to the proposed PDE model. Section 4 shows an appropriate numerical realization in terms of higher order accurate
implicit finite difference scheme. Numerical validation of the proposed approach using experimental study is carried
out in section 5. At the end, we summarize our observations in section 6.
2 Proposed Coupled Diffusion Model
The main aim of this work is to present a novel approach for image smoothing by the evolution of coupled non-linear
PDEs. Among the existing CPDE based models mentioned in the previous section, the choices of diffusion function
and fidelity term include the space and time regularization. Numerical simulations and experimental results depict that
the space-time regularization plays a important role in the quality of recovered images degraded by noise. We note
that all of these existing coupled PDEs described in the last section utilizes the time-delay regularization to calculate
either diffusion coefficient or fidelity term, and sometimes fail to preserve significant contents of the image. Thus, one
issue of these smoothing algorithms is how to design a non-linear diffusion algorithm, which can remove the noise and
preserve the small features, simultaneously. For this purpose, inspired by the impressive performance of time-delay
regularization in image smoothing, we proposed the following energy functional:
Z
1
minEu (I) =
g(u)|∇I|2 dΩ + λ||I0 − I||2H −1 (Ω) ,
(5)
2
Ω
where I0 represents the high pass filtered version of image J and λ is weight parameter. Generally, λ is chosen to be
inversely proportional to the variance of the noise in the given image. The first term represents a regularizing term
producing a smooth varying variable function. The term g(u) serves the purpose of detecting the edges in the image,
and u is the edge variable. The second term forces the denoised solution I to be a close approximation to I0 and
referred to as data fidelity term. In the data fidelity term, we have replaced L2 norm of (I0 − I) by H −1 norm to
preserve the oscillatory and texture information more appropriately in a denoised image. The solution of the above
minimization problem (5) can be given by the steady-state of the Euler-Lagrange equation for Eu (I):
1
∇(g(u)∇I) = ∆−1 (I0 − I).
2λ
(6)
This is equivalent to following equations,
∇(g(u)∇I) − 2λv = 0,
∇(∇v) − (I0 − I) = 0.
(7)
To this end, we consider the function u as the following,
ut = α(|∇Iξ |2 − u +
β2
∆u),
2
(8)
where α > 0 and β > 0 are parameters to be specified; average time delay (α−1 ) and amount of the spatial smoothing
(β) control the deformation and smoothness of the edge information and chosen as described in [27]. And, Iξ =
Gξ (X) ∗ I where two dimensional Gaussian function Gξ (.) at each pixel, X ∈ Ω is adopted as,
√ −1
−|X|2
Gξ (X) = ξ 2π
,
(9)
exp
2ξ 2
3
A
PREPRINT
- AUGUST 8, 2019
where “ ∗ ” denote the convolution operator and ξ > 0 is the standard deviation of the function. Finally, our proposed
model in which a set of PDEs interacting with each other, can be expressed as follows,
∂I
= ∇(g(u)∇I) − 2λv,
in ΩT := (0, T ) × Ω,
∂t
β2
∂u
= α(h(|∇Iξ |2 ) − u +
△ u),
in ΩT ,
∂t
2
∂v
= ∇(∇v) − (I0 − I),
in ΩT ,
∂t
I(0, x) = I0 (x),
v(0, x) = 0,
u(0, x) = Gξ ∗ |∇I0 |2 ,
in Ω.
(10)
(11)
(12)
(13)
In the above model, diffusion coefficient g(u) is chosen as
g(u) =
1
1+
|Gξ ∗u|
k2
,
(14)
where k > 0 is a threshold parameter, and u represents edge strength at each iteration. Note that we replace |∇Iξ |2
in (8) by h(|∇Iξ |2 ), where h being a sort of truncation, see [24]. Our proposed model is different from the models
reported earlier in [29], as it yields two separate evolution equations (PDEs) to handle the diffusion coefficient and data
fidelity term. These two evolution equations are responsible for stopping the image smoothing at the edges as well as
textures and forces the resulting solution to be a close approximation to the given initial image. The proposed model
is in a spirit similar to the reaction-diffusion model equation in [25]; however, the diffusion function is very different
in our case. In (10)-(13), the smoothing equation (10) has different edge variable u and fidelity term v, obtained from
two different PDEs. The data fidelity term between I and I0 can be handled by function v, which can be obtained
from equation (12). Whereas the function u in the diffusion coefficient can be calculated from equation (11). In
summary, the proposed model achieves a suitable edge map and fidelity between noisy image and restored image at
each iteration, which ultimately leads to quality denoising results. Hence, the proposed model provides a potential
approach for image noise removal and enhancement.
3 Existence and uniqueness of weak solution
In this section, we establish the well-posedness of (10)-(13) with Dirichlet’s boundary condition for the fidelity variable
and Neumann’s boundary condition for other two variables. Furthermore, for simplicity we choose all the constants
involved in the equations (10)-(13) equals to 1. Note that, h : R+ → R+ is a bounded, Lipschitz continuous function
with Lipschitz constant ch such that
δ ≤ h(ũ) ≤ 1
∀ ũ ∈ R+ ,
(15)
for some δ > 0. Also from (14) we observe that, g : R → R+ is a bounded, decreasing and Lipschitz continuous
C
function with Lipschitz constant k2ξ . Moreover, g(0) = 1 and lim g(u) = 0.
u→+∞
3.1 Technical framework
We denote by H k (Ω), k is a positive integer, the set of all functions I : Ω → R such that I and its distributional
P2
∂mI
2
k
derivatives ∂x
m of order |m| =
j=1 mj ≤ k all belongs to L (Ω). H (Ω) is a Hilbert space endowed with the
norm
X Z ∂ m I 2 1/2
kIkH k := ||I||H k (Ω) =
dx
.
m
Ω ∂x
|m|≤k
We denote by Lp (0, T ; H k (Ω)), p > 1 and k is a positive integer, the set of all measurable functions I : [0, T ] →
H k (Ω) such that
Z T
1/p
kIkLp(0,T ;H k (Ω)) :=
||I(t)||pH k (Ω) dt
< ∞.
0
We denote by H (Ω) the dual of H (Ω), and H (Ω) the dual of H01 (Ω). For any f ∈ H 1 (Ω)′ , we define a norm as
n
o
||f ||H 1 (Ω)′ = sup hf, ui : u ∈ H 1 (Ω) , ||u||H 1 (Ω) ≤ 1 .
1
′
1
−1
4
A
PREPRINT
- AUGUST 8, 2019
We introduce the solution space W (0, T ) := W1 (0, T ) × W1 (0, T ) × W2 (0, T ) for the problem (10)-(13), where
n
o
∂w
W1 (0, T ) = w : w ∈ L2 (0, T ; H 1(Ω)) ,
∈ L2 (0, T ; (H 1 (Ω))′ ) ,
∂t
o
n
∂w
2
1
∈ L2 (0, T ; H −1 (Ω)) .
W2 (0, T ) = w : w ∈ L (0, T ; H0 (Ω)) ,
∂t
Note that, the space Wi (0, T ) (i = 1, 2) is a Hilbert space for the graph norm, see [34].
Definition 3.1 (Weak solution) Let I0 ∈ H 1 (Ω). We say that a triple (I, u, v) is a weak solution of (10)-(13) if
(i) I, u ∈ W1 (0, T ), and v ∈ W2 (0, T ).
(ii) For all Ψ ∈ H 1 (Ω) and Φ ∈ H01 (Ω), there holds
Z T
Z
Z
∂I
vΨ dx dt = 0 ,
g(u) ∇I · ∇Ψ dx dt + 2
, Ψ dt +
∂t
0
ΩT
ΩT
Z T
Z
Z
∂u
uΨ dx dt
h(|∇Iξ |2 )Ψ dx dt +
, Ψ dt −
∂t
0
ΩT
ΩT
Z
∇u · ∇Ψ dx dt = 0 ,
+
Z
0
T
∂v
, Φ dt +
∂t
Z
ΩT
∇v · ∇Ψ dx dt +
Z
(16)
(17)
ΩT
ΩT
(I0 − I)Ψ dx dt = 0 .
(18)
(iii) (13) holds.
3.2 Discretized system and existence of its weak solution
To prove existence of a weak solution of (10)-(13), we use semi discretization in time with a semi-implicit Euler
T
method. Let 0 = t0 < t1 < t2 · · · < tN = T be a uniform partition of [0, T ] with time-step size τ = N
for some
+
2
1
1
N ∈ N , i.e, tn = τ n for 0 ≤ n ≤ N . For simplicity, in the following we write L , H , H0 instead of L2 (Ω),
H 1 (Ω), H01 (Ω) respectively. We may then consider time discretization of (10)-(13) as follows: for 0 ≤ n ≤ N − 1 ,
∀ Ψ ∈ H 1 , and ∀ Φ ∈ H01 iterate:
(i) compute un+1 ∈ H 1 such that
1
(un+1 − un ), Ψ L2 + ∇un+1 , ∇Ψ L2 + un+1 , Ψ L2
τ
= h(|∇Gξ ∗ In |2 ), Ψ L2 ,
(ii) compute vn+1 ∈
H01
(iii) compute In+1 ∈ H 1
such that
1
1
vn+1 , Φ L2 + ∇vn+1 , ∇Φ L2 =
vn − (I0 − In ), Φ L2
τ
τ
1
1
In+1 , Ψ L2 + g(un )∇In+1 , ∇Ψ L2 =
In − 2vn , Ψ L2 .
τ
τ
(19)
(20)
(21)
Thanks to (15), the existence of un+1 resp. vn+1 in step (i) resp. in step (ii) easily follows from Lax-Milgram lemma.
Moreover, these solutions are unique. To show the well-posedness of (21), we use standard Lax-Milgram lemma.
Define a continuous bilinear form C : H 1 × H 1 → R, and a linear bounded functional ℓ : H 1 → R as
1
1
C(Θ, Ψ) =
Θ, Ψ L2 + g(un )∇Θ, ∇Ψ L2 ; ℓ(Ψ) =
In − 2vn , Ψ L2 .
(22)
τ
τ
Then, by Lax-Milgram lemma, there exists a unique In+1 ∈ H 1 satisfying (21) provided C is H 1 -coercive. Indeed,
thanks to Gagliardo-Nirenberg inequality,there exists a constant C > 0, depending only on Ω and Gξ , such that
1
≤ g(un ) ,
(23)
C||un ||H 1
1+
k2
o
n
1
1
2
2
e
e := min 1 ,
yielding
||∇Ψ||
≥
C||Ψ||
where
C
and hence C(Ψ, Ψ) ≥ τ1 ||Ψ||2L2 +
2
1
L
H
τ
C||un ||H 1
C||un ||H 1
1+
1+
k2
k2
the H 1 -coercivity of C.
5
A
PREPRINT
- AUGUST 8, 2019
3.3 A-priori estimates
Let 0 ≤ n ≤ N − 1. Choose Ψ = un+1 as test function in (19) and use the algebraic identity
ha, a − bi =
1
|a|2 − |b|2 + |a − b|2 ,
2
∀ a, b ∈ Rp (p ≥ 1) ,
(24)
Cauchy-Schwarz and Young’s inequalities, (15), and then sum over n = 0, 1, 2, . . . ,
j − 1 with 1 ≤ j ≤ N . The result is
||uj ||2L2
+
j−1
X
n=0
||un+1 −
un ||2L2
j−1
X
≤ ||u0 ||2L2 + |Ω|
n=0
+ 2τ
j−1
X
n=0
||∇un+1 ||2L2
τ ≤ ||u0 ||2L2 + T |Ω| .
(25)
Hence, there exists a constant C1 > 0, independent of τ > 0, such that
max ||uj ||2L2 +
1≤j≤N
N
−1
X
n=0
||un+1 − un ||2L2 + 2τ
N
−1
X
n=0
||∇un+1 ||2L2 ≤ C1 .
(26)
Similarly, by choosing (formally) Ψ = −∆un+1 in (19) and using the integration by parts formula, (24), CauchySchwarz and Young’s inequalities, (15), and then summing over n = 0, 1, 2, . . . , j − 1 with 1 ≤ j ≤ N , we obtain
max
1≤j≤N
||∇uj ||2L2
+
N
−1
X
n=0
||∇(un+1 −
un )||2L2
+τ
N
−1
X
n=0
||∆un+1 ||2L2 ≤ C2 ,
(27)
where C2 > 0 is a constant, independent of τ > 0. In view of (26) and (27), there exists a constant C3 > 0,
independent of τ , such that
(28)
max kuj kH 1 (Ω) ≤ C3 .
0≤j≤N
Therefore, from (23), the positive lower bound for g(un ) is given by
ν :=
1
1+
CC3
k2
(29)
≤ g(un ) (0 ≤ n ≤ N ) .
Again, one can use the test function Φ = vn+1 and Φ = −∆vn+1 (formally) in (20), and proceed as above (under the
cosmetic changes) to obtain (1 ≤ j ≤ N ):
||vj ||2L2 +
j−1
X
n=0
||vn+1 − vn ||2L2 + τ
≤
||∇vj ||2L2 +
j−1
X
n=0
j−1
X
n=0
kv0 k2L2
||∇vn+1 ||2L2
+
||∇(vn+1 − vn )||2L2 + τ
CP2 T kI0 k2L2
j−1
X
n=0
+
τ CP2
j−1
X
n=0
kIn k2L2 ,
(30)
||∆vn+1 ||2L2
≤ ||∇v0 ||2L2 + T kI0 k2L2 + τ
j−1
X
n=0
||In ||2L2 ,
(31)
where CP is the Poincaré constant. Note that, because of coupled system, we are not able to find the bound of vn ’s
from (30) and (31). May be we need to combine with the estimates coming from (21) by choosing appropriate test
6
A
PREPRINT
- AUGUST 8, 2019
functions. Taking Ψ = In+1 in (21) and using (24), (29), and (30), we have
||Ij ||2L2 +
j−1
X
n=0
||In+1 − In ||2L2 + 2τ ν
≤ ||I0 ||2L2 + 2τ
j−1
X
n=0
j−1
X
n=0
≤ C4 + C5 τ
k=0
n=0
||∇In+1 ||2L2
||In+1 ||2L2
+ 2τ
j−1
X
j−1
X
kv0 k2L2 + CP2 T kI0 k2L2 + τ CP2
n
X
k=0
kIk+1 k2L2
kIk+1 k2L2 ,
(32)
for some constants C4 , C5 > 0, only depend on T, CP , v0 , I0 . We now apply discrete Gronwall’s lemma (implicit
form) and obtain the following: there exists τ1 > 0 such that for all time step sizes 0 < τ < τ1 ,
max ||Ij ||2L2 +
0≤j≤N
N
−1
X
n=0
||In+1 − In ||2L2 + 2τ ν
N
−1
X
n=0
||∇In+1 ||2L2 ≤ C6 ,
(33)
where C6 > 0 is independent of τ > 0. We now use (33) in (30) and (31), and get
max ||vj ||2H 1 +
0≤j≤N
N
−1
X
n=0
||vn+1 − vn ||2H 1 + τ
||∇vn+1 ||2L2 + k∆vn+1 k2L2 ≤ C7 ,
N
−1
X
n=0
(34)
for some constant C7 > 0, independent of τ with 0 < τ < τ1 . Again one may choose Ψ = −∆In+1 (formally) in (21)
and use the integration by partsformula, (24), Cauchy-Schwarz and Young’s inequalities, (29), (28), (34) and the fact
that k∇g(un )kL∞ ≤ C ξ, Ω, k kun kH 1 to get, after the application of discrete Gronwall’s lemma (implicit form),
max ||∇Ij ||2L2 +
1≤j≤N
N
−1
X
n=0
||∇(In+1 − In )||2L2 + τ
j−1
X
n=0
||∆In+1 ||2L2 ≤ C8 ,
(35)
with C8 > 0, independent of τ with 0 < τ < τ2 for some τ2 > 0. Note that, for τ0 = min{τ1 , τ2 }, (33), (34) and (35)
all hold true. Putting things together, we arrive at the following theorem.
Theorem 3.1 Let I0 ∈ H 1 . Then there exist constants τ0 > 0 and C > 0 such that for all time step sizes τ > 0 with
0 < τ < τ0 , there holds
max kuj k2H 1 + kvj k2H 1 + kIj k2H 1
1≤j≤N
+
N
−1
X
n=0
+τ
||un+1 − un ||2H 1 + ||vn+1 − vn ||2H 1 + ||In+1 − In ||2H 1
k∆un+1 k2L2 + k∆vn+1 k2L2 + k∆In+1 k2L2 ≤ C ,
N
−1
X
n=0
(36)
where uj , vj and Ij solves (19), (20) and (21) respectively.
3.4 Continufication and existence of weak solution
Let {tn }N
n=0 be a uniform partition of [0, T ] with time-step size τ as described in Subsection 3.2. For any sequence
xn+1 − xn
{xn } ⊂ X, where X is a Banach space, we define the difference quotient dt xn+1 =
for 0 ≤ n ≤ N − 1.
τ
The globally time interpolant Xτ ∈ C(X) of {xn } is defined via
Xτ (t) :=
t − tn
tn+1 − t
xn+1 +
xn ,
τ
τ
7
∀ t ∈ (tn , tn+1 ] .
A
PREPRINT
- AUGUST 8, 2019
−
Moreover, we define the piecewise constant in time interpolants X+
τ (t) and Xτ (t) as follows:
X+
τ (t) := xn+1 ,
X−
τ (t) := xn ,
∀ t ∈ (tn , tn+1 ] ,
with x−1 = x0 and xN +1 = 0. With the above notation, we now show the boundedness of the sequence {∂t Uτ }τ >0
in L2 (H 1 (Ω)′ ). To do so, let Ψ ∈ H 1 (Ω)′ \ {0}, and t ∈ (tn , tn+1 ] where 0 ≤ n ≤ N − 1. Then from (19), we have,
by using Cauchy-Schwarz inequality and (15),
1
∂t Uτ , Ψ ≤ |Ω| 2 + max kun+1 kH 1 kΨkH 1 ,
(37)
0≤n≤N −1
and hence, thanks to Theorem 3.1,
−1
NX
|Ω| + T
||∂t Uτ ||2L2 (H 1 (Ω)′ ) ≤C τ
n=0
max
0≤n≤N −1
kun+1 k2H 1 ≤ C .
(38)
Similarly, the sequences {∂t Vτ }τ >0 and {∂t Iτ }τ >0 are bounded in L2 (H −1 (Ω)) and L2 (H 1 (Ω)′ ) respectively. Moreover, there exists a constant C > 0 such that for all 0 < τ < τ0 ,
||∂t Vτ ||2L2 (H −1 (Ω)) + ||∂t Iτ ||2L2 (H 1 (Ω)′ ) ≤ C .
(39)
Thanks to (38), (39) and Theorem 3.1, by using classical results of compact inclusion in Sobolev spaces [31], there
exists (I, u, v) ∈ W1 (0, T ) × W1 (0, T ) × W2 (0, T ) such that along a subsequence (still we denote it by the same
index) the following hold:
+ − ∗
Iτ , Iτ ⇀ I in L∞ (H 1 ) ; ∂t Iτ ⇀ ∂t I in L2 (H 1 (Ω)′ ) ,
Iτ , Iτ± ⇀ I in L2 (H 1 ) ; Iτ , Iτ± → I in L2 (L2 ) ,
|∇Gξ ∗ Iτ− |2 → |∇Gξ ∗ I|2 in L2 (L2 ) and a.e. (0, T ) × Ω ,
2
2
2
2
h(|∇Gξ ∗ Iτ | ) → h(|∇Gξ ∗ I| ) in L (L ) ,
+
−
∗
∞
1
(40)
Uτ , Uτ ⇀ u in L (H ) ; ∂t Uτ ⇀ ∂t u in L2 (H 1 (Ω)′ ) ,
(τ → 0) ,
±
2
1
±
2
2
U
,
U
⇀
u
in
L
(H
)
;
U
,
U
→
u
in
L
(L
)
,
τ
τ
τ
τ
g(Uτ+ ) → g(u) in L2 (L2 ) ,
V + , V − ⇀∗ v in L∞ (H01 ) ; ∂t Vτ ⇀ ∂t v in L2 (H −1 (Ω)) ,
τ ±τ
Vτ , Vτ ⇀ v in L2 (H01 ) ; Vτ , Vτ± → v in L2 (L2 ) .
For any Ψ ∈ H 1 and Φ ∈ H01 , we rewrite (19), (20) and (21) in terms of Uτ , Uτ+ , Uτ− , Vτ , Vτ+ , Vτ− , Iτ , Iτ+ , and Iτ− ,
integrate with respect to the time variable and make use of (40) to pass to the limit as τ → 0 in the resulting variational
formulation to arrive at (16), (17) and (18). In other words, the triple (I, u, v) is a weak solution of (10)-(13).
3.5 Uniqueness of weak solution
We use standard methodology [32] to prove the uniqueness of weak solution of (10)-(13). Let (I1 , u1 , v1 ) and
(I2 , u2 , v2 ) be two sets of solution for the system (10)-(13) with I1 6= I2 , u1 6= u2 and v1 6= v2 . Let
e := u1 − u2 , Ve := v1 − v2 and gi = g(ui ) (i = 1, 2). Then, the following equations hold in
I˜ := I1 − I2 , U
the sense of distribution:
∂ I˜
˜ = ∇((g1 − g2 )∇I2 ) − 2Ve ,
− ∇(g1 ∇I)
in ΩT ,
(41)
∂t
e
∂U
e = h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 ) − U
e,
− ∆U
in ΩT ,
(42)
∂t
∂ Ve
− ∆Ve = I˜ ,
in ΩT ,
(43)
∂t
˜ x) = 0 = U
e (0, x) = Ve (0, x) ,
I(0,
in Ω .
(44)
Note that
1
ν1 :=
1+
Cku1 kL∞ (H 1 )
k2
≤ g(u1 ) ;
||(g1 − g2 )||L∞ ≤ Cξ ||(u1 − u2 )||L2 .
8
A
PREPRINT
- AUGUST 8, 2019
Similar to a-priori estimates, we have for a.e. t ∈ (0, T ) with C̃I := k∇Ik2L∞ (L2 ) ,
1 d ˜ 2
˜ 2 2 + kVe k2 2
kIkL2 ≤ C̃I C(ν1 )||(g1 − g2 )||2L∞ + kIk
L
L
2 dt
˜ 2 2 + ||Ve ||2 2 .
e U
e ||2 2 + ||I||
≤ C||
L
L
L
e and integrating over Ω, we get
Similarly, by multiplying (42) by U
d e 2
e ||2 2 .
e ||2 2 ≤ C||h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 )||2 2 + ||U
||U ||L2 + 2||∇U
L
L
L
dt
By using Lipschitz continuity of h along with Gagliardo-Nirenberg inequality, we have
and hence
˜ 22 ,
||h(|∇Gξ ∗ I1 |2 ) − h(|∇Gξ ∗ I2 |2 )||2L2 ≤ C(ch , ξ, I1 , I2 )kIk
L
In a similar way, one can easily deduce that
d e 2
˜ 22 .
||U ||L2 ≤ C||I||
L
dt
(45)
(46)
(47)
(48)
d e 2
˜ 2 2 + ||Ve ||2 2 .
(49)
||V ||L2 ≤ C ||I||
L
L
dt
e ≡ 0 and Ve ≡ 0. In other words,
We now add (45)-(49) and then apply Gronwall’s inequality to infer that I˜ ≡ 0, U
weak solution of the proposed model (10)-(13) is unique.
4 Numerical Approximation
In this section, we discuss the numerical implementation of the proposed model (10)-(13). To solve the model numerically, we have applied the generalized weighted average finite difference scheme which has a combined nature of the
forward Euler method and the backward Euler method, at nth and (n + 1)th step, respectively. The image variable I,
edge variable u and fidelity term v are calculated at each iteration.
In the derivation of finite difference formulas, h and τ are considered as spatial step size and time step size, respectively.
n
Let I(xi , yj , tn ) express the gray level of the image plane Ii,j
, where xi = ih, yj = jh and tn = nτ . Also, we consider
0
0 2
vi,j
= 0 and u0i,j = Gξ ∗ |∇Ii,j
| . Derivatives are approximated by central difference formula as,
n
n
n
n
Ii+1,j
− Ii−1,j
Ii,j+1
− Ii,j−1
∂I
∂I
≈
,
≈
,
∂x
2h
∂y
2h
n
n
n
n
n
n
Ii+1,j
− 2Ii,j
+ Ii−1,j
Ii,j+1
− 2Ii,j
+ Ii,j−1
∂2I
∂2I
≈
,
≈
.
∂2x
h2
∂2y
h2
The discrete form of (10)-(13) can be written as;
τ
τ
n+1
n+1
n
n
Ii,j
−
,
(∇(g(u)∇I)n+1
= Ii,j
+
(∇(g(u)∇I)ni,j − 2λvi,j
i,j − 2λvi,j
2
2
n
un+1
ψ2
i,j − ui,j
n+1
,
+
= ϕ h(|∇Iξni,j |2 ) − un+1
∆u
i,j
i,j
τ
2
n+1
n
vi,j
− vi,j
n+1
n
0
,
= ∆vi,j
− Ii,j
− Ii,j
τ
where,
n
n
n
n
n
n
wi+1,j
− 2wi,j
+ wi−1,j
wi,j+1
− 2wi,j
+ wi,j−1
+
,
h2
h2
for w = u, v. Moreover the diffusion term in (10)-(13) approximated using the central difference scheme as,
n
∆wi,j
=
∇(g(u)∇I)i,j =
1
(g 1 (Ii+1,j − Ii,j ) − gi− 21 ,j (Ii,j − Ii−1,j )+
h2 i+ 2 ,j
gi,j+ 12 (Ii,j+1 − Ii,j ) − gi,j− 12 (Ii,j − Ii,j−1 )) ,
9
(50)
(51)
(52)
(53)
(54)
A
PREPRINT
- AUGUST 8, 2019
with the Neumann boundary conditions:
Ii,0 − Ii,1 = 0,
Ii,N −1 − Ii,N = 0,
for 0 ≤ i ≤ M,
I0,j − I1,j = 0,
IM−1,j − IM,j = 0,
for 0 ≤ j ≤ N.
Further, to solve the algebraic system of the the form AI n+1 = B, generated from numerical discretization, hybrid
Bi-Conjugate Gradient Stabilized solver has been used [33, 37].
Apart from the numerical discretization of (10)-(13), a stopping criterion is needed to terminate the diffusion process.
To achieve this goal, we start with an initial image I0 and apply the system (10)-(13) repeatedly. This results in a
family of smoother images I(t, x); t > 0, which depicts refined versions of I0 . After sufficient iterations, changes
between two consecutive iterations become redundant. At this point, the convergence of the iterative procedure has
been achieved. Now, we use the following measure as a stopping criterion,
2
||I k+1 − I k ||2
2
||I k ||2
(55)
≤ ε.
th
Here, we use ε = 10−4 as a fixed threshold. I k and I (k+1) depict the image planes at the k th and (k + 1) iteration,
respectively. Parameter values used for the numerical experiments are mentioned in the caption of each figure. For
Algorithm 1 Algorithm of CPDE model for Image Denoising
1: Inputs: I (0) = I0 as initial solution or initial noisy image,
v (0) = 0 as initial edge variable and u(0) = Gξ ∗ |∇I (0) |2 , τ as time step, k = 0.
2: Result: Denoised image (I)
3: Calculate I (k+1) , u(k+1) and v (k+1) from
k+1 k τ
0
τλ
1 − τ2 A(I k+1 )
I + 2 (A(I k ) − 2λv k )
I
2
k+1
= uk + τ ϕh(||∇Iξk ||2 )
0
1 + τ (1 − ψ2 ∆)
0 u
v k − τ (I k − I 0 )
v k+1
0
0
1 − τ∆
.
4: Solve the algebraic system (obtained using step 3), using Hybrid BiCGStab solver .
5: Check if
2
||I k+1 − I k ||2
≤ ε = 10−4
||I k ||22
then stop or go to Step 3.
6: end
the TV and NLM model, we have reproduced the results by utilizing the fact as mention in [19] and [39] respectively.
Apart from the parameters displayed in the captions, we have chosen ϕ = 1 and ξ = 1 for the present model and a
uniform time step size τ = 0.1 for the current as well as the other discussed approaches.
5 Results and Discussion
The proposed CPDE model was employed using the finite difference scheme given in the previous section. Image denoising using (10)-(13) was compared with the results of other state-of-art models available in the literature. Especially,
TV model [19], NS model [27], Luo model [26], RD model [25], PV model [28], SYS model [38] are considered for
the comparison. Apart from the diffusion-based approaches, non-iterative approach (Non-Local Means method) [39]
is also used for comparison with the present method. Since the space-time regularization based proposed model is
claimed to be an improvement over the existing CPDE models, our main aim was to compare the edge detection and
denoising results with these models (NS Model, and Luo model). In this process, the considered non-linear diffusion
models are solved by the existing numerical schemes. And also, to stop the iterative process of each smoothing algorithm discussed stopping criteria equation (55) is employed. The effectiveness of results was evaluated through several
standard gray level test images which are degraded with additive Gaussian noise of zero mean and different levels of
standard deviations. We have artificially added additive Gaussian noise of different standard deviations ranging from
20 to 50 by using our MATLAB program. Also, to evaluate the efficiency of our model, quantitative comparisons
in terms of Peak signal to noise ratio (PSNR)[2], mean structural similarity index measure (MSSIM)[40], Gradient
of peak signal to noise ratio (PSNRGrad )[41], Improvement in signal to noise ratio (ISNR)[2] have been shown with
existing models. A higher value of quantity metrics suggests that the filtered image is closer to the noise-free image.
10
A
PREPRINT
- AUGUST 8, 2019
ˆ
(a). Peak signal to noise ratio (PSNR)[2] can measure the match between the clean (I) and denoised image (I),
ˆ = 10log
PSNR(I, I)
10
MN|max(I) − min(I)|2
.
2
||Iˆ − I|| 2
L
(b). In addition to PSNR, gradient of peak signal to noise ratio (PSNRGrad ) is also used to measure the match between
the derivatives of reconstructed and true image, defined as
PSNRGrad =
1
(PSNR(Ix , Iˆx ) + PSNR(Iy , Iˆy )) ,
2
ˆ
where (Ix , Iy ) and (Iˆx , Iˆy ) are the derivatives of ground truth image (I) and denoised image (I).
(c). Structural similarity index (SSIM) [40], is used to calculate the similarity between structure of clean and reconstructed images, and can be given as,
SSIM(X, Y ) =
(2µx µy + c1 )(2σxy + c2 )
,
(µ2x + µ2y + c1 )(σx2 + σy2 + c2 )
here µx , µy , σx2 , σy2 , σxy are the average, variance and covariance of X and Y , respectively. The variables c1 and c2
are used to stabilize the division with weak denominator. To find overall performance of the image, we used a mean
SSIM index (MSSIM),
N
MSSIM(X, Y ) =
1X
SSIM(Xi , Yi ) ,
N i=1
where X and Y are used to represent the original and the denoised images, respectively. N is the number of local
windows in the image, and Xi and Yi are the image contents at the ith local window.
(d). Improvement in signal to noise ratio (ISNR) [2], is simply the difference between the improved and the original
signal to noise ratio,
M N
PP
2
(I
−
I
)
0
i=1j=1
ISNR = 10log10 M N
.
P P
2
(I − It )
i=1 j=1
To test the effectiveness of the proposed model, figure[1] presents the denoising results for Walkbridge test image
that contain additive Gaussian noise of σ = 40, which make the features hard to visualize. The restored outputs
obtained from the state-of-the-art diffusion models depict that the denoised images are not satisfactory and the texture
information of the image is ruined. From the figure[1e], it is easy to perceive that the proposed model can effectively
eliminate all noise particles and preserve the original structure of the image.
Further, to showcase the efficiency of our model, in figure 2, we compare the visual quality of the mosaic image
computed by the present model along with the other existing models. Here we present the denoised images and 3D
surface plots of the denoised images side by side for better comparison. From the figure, it can be observed that
our model works better in terms of noise removal as compared to other models. From the 3D surface plots of the
denoised images, one can see that our model leave fewer fluctuations with compare to the other models, it indicates
that the proposed model not only removes noise efficiently but also preserves the fine structure as compared to the
other models.
To further confirm the ability of the proposed model, in figure[3], we present the comparison of the denoised Pirate
image which was initially degraded with additive Gaussian noise with very low SNR value 3.70. From figure[3], it
is clear that the proposed algorithm retains more texture information in addition to noise reduction, which makes our
filtered image better in comparison to other models considered here.
Along with the qualitative analysis using full image surface, we also study the quality of the resultant images considering a slice of the Boat image and the Livingroom image with different noise levels. Figure 4 shows the signals of the
original, noisy and restored images obtained by the proposed CPDE model and other discussed CPDE models. From
these figures, it is easy to conclude that the restored signals computed by the proposed model are more closure to the
clean signals in comparison to other discussed models.
In addition to qualitative comparisons, the quantitative results, in terms of PSNR and MSSIM values for different test
images as well as noise levels (σ = 20, 30, 40, 50), are shown in tables [1]-[4]. To ease the comparison, the highest
11
A
(a) Clear
PREPRINT
- AUGUST 8, 2019
(b) Noisy
(c) PV
(d) SYS
(e) Proposed
Figure 1: (a) Clear image; (b) Noisy image with Gaussian noise of mean 0.0 and σ = 40; Denoised image using (c)
PV Model, λ = 0.01, K = 8; (d) SYS Model, λ = 4, K = 8; (e) Proposed model, ψ = 0.1, k = 4.45.
Table 1: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is
degraded by additive Gaussian noise of mean 0.0 and σ = 20.
Images
Measure
TV[19]
NS[27]
Luo[26]
RD[25]
NLM[39]
PV[28]
SYS[38]
Proposed
Boat
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSI
PSNR
MSSIM
PSNR
MSSIM
PSNR
0.8232
28.10
0.8539
28.37
0.8790
28.51
0.8958
27.51
0.8615
28.54
0.9077
26.79
0.8007
30.35
0.9699
34.28
0.8883
28.81
0.8996
28.07
0.8698
28.62
0.9001
27.71
0.9053
29.40
0.9141
26.83
0.9191
33.83
0.9091
32.71
0.8774
29.40
0.8775
28.65
0.8855
28.68
0.8909
27.77
0.8764
29.20
0.9070
27.05
0.8257
31.31
0.9270
33.66
0.8886
29.61
0.8788
28.61
0.8849
28.80
0.8997
27.95
0.8714
29.07
0.9101
26.99
0.8206
31.29
0.9600
34.44
0.9105
30.21
0.9140
29.27
0.8870
29.11
0.8684
27.62
0.8908
29.69
0.8747
26.72
0.9243
33.64
0.9588
34.12
0.8900
29.69
0.9060
28.87
0.8903
28.87
0.9015
28.24
0.8936
29.51
0.9093
26.99
0.9262
34.47
0.9675
33.38
0.8943
29.84
0.9074
29.08
0.8912
28.93
0.9029
28.37
0.8943
29.53
0.9089
27.02
0.9225
34.59
0.9290
32.54
0.9219
30.24
0.9268
29.57
0.9019
29.22
0.9034
28.46
0.9055
29.90
0.9143
27.24
0.9396
35.03
0.9780
34.74
Lake
Livingroom
Mandril
Pirate
Walkbridge
Woman
Mosaic
Table 2: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is
degraded by additive Gaussian noise of mean 0.0 and σ = 30.
Images
Measure
TV[19]
NS[27]
Luo[26]
RD[25]
NLM[39]
PV[28]
SYS[38]
Proposed
Boat
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
0.7157
25.28
0.7644
25.27
0.7859
25.91
0.8188
24.94
0.7745
26.18
0.8391
24.60
0.6834
27.58
0.9558
31.63
0.6923
25.71
0.7676
25.41
0.8510
27.28
0.7633
25.46
0.8348
27.05
0.8194
24.07
0.8559
31.16
0.8827
30.38
0.7324
25.83
0.7751
25.72
0.7922
26.14
0.8262
25.55
0.7765
26.26
0.8431
24.91
0.6805
27.47
0.8800
30.68
0.7327
26.01
0.7820
26.01
0.7969
26.40
0.8268
25.52
0.7905
26.86
0.8453
24.99
0.7108
28.58
0.9467
31.65
0.8568
28.02
0.8705
27.38
0.8229
26.90
0.7873
25.43
0.8346
27.72
0.8015
24.96
0.8861
31.38
0.9477
31.69
0.8360
27.79
0.8576
27.12
0.8273
26.97
0.8309
26.13
0.8392
27.86
0.8462
25.19
0.8963
32.58
0.9520
31.10
0.8400
27.96
0.8637
27.33
0.8321
27.02
0.8395
26.27
0.8400
27.81
0.8508
25.22
0.8922
32.64
0.8764
30.24
0.8841
28.47
0.8956
27.87
0.8515
27.39
0.8454
26.47
0.8575
28.20
0.8566
25.47
0.9200
33.25
0.9659
31.97
Lake
Livingroom
Mandril
Pirate
Walkbridge
Woman
Mosaic
12
A
(a) Clear
(c) PV
PREPRINT
- AUGUST 8, 2019
(b) Noisy
(d) SYS
(e) Proposed
(f) Noisy
(g) PV
(h) SYS
(i) Proposed
Figure 2: (a) Clear image (b) Noisy image with Gaussian noise of mean 0.0 and s.d. 40 (c) PV Model, λ = 0.01, K =
1; (d) SYS Model, λ = 0.1, K = 4; (e) Proposed Model, ψ = 1, k = 3.75. (f-i) 3D surfaces of images
13
A
(a) Original
PREPRINT
- AUGUST 8, 2019
(b) Noisy
(c) PV Model
(d) SYS Model
(e) Proposed Model
Figure 3: (a) Clear image (b)Heavily noised image with low SNR=3.70; Improved SNR values by different models:
(c) SNR=15.41; λ = 0.01, K = 11 (d) SNR=15.37; λ = 1, K = 10 (e) SNR=15.53; ψ = 1, k = 14.
Table 3: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is
degraded by additive Gaussian noise of mean 0.0 and σ = 40.
Images
Measure
TV[19]
NS[27]
Luo[26]
RD[25]
NLM[39]
PV[28]
SYS[38]
Proposed
Boat
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
0.6305
23.43
0.6964
23.67
0.6941
23.85
0.7409
23.14
0.6882
24.32
0.7728
23.05
0.5976
25.75
0.9384
29.43
0.6146
23.41
0.7052
23.70
0.6989
23.36
0.7456
23.24
0.6913
24.23
0.7883
23.48
0.7317
27.59
0.8892
28.51
0.6290
23.35
0.6927
23.62
0.6855
23.56
0.7462
23.41
0.6800
24.04
0.7722
23.08
0.5671
25.93
0.8315
28.38
0.6481
24.19
0.7117
24.33
0.7095
24.50
0.7525
23.82
0.7067
25.06
0.7821
23.54
0.6235
26.61
0.9328
29.46
0.8022
26.40
0.8245
25.93
0.7600
25.37
0.7164
24.01
0.7849
26.35
0.7371
23.76
0.8513
29.75
0.9388
29.66
0.7862
26.51
0.8154
25.94
0.7699
25.68
0.7692
24.75
0.7869
26.64
0.7922
24.01
0.8690
31.00
0.9303
29.04
0.7798
26.60
0.8173
26.05
0.7731
25.70
0.7777
24.90
0.7870
26.58
0.7944
24.05
0.8641
30.96
0.8373
28.15
0.8470
27.21
0.8644
26.62
0.7991
26.11
0.7844
25.10
0.8132
27.06
0.7985
24.24
0.9037
31.65
0.9508
29.77
Lake
Livingroom
Mandril
Pirate
Walkbridge
Woman
Mosaic
Table 4: Comparison of MSSIM and PSNR values for various approaches and proposed model. Clean image is
degraded by additive Gaussian noise of mean 0.0 and σ = 50.
Images
Measure
TV [19]
NS [27]
Luo [26]
RD [25]
NLM[39]
PV[28]
SYS[38]
Proposed
Boat
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
MSSIM
PSNR
0.5648
22.02
0.6404
22.47
0.6265
22.51
0.6707
21.77
0.6275
23.11
0.7150
22.00
0.5318
24.30
0.9203
27.61
0.5793
21.70
0.6657
22.09
0.6424
21.57
0.6177
21.48
0.6486
21.66
0.7046
21.75
0.5702
22.02
0.8404
26.93
0.5519
21.50
0.6232
21.98
0.5981
21.54
0.6765
22.06
0.5969
22.15
0.7025
21.52
0.5626
22.56
0.8662
26.55
0.5813
22.75
0.6540
23.06
0.6429
23.17
0.6846
22.49
0.6450
23.78
0.7243
22.44
0.6143
25.02
0.9189
27.66
0.7523
25.14
0.7828
24.75
0.7090
24.28
0.6505
22.99
0.7423
25.24
0.6838
22.84
0.8150
28.23
0.9144
27.71
0.7425
25.48
0.7747
24.89
0.7245
24.75
0.7102
23.69
0.7483
25.68
0.7414
23.10
0.8430
29.31
0.9078
27.34
0.7333
25.46
0.7710
24.94
0.7229
24.72
0.7173
23.87
0.7454
25.60
0.7428
23.16
0.8354
29.25
0.7730
26.51
0.8159
26.16
0.8334
25.53
0.7562
25.17
0.7261
24.06
0.7762
26.10
0.7483
23.34
0.8830
29.81
0.9333
27.87
Lake
Livingroom
Mandril
Pirate
Walkbridge
Woman
Mosaic
14
A
300
PREPRINT
300
Noisy Signal
250
250
200
200
Intensity
Intensity
Noisy Signal
150
150
100
100
50
50
0
0
0
50
100
150
200
250
300
0
50
100
Pixel Coordinate
250
300
300
Clean signal
PV Model
250
Clean signal
PV Model
250
200
Intensity
Intensity
200
(b)
300
150
200
150
100
100
50
50
0
0
0
50
100
150
200
250
300
0
50
100
Pixel Coordinate
150
200
250
300
Pixel Coordinate
(c) λ = 0.1, K = 5
(d) λ = 0.1, K = 15
300
300
Clean signal
SYS Model
250
Clean signal
SYS Model
250
200
Intensity
Intensity
150
Pixel Coordinate
(a)
150
100
200
150
100
50
50
0
0
0
50
100
150
200
250
300
0
50
100
Pixel Coordinate
150
200
250
300
Pixel Coordinate
(e) λ = 2, K = 6
(f) λ = 4, K = 9
300
300
Clean Signal
Proposed Model
250
Clean Signal
Proposed Model
250
200
Intensity
Intensity
- AUGUST 8, 2019
150
100
200
150
100
50
50
0
0
0
50
100
150
200
250
300
Pixel Coordinate
0
50
100
150
200
250
300
Pixel Coordinate
(g) ψ = 1, k = 4
(h) ψ = 1, k = 4.45
Figure 4: (a) The 219th slice of Boat image corrupted with Gaussian noise of mean 0.0 and σ = 40; (b) The 95th
slice of Livingroom image corrupted with additive Gaussian noise of mean 0.0 and σ = 50; denoised using (c-d) PV
model; (e-f) SYS model, and (g-h) Proposed model.
15
Table 5: Comparison of ISNR and PSNRgrad values for various approaches and proposed model. Clean image is
degraded by additive Gaussian noise of mean 0.0 and σ = 40.
Images
Measure
TV[19]
NS[27]
Luo[26]
RD[25]
NLM[39]
PV[28]
SYS[38]
Proposed
Boat
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
ISNR
PSNRgrad
7.02
27.04
7.05
27.23
7.51
27.68
7.12
26.59
7.92
28.11
6.59
26.55
9.11
30.05
12.57
37.34
7.01
27.01
7.12
27.30
7.02
27.33
7.19
26.70
7.83
27.57
6.88
26.65
10.95
32.36
11.64
35.68
6.93
26.86
7.00
27.12
7.22
27.51
7.31
26.83
7.64
27.41
6.62
26.59
9.29
30.72
11.51
35.48
7.78
28.07
7.71
28.05
8.16
28.62
7.58
27.44
8.66
29.14
7.07
27.14
9.97
31.46
12.60
37.37
9.90
30.84
9.12
30.11
9.03
29.84
7.77
27.91
9.96
30.59
7.30
27.68
13.11
34.77
12.74
37.52
9.99
30.90
9.19
30.19
9.20
30.04
8.51
28.48
10.14
30.62
7.55
27.71
13.94
37.37
12.17
36.01
10.19
31.28
9.43
30.21
9.36
30.25
8.66
28.58
9.99
30.60
7.59
27.73
13.87
37.30
11.29
34.92
10.80
31.64
10.00
30.62
9.77
30.62
8.86
28.88
10.66
31.26
7.78
27.98
15.01
38.80
12.86
37.56
Lake
Livingroom
Mandril
Pirate
Walkbridge
Woman
Mosaic
values of both MSSIM and PSNR measures are highlighted in each table. The highest values of PSNR and MSSIM
for each image, clearly shows that the proposed model is better than all other models considered in the table. We
have a significant observation that proposed model performs better for higher noise densities, as can be seen from
tables[3]-[4]. For low noise densities, such as σ = 20, our results are still batter or very closed to the best values.
Hence, analysis of the results presented in the tables, reveals that the performance of the present model increases
with increasing noise density. In table [5], apart from the MSSIM and PSNR, we show the quantitative comparison
with gradient PSNR and ISNR values of the proposed model and the alternative approaches. Finally, we note that the
different quantitative measures (in terms of PSNR, MSSIM, gradient PSNR, and ISNR values) confirm the superiority
of the proposed approach in comparison with other models. Our comparison study shows that the proposed CPDE
model preserves image structure efficiently in comparison to other earlier reported methods when signal-to-noise ratio
is low. On the other hand, when signal-to-noise ratio is high, the present model restores image structure better than the
other existing models. From the above discussion, it is confirmed that the proposed model is robust and more efficient
than the available techniques considered here.
6 Conclusion
In this work, we propose a new CPDE based denoising framework. The proposed space-time regularization based
image denoising approach highlights the choice of diffusion function and data fidelity term. For this purpose, we use
the evolution equations to determine both the terms in the proposed model. It may be treated as the non-linear diffusion
model which yield two separate PDEs to remove the noise with preservation of significant edges and fine structures.
First, we establish the well-posedness of the proposed model using a time discretization method. Then, to solve the
proposed model numerically, an implicit finite difference scheme along with advanced iterative solver has been used.
We validate our model with different standard test images. Experimental studies confirm that the proposed model is
more efficient than the existing models, in terms of quantitative measures and visual quality. Especially, the present
model produces best results for high noise levels. We believe that present model has multidimensional applications
which will ultimately benefit the society.
Competing interests
The authors declare that they have no competing interests.
Acknowledgments
Authors sincerely acknowledge Dr. Arnav Bhavsar, School of Computing and Electrical Engineering, Indian Institute
of Technology Mandi, Himachal Pradesh, India for some valuable discussions and suggestions during the preparation
of this manuscript.
A
PREPRINT
- AUGUST 8, 2019
References
[1] L. Alvarez, F. Guichard, P.-L. Lions, J.-M. Morel, Axioms and fundamental equations of image processing,
Archive for rational mechanics and analysis 123 (3) (1993) 199-257.
[2] R. C. Gonzalez, R. E. Woods, Digital image processing (2002).
[3] G. Aubert, P. Kornprobst, Mathematical problems in image processing: partial differential equations and the
calculus of variations, Vol. 147, Springer, 2006.
[4] K. Mikula, Image processing with partial differential equations, in: Modern methods in scientifc computing and
applications, Springer, 2002, pp. 283-321.
[5] J. Nolen, Partial differential equations and difusion processes, Tech. rep., Technical report, Stanford University.
Department of Mathematics (2009).
[6] J. Weickert, Anisotropic diffusion in image processing, Vol. 1, Teubner Stuttgart, 1998.
[7] J. Weickert, Applications of nonlinear diffusion in image processing and computer vision, Acta Math. Univ. Comenianae 70 (1) (2001) 33-50.
[8] A. P. Witkin, Scale-space filtering: A new approach to multi-scale description,in: Acoustics, Speech, and Signal
Processing, IEEE International Conference on ICASSP’84., Vol. 9, IEEE, 1984, pp. 150-153.
[9] F. Catté, P.-L. Lions, J.-M. Morel, T. Coll, Image selective smoothing and edge detection by nonlinear diffusion,
SIAM Journal on Numerical analysis 29 (1) (1992) 182-193.
[10] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, Pattern Analysis and Machine
Intelligence, IEEE Transactions on 12 (7) (1990) 629-639.
[11] J. Weickert, A review of nonlinear diffusion filtering, in: Scale-space theory in computer vision, Springer, 1997,
pp. 1-28.
[12] X. Liu, L. Huang, Z. Guo, Adaptive fourth-order partial differential equation filter for image denoising, Applied
Mathematics Letters 24 (8) (2011) 1282-1288.
[13] A. Siddig, Z. Guo, Z. Zhou, B. Wu, An image denoising model based on a fourth-order nonlinear partial differential equation, Computers & Mathematics with Applications 76 (5) (2018) 1056-1074.
[14] Y.-L. You, M. Kaveh, Fourth-order partial differential equations for noise removal, IEEE Transactions on Image
Processing 9 (10) (2000) 1723-1730.
[15] A. Araújo, S. Barbeiro, P. Serranho, Stability of finite difference schemes for complex diffusion processes, SIAM
Journal on Numerical Analysis 50 (3) (2012) 1284-1296.
[16] A. Araújo, S. Barbeiro, P. Serranho, Stability of finite difference schemes for nonlinear complex reactiondiffusion processes, IMA Journal of Numerical Analysis 35 (3) (2015) 1381-1401.
[17] A. Chambolle, P.-L. Lions, Image recovery via total variation minimization and related problems, Numerische
Mathematik 76 (2) (1997) 167-188.
[18] C. Elliott, S. Smitheman, Numerical analysis of the tv regularization and H −1 fidelity model for decomposing
an image into cartoon plus texture, IMA Journal of Numerical Analysis 29 (3) (2009) 651-689.
[19] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear
Phenomena 60 (1) (1992) 259-268.
[20] V. S. Prasath, A. Singh, A hybrid convex variational model for image restoration, Applied Mathematics and
Computation 215 (10) (2010) 3655-3664.
[21] Y.-H. R. Tsai, S. Osher, Total variation and level set methods in image science, Acta Numerica 14 (2005) 509573.
[22] R. Zanella, F. Porta, V. Ruggiero, M. Zanetti, Serial and parallel approaches for image segmentation by numerical
minimization of a second-order functional, Applied Mathematics and Computation 318 (2018) 153-175.
[23] H. Amann, Time-delayed perona-malik type problems, Acta Math. Univ. Comenianae 76 (1) (2007) 15-38.
[24] A. Belahmidi, A. Chambolle, Time-delay regularization of anisotropic diffusion and image processing, ESAIM:
Mathematical Modelling and Numerical Analysis 39 (2) (2005) 231-251.
[25] Z. Guo, J. Yin, Q. Liu, On a reaction-diffusion system applied to image decomposition and restoration, Mathematical and Computer Modelling 53 (5) (2011) 1336-1350.
[26] H. Luo, L. Zhu, H. Ding, Coupled anisotropic diffusion for image selective smoothing, Signal Processing 86 (7)
(2006) 1728-1736.
17
A
PREPRINT
- AUGUST 8, 2019
[27] M. Nitzberg, T. Shiota, Nonlinear image filtering with edge and corner enhancement, IEEE transactions on
pattern analysis and machine intelligence 14 (8) (1992) 826-833.
[28] V. S. Prasath, D. Vorotnikov, On a system of adaptive coupled pdes for image restoration, Journal of mathematical
imaging and vision 48 (1) (2014) 35-52.
[29] P. Guidotti, Anisotropic diffusions of image processing from perona-malik on, Advanced Studies in Pure Mathematics 99 (2015) 20XX.
[30] S. Osher, A. Solé, L. Vese, Image decomposition and restoration using total variation minimization and the H −1
norm, Multiscale Modeling and Simulation 1 (3) (2003) 349-370.
[31] R. Adam, Sobolev spaces, in: Pure and Applied Mathematics Series of Monographs and Textbooks, Vol. 65,
Academic Press, Inc., New York, San Francisco, London„ 1975.
[32] L. Evans, Partial Differential Equations, in: Graduate Studies in Mathematics, Vol. 19, American Mathematical
Society, Providence, Rhode Island, 1998.
[33] S. K. Jain, R. K. Ray, A. Bhavsar, Iterative solvers for image denoising with diffusion models: A comparative
study, Elsevier, 2015.
[34] J. L. Lions, Contrôle Optimal de Systèmes Gouvernés par des Équations aux Dérivées Partielles, Dunod, Paris,
1968.
[35] M. Chipot, I. Shafrir, V. Valente, G. V. Caffarelli, On a hyperbolic- parabolic system arising in magnetoelasticity,
Journal of Mathematical Analysis and Applications 352 (1) (2009) 120-131.
[36] S. Zheng, Nonlinear parabolic equations and hyperbolic-parabolic coupled systems, CRC Press, 1995.
[37] S. K. Jain, R. K. Ray, A. Bhavsar, A comparative study of iterative solvers for image denoising, in: Proceedings
of the 3rd International Conference on Frontiers of Intelligent Computing: Theory and Applications (FICTA)
2014, Springer, 2015, pp. 307-314.
[38] J. Sun, J. Yang, L. Sun, A class of hyperbolic-parabolic coupled systems applied to image restoration, Boundary
Value Problems 2016 (1) (2016) 187.
[39] A. Buades, B. Coll, J.-M. Morel, A non-local algorithm for image denoising, in: Computer Vision and Pattern
Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, Vol. 2, IEEE, 2005, pp. 60-65.
[40] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, Image quality assessment: from error visibility to structural
similarity, Image Processing, IEEE Transactions on 13 (4) (2004) 600-612.
[41] X.-C. Tai, K.-A. Lie, T. F. Chan, S. Osher, Image Processing Based on Partial Differential Equations: Proceedings of the International Conference on PDE-Based Image Processing and Related Inverse Problems, CMA,Oslo,
August 8-12, 2005, Springer Science & Business Media, 2006.
18