# A Total Variation Model Based on the Strictly Convex Modification for Image Denoising

Authors: Boying Wu; Elisha Achieng Ogada; Jiebao Sun; Zhichang Guo
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a strictly convex functional in which the regular term consists of the total variation term and an adaptive logarithm based convex modification term. We prove the existence and uniqueness of the minimizer for the proposed variational problem. The existence, uniqueness, and long-time behavior of the solution of the associated evolution system is also established. Finally, we present experimental results to illustrate the effectiveness of the model in noise reduction, and a comparison is made in relation to the more classical methods of the traditional total variation (TV), the Perona-Malik (PM), and the more recent D-α-PM method. Additional distinction from the other methods is that the parameters, for manual manipulation, in the proposed algorithm are reduced to basically only one.


---

## Body

## 1. Introduction

Noise removal, edge detection, contrast enhancement, inpainting, and segmentation have been the subject of intense mathematical image analysis and processing research for nearly three decades. Several methods have been pursued over the passage of time. These include wavelet transform [1, 2], curvelet shrinkage methods [3–6], and variational partial differential equation (PDE) based methods [7, 8]. These methods generate processes that can easily be divided into either linear and nonlinear processes or isotropic and nonisotropic processes [9].Due its ability to preserve crucial image features, such as edges, nonlinear anisotropic diffusion is favored over isotropic diffusion [10]. Much interest, therefore, has focused on understanding operations and mathematical properties of the nonlinear anisotropic diffusion and associated variational formulations [11, 12], formulation of well-posed and stable equations [8, 11], extending and modifying anisotropic diffusion [11, 13, 14], and studying the relationships that exist between the various image processing techniques [11, 15, 16].The objective of any image denoising process should not focus only on the removal of noise, but it should also ensure that no spurious details are created on the restored image and that the edges are preserved or sharpened [7, 17, 18]. It is, therefore, necessary to develop formulations which are sensitive to the local image structure, especially edges/contours [19]. Consequently, a number of edge indicators have been proposed and logically grafted into the partial differential equation (PDE) based evolution equations [7, 11, 20].Some of these PDEs originate from variational problems. For instance, Rudin et al. [8] proposed a minimization functional, widely referred to as the total variation (TV) functional, of the form
(1)
min
⁡
u
{
F
(
u
)
=
∫
Ω
‍
|
∇
u
|
d
x
+
λ
∫
Ω
‍
(
u
-
f
)
2
d
x
}
,
where λ is the fidelity parameter, f
=
f
(
x
) denotes the noise image, and Ω is an open bounded subset of R
2. TV functionals are defined in the space of functions of bounded variation (BV) and, therefore, do not necessarily require image functions to be continuous and smooth. This fact makes them allow jumps or discontinuities, and hence they are able to preserve edges.The original TV formulation has certain weaknesses. Firstly, the formulation is susceptible to backward diffusion since it is not strictly convex. Secondly, the numerical implementation cannot be accomplished without additional small perturbation quantity, sayϵ, at the denominator [21–23]. Otherwise a spike/singularity is suddenly generated when |
∇
u
|
=
0 in the homogeneous regions. This perturbation phenomenon is believed to contribute to some loss of accuracy in the results of the restoration process. Moreover, the additional parameter unnecessarily increases the number of parameters, thereby making it difficult to determine which permutation of parameter values will give optimal result.Additionally, given that the method is more efficient in preserving edges of uniform and small curvature, it may excessively smoothen and possibly destroy small scale features having more pronounced curvature edges [24]. TV regularization approach may also result in a loss of contrast and geometry of the final output images, even in noise free observed images [9, 25]. Furthermore, TV regularization has difficulties recovering texture, and there is also evidence of enhanced noise when the fidelity parameter is chosen so that texture is not removed [26]. Lastly, the formulation favors piecewise constant solutions. This has the material effect of causing staircases (false edges) on the resultant image, especially from images severely degraded by noise [25].However, given the strength of TV based techniques, especially in edge preservation, various modifications have been proposed. For instance, Vogel in [22, 27] proposed a total variation penalty method of the form
(2)
min
⁡
u
∫
Ω
(
A
u
-
z
)
2
d
x
‍
+
α
∫
Ω
‍
(
|
∇
u
|
2
+
β
)
d
x
,
where β
≥
0, A is a linear operator, and α
>
0 is the penalty parameter. The formulation becomes total variation formulation of the Rudin et al. form [8] when β
=
0 and therefore largely suffers the same shortcoming of the ordinary TV formulation. Hence, it has no significant practical advantage over TV.Strong and Chan [28] proposed an adaptive total variation based regularization model of the form
(3)
min
⁡
u
∈
BV
(
Ω
)
∫
Ω
‍
α
(
x
)
|
∇
u
|
,
where 0
≤
α
(
x
)
≤
1 is a control factor which controls the speed of diffusion depending on whether the region is homogeneous or an edge. This model demonstrated fairly good results. However, since it is also not strictly convex, it is still susceptible to backward diffusion, which has the potential of introducing blurs in the restored image.Chambolle and Lions [29] proposed to minimize a combination of total variation and the integral of the squared norm of the gradient and thus have
(4)
min
⁡
u
∈
BV
(
Ω
)
1
2
ε
∫
|
∇
u
|
≤
ε
|
∇
u
|
2
+
∫
|
∇
u
|
≥
ε
|
∇
u
|
-
ε
2
‍
,
where ε is a parameter. In this formulation |
∇
u
|
≥
ε signals edges, while |
∇
u
|
≤
ε signals homogeneous regions. This model is successful in restoring images where homogeneous regions are separated by distinct edges but may become sensitive to the thresholding parameter ε in the event of nonuniform image intensities or heavy degradation [24].A variable exponent adaptive model was proposed by Chen et al. in [24]. It has the form
(5)
min
⁡
u
∈
BV
(
Ω
)
∫
Ω
‍
|
D
u
|
q
(
x
)
+
λ
2
(
u
-
I
)
2
,
λ
>
0
,
where 1
<
α
≤
q
(
x
)
≤
2, with q
(
x
) proposed as q
(
x
)
=
1
+
1
/
(
1
+
k
|
∇
G
σ
*
I
(
x
)
|
2
), G
σ is the Gaussian filter, and k
>
0, σ
>
0 are fixed parameters. It is observed that the method exploits the benefits of Gaussian smoothing when q
(
x
)
=
2 and the strength of TV regularization when q
(
x
)
=
1. This method also demonstrates good results and is indeed a great improvement of the earlier models. However, the convolution of the image with a Gaussian kernel before the evolution process diminishes the accuracy of these models. This is because of the introduction of the scale variance of the Gaussian; the scale variance is itself an additional parameter that is subject to manual manipulation [12]. Furthermore, the diffusion process becomes ill-posed if scale variance is too small, while image features become smeared if the Gaussian variance is too large [30]. Therefore, optimal selection of the scale variance remains a challenge. Parameter permutation giving optimal result is a challenge due to the number of such parameters involved.Further literature surveys attest to the fact that research in effective regularization functionals which have the ability to generate diffusion processes that restore images, while simultaneously preserving critical images features, the analysis, and practical implementation of such models, is still an extremely active concern.Consequently, in this paper, we propose a new adaptive total variation (TV) formulation for image denoising, which is strictly convex. The only parameterK, which is a thresholding parameter to be tweaked, is such that it depends on the evolution parameter t and is therefore not as grossly limited as the choice of numerical values. The other parameter λ is dynamically obtained and is therefore not manually tuned. With all this the number of parameters is reduced to basically only one. Our method approximates TV model, for higher values of the thresholding parameter K. Experimental results, presented here, demonstrate the effectiveness of our model over the classical models of TV, PM, and D-α-PM.The structure of this paper is as follows. In Section2, we present the proposed model (6), its properties, justifications of the model, and the associated evolution equation. In Section 3, we give certain preliminary definitions we rely on from time to time in this paper; we also prove the existence and uniqueness of the solution to the minimization problem (6). In Section 4, we discuss the associated evolution equation to the minimization problem (6). Also, we define the weak solution to the evolution problem (12)–(14), derive the estimates for the solution of an approximating problem (32), prove the existence and uniqueness of the weak solution of the evolution problem (12)–(14), and discuss the asymptotic behaviour of the weak solution as t
→
∞. In Section 5, we give the numerical schemes and experimental results to demonstrate the strength and effectiveness of our method. We have also presented a brief comparative discussion of our results with the PM, original TV methods, and D-α-PM methods. A brief summary concludes the paper in Section 6.


## 2. Proposed Model

In this section, motivated by the methods of [11, 21, 24], we propose a modified version of total variation denoising model. The model is based on minimization of a total variation functional with a logarithm-based strictly convex modification. First, we present the model and certain important properties of it.


### 2.1. The New Energy Functional

The new strictly convex energy functional is given as follows:(6)
min
⁡
u
∈
BV
(
Ω
)
∩
L
2
(
Ω
)
{
I
(
u
)
=
∫
Ω
‍
Φ
(
|
∇
u
|
)
d
x
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
}
,
where
(7)
Φ
(
s
)
=
s
-
1
K
ln
⁡
(
1
+
K
s
)
,
where K is a positive parameter, and f is the noise image.Calculating directly, we can obtain the following proposition.Proposition 1.
The functionΦ
(
s
) satisfies the following properties: (i)
Φ
(
0
)
=
0
,
Φ
′
(
s
)
=
K
s
/
(
1
+
K
s
)
, and therefore
Φ
(
s
)
is a strictly nondecreasing function for
s
>
0;
(ii)
Φ
′′
(
s
)
=
K
/
(
1
+
K
s
)
2
>
0
, and therefore
Φ
(
s
)
is strictly convex;
(iii)
Φ
(
s
)
=
(
K
/
2
)
s
2
+
o
(
s
2
)
as
s
→
0
;
(iv)
lim
⁡
s
→
+
∞
⁡
(
Φ
(
s
)
/
s
)
=
1
,
α
|
s
|
≤
Φ
(
|
s
|
)
<
|
s
|
,
0
<
α
<
1
, and therefore
Φ
(
s
)
is the linear growth.Compared with the work in [31], the new model satisfies the linear growth condition, which is much more natural. And the existence result for the flow associated with the minimization problem in [31] is only in one and two dimensions because the methods employed there use general results on maximal monotone operators and evolution operators in Hilbert spaces.Remark 2.
In this method,Φ satisfies some minimal hypotheses; namely: (H1)
Φ is a strictly convex, nondecreasing function from R
+ to R
+, with Φ
(
0
)
=
0;
(H2)
because in the homogeneous regions weak gradients should be smoothed, there should be weak penalization in these region. Thus as|
∇
u
|
=
s
→
0, Φ
(
s
)
≈
K
s
2; this signals isotropic diffusion;
(H3)
because edges, signaled by regions of strong gradients, should be protected, regularization near the edges will be penalized strongly. Thus, as|
∇
u
|
=
s
→
+
∞, α
|
s
|
<
Φ
(
|
s
|
)
<
Λ
|
s
|, 0
<
α
≤
Λ; this not only demonstrates the linear growth nature of the model but also signals the TV edge preservation behaviour of the model.In [27], the authors obtained
(8)
lim
⁡
β
→
0
J
β
(
u
)
=
J
0
(
u
)
,
where J
β
(
u
)
=
∫
Ω
‍
β
+
|
∇
u
|
2
d
x. Actually, the perturbation β is always used to the eliminate singularity of the term div
⁡
(
D
u
/
|
D
u
|
) in the numerical experiments. With the inclusion or addition of β, what is then implemented is not actually the TV formulation, but an approximation of it. And, so, because of the approximation, the time step in the discrete scheme becomes limited only to smaller values, as β diminishes [32]. The proposed model beats all these challenges and therefore has an easier design for numerical implementation without the need for lifting parameters. Moreover, it is noticed that, for s
∈
R,
(9)
lim
⁡
K
→
+
∞
Φ
(
|
s
|
)
=
|
s
|
,
which is the form of the TV model. This means that any merits accruing from the TV model could be still obtained as K becomes larger and larger.Remark 3.
(1) Since TV potential (i.e.,Φ
(
s
)
=
s) is only convex (i.e., Φ
′′
(
s
)
=
0), it gives a local minimum, but it cannot guarantee uniqueness of the minimum energy. The proposed energy potential (i.e., Φ
(
s
)
=
s
-
(
1
/
K
)
ln
⁡
(
1
+
K
s
)), however, is strictly convex (i.e., Φ
′′
(
s
)
=
K
/
(
1
+
K
s
)
2
>
0). And, with additional strict convexity in the fidelity part, H
(
z
)
=
(
z
-
f
)
2, the resultant energy functional is strictly convex in both |
∇
u
| and u, thereby giving a global minimum energy and hence guaranteeing uniqueness of the results.
(2) The proposed method seeks to reduce the number of parameters subject to manual manipulation. In the implementation,K is made to depend on time, thereby leaving time (evolution parameter) as the only parameter to be tweaked.Other modifications of the TV model are only convex in|
∇
u
| and therefore do not guarantee uniqueness of solutions; hence uniqueness of results is not easy to assure. Besides, the evolution system of the models contains 1
/
|
∇
u
| component, which runs into a spike when |
∇
u
|
=
0, in smooth regions. This then makes only the implementation of the approximations of the corresponding formulations possible (see [13, 20, 24, 27, 29]).


### 2.2. The Associated Evolution Equation

The Euler-Lagrange equation for the energy functional (6) is
(10)
0
=
-
div
⁡
(
K
∇
u
1
+
K
|
∇
u
|
)
+
λ
(
u
-
f
)
,
x
∈
Ω
,
with
(11)
∂
u
∂
n
→
=
0
,
x
∈
∂
Ω
.
To compute a solution of (10) numerically, it is embedded into a dynamical scheme, where t is used as the evolution parameter. Hence, corresponding to the proposed minimization model (6), we have the following evolution system:
(12)
u
t
=
div
⁡
(
K
∇
u
1
+
K
|
∇
u
|
)
-
λ
(
u
-
f
)
,
(
x
,
t
)
∈
Ω
×
(
0
,
T
)
,
(13)
u
(
x
,
0
)
=
f
(
x
)
,
x
∈
Ω
,
(14)
∂
u
∂
n
→
=
0
,
(
x
,
t
)
∈
∂
Ω
×
[
0
,
T
]
.Note 1.
The modified formulation also does not have the disadvantage of running into a singularity, since the denominator does not become zero. This particular observation makes it more convenient to design a numerical scheme for the model.Furthermore, to see more clearly the action of the diffusion operator (kernel), we decompose the divergence term on the premise of the local image structures. That is, we break it into the tangential (T) and normal (N) directions to the isophote lines. Consequently we have
(15)
div
⁡
(
∇
u
1
+
K
|
∇
u
|
)
=
1
(
1
+
K
|
∇
u
|
)
2
u
N
N
+
1
1
+
K
|
∇
u
|
u
T
T
,
where
(16)
u
N
N
=
1
|
∇
u
|
2
(
u
x
1
2
u
x
1
x
1
+
u
x
2
2
u
x
2
x
2
+
2
u
x
1
u
x
2
u
x
1
x
2
)
,
u
T
T
=
1
|
∇
u
|
2
(
u
x
1
2
u
x
2
x
2
+
u
x
2
2
u
x
1
x
1
-
2
u
x
1
u
x
2
u
x
1
x
2
)
.
The divergence term is visibly a weighted sum of the directional derivatives, along the normal and tangential directions. Since we are also seeking to preserve the edges and other important features of the image, smoothing in the tangential direction should be encouraged more than that in the normal direction, in the regions near the boundaries (edges). Observe in the above formulation that, as |
∇
u
| increases, the coefficient of u
N
N diminishes faster, reducing diffusion in the normal direction, thus preserving edges. The edges are signalled by higher values of the magnitude of gradient of u. However, as |
∇
u
| decreases, there is relatively uniform diffusion in both u
N
N and u
T
T directions, thereby achieving isotropic diffusion in homogeneous regions. The homogeneous regions are signaled by diminishing values of the magnitude of gradient of image u. The proposed model is, therefore, sensitive to the local image structure. However, TV based denoising model and most other modifications do not have that u
N
N component. This leads to a situation where the homogeneous parts are processed into piecewise constant regions, whose boundaries reflect staircases or false edges in the image. This situation, for instance, affects the results of a modification by Chen et al. [24] when p
=
1, and the model becomes the traditional TV model.


## 2.1. The New Energy Functional

The new strictly convex energy functional is given as follows:(6)
min
⁡
u
∈
BV
(
Ω
)
∩
L
2
(
Ω
)
{
I
(
u
)
=
∫
Ω
‍
Φ
(
|
∇
u
|
)
d
x
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
}
,
where
(7)
Φ
(
s
)
=
s
-
1
K
ln
⁡
(
1
+
K
s
)
,
where K is a positive parameter, and f is the noise image.Calculating directly, we can obtain the following proposition.Proposition 1.
The functionΦ
(
s
) satisfies the following properties: (i)
Φ
(
0
)
=
0
,
Φ
′
(
s
)
=
K
s
/
(
1
+
K
s
)
, and therefore
Φ
(
s
)
is a strictly nondecreasing function for
s
>
0;
(ii)
Φ
′′
(
s
)
=
K
/
(
1
+
K
s
)
2
>
0
, and therefore
Φ
(
s
)
is strictly convex;
(iii)
Φ
(
s
)
=
(
K
/
2
)
s
2
+
o
(
s
2
)
as
s
→
0
;
(iv)
lim
⁡
s
→
+
∞
⁡
(
Φ
(
s
)
/
s
)
=
1
,
α
|
s
|
≤
Φ
(
|
s
|
)
<
|
s
|
,
0
<
α
<
1
, and therefore
Φ
(
s
)
is the linear growth.Compared with the work in [31], the new model satisfies the linear growth condition, which is much more natural. And the existence result for the flow associated with the minimization problem in [31] is only in one and two dimensions because the methods employed there use general results on maximal monotone operators and evolution operators in Hilbert spaces.Remark 2.
In this method,Φ satisfies some minimal hypotheses; namely: (H1)
Φ is a strictly convex, nondecreasing function from R
+ to R
+, with Φ
(
0
)
=
0;
(H2)
because in the homogeneous regions weak gradients should be smoothed, there should be weak penalization in these region. Thus as|
∇
u
|
=
s
→
0, Φ
(
s
)
≈
K
s
2; this signals isotropic diffusion;
(H3)
because edges, signaled by regions of strong gradients, should be protected, regularization near the edges will be penalized strongly. Thus, as|
∇
u
|
=
s
→
+
∞, α
|
s
|
<
Φ
(
|
s
|
)
<
Λ
|
s
|, 0
<
α
≤
Λ; this not only demonstrates the linear growth nature of the model but also signals the TV edge preservation behaviour of the model.In [27], the authors obtained
(8)
lim
⁡
β
→
0
J
β
(
u
)
=
J
0
(
u
)
,
where J
β
(
u
)
=
∫
Ω
‍
β
+
|
∇
u
|
2
d
x. Actually, the perturbation β is always used to the eliminate singularity of the term div
⁡
(
D
u
/
|
D
u
|
) in the numerical experiments. With the inclusion or addition of β, what is then implemented is not actually the TV formulation, but an approximation of it. And, so, because of the approximation, the time step in the discrete scheme becomes limited only to smaller values, as β diminishes [32]. The proposed model beats all these challenges and therefore has an easier design for numerical implementation without the need for lifting parameters. Moreover, it is noticed that, for s
∈
R,
(9)
lim
⁡
K
→
+
∞
Φ
(
|
s
|
)
=
|
s
|
,
which is the form of the TV model. This means that any merits accruing from the TV model could be still obtained as K becomes larger and larger.Remark 3.
(1) Since TV potential (i.e.,Φ
(
s
)
=
s) is only convex (i.e., Φ
′′
(
s
)
=
0), it gives a local minimum, but it cannot guarantee uniqueness of the minimum energy. The proposed energy potential (i.e., Φ
(
s
)
=
s
-
(
1
/
K
)
ln
⁡
(
1
+
K
s
)), however, is strictly convex (i.e., Φ
′′
(
s
)
=
K
/
(
1
+
K
s
)
2
>
0). And, with additional strict convexity in the fidelity part, H
(
z
)
=
(
z
-
f
)
2, the resultant energy functional is strictly convex in both |
∇
u
| and u, thereby giving a global minimum energy and hence guaranteeing uniqueness of the results.
(2) The proposed method seeks to reduce the number of parameters subject to manual manipulation. In the implementation,K is made to depend on time, thereby leaving time (evolution parameter) as the only parameter to be tweaked.Other modifications of the TV model are only convex in|
∇
u
| and therefore do not guarantee uniqueness of solutions; hence uniqueness of results is not easy to assure. Besides, the evolution system of the models contains 1
/
|
∇
u
| component, which runs into a spike when |
∇
u
|
=
0, in smooth regions. This then makes only the implementation of the approximations of the corresponding formulations possible (see [13, 20, 24, 27, 29]).


## 2.2. The Associated Evolution Equation

The Euler-Lagrange equation for the energy functional (6) is
(10)
0
=
-
div
⁡
(
K
∇
u
1
+
K
|
∇
u
|
)
+
λ
(
u
-
f
)
,
x
∈
Ω
,
with
(11)
∂
u
∂
n
→
=
0
,
x
∈
∂
Ω
.
To compute a solution of (10) numerically, it is embedded into a dynamical scheme, where t is used as the evolution parameter. Hence, corresponding to the proposed minimization model (6), we have the following evolution system:
(12)
u
t
=
div
⁡
(
K
∇
u
1
+
K
|
∇
u
|
)
-
λ
(
u
-
f
)
,
(
x
,
t
)
∈
Ω
×
(
0
,
T
)
,
(13)
u
(
x
,
0
)
=
f
(
x
)
,
x
∈
Ω
,
(14)
∂
u
∂
n
→
=
0
,
(
x
,
t
)
∈
∂
Ω
×
[
0
,
T
]
.Note 1.
The modified formulation also does not have the disadvantage of running into a singularity, since the denominator does not become zero. This particular observation makes it more convenient to design a numerical scheme for the model.Furthermore, to see more clearly the action of the diffusion operator (kernel), we decompose the divergence term on the premise of the local image structures. That is, we break it into the tangential (T) and normal (N) directions to the isophote lines. Consequently we have
(15)
div
⁡
(
∇
u
1
+
K
|
∇
u
|
)
=
1
(
1
+
K
|
∇
u
|
)
2
u
N
N
+
1
1
+
K
|
∇
u
|
u
T
T
,
where
(16)
u
N
N
=
1
|
∇
u
|
2
(
u
x
1
2
u
x
1
x
1
+
u
x
2
2
u
x
2
x
2
+
2
u
x
1
u
x
2
u
x
1
x
2
)
,
u
T
T
=
1
|
∇
u
|
2
(
u
x
1
2
u
x
2
x
2
+
u
x
2
2
u
x
1
x
1
-
2
u
x
1
u
x
2
u
x
1
x
2
)
.
The divergence term is visibly a weighted sum of the directional derivatives, along the normal and tangential directions. Since we are also seeking to preserve the edges and other important features of the image, smoothing in the tangential direction should be encouraged more than that in the normal direction, in the regions near the boundaries (edges). Observe in the above formulation that, as |
∇
u
| increases, the coefficient of u
N
N diminishes faster, reducing diffusion in the normal direction, thus preserving edges. The edges are signalled by higher values of the magnitude of gradient of u. However, as |
∇
u
| decreases, there is relatively uniform diffusion in both u
N
N and u
T
T directions, thereby achieving isotropic diffusion in homogeneous regions. The homogeneous regions are signaled by diminishing values of the magnitude of gradient of image u. The proposed model is, therefore, sensitive to the local image structure. However, TV based denoising model and most other modifications do not have that u
N
N component. This leads to a situation where the homogeneous parts are processed into piecewise constant regions, whose boundaries reflect staircases or false edges in the image. This situation, for instance, affects the results of a modification by Chen et al. [24] when p
=
1, and the model becomes the traditional TV model.


## 3. The Minimization Problem

In this section, we prove the existence and uniqueness of the minimization problem (6). But, first, we give the following preliminary presentations which will guide our reasoning in the subsequent sections of this paper.


### 3.1. Preliminaries

Definition 4.
LetΩ be an open subset of R
n. A function u
∈
L
1 has a bounded variation in Ω if
(17)
sup
⁡
{
∫
Ω
‍
u
div
⁡
φ
d
x
:
φ
∈
C
0
1
(
Ω
;
R
n
)
,
|
φ
|
≤
1
}
<
∞
,
where BV
(
Ω
) denotes the space of such functions. Then BV-norm is given by
(18)
∥
∇
u
∥
BV
=
∫
Ω
‍
|
∇
u
|
+
|
u
|
L
1
(
Ω
)
.Definition 5 (see [33]).
Ifu
∈
BV
(
Ω
), then
(19)
D
u
=
∇
u
·
L
n
+
D
s
u
and D
u is a radon measure, where ∇
u is the density of the absolutely continuous part of D
u with respect to the Lebesgue measure, L
n, and D
s
u is the singular part.Lemma 6 (see [34]).
Letϕ
:
R
→
R
+ be convex, even, nondecreasing on R
+ with linear growth at infinity. Also let Φ
∞ be the recession function of Φ defined by
(20)
Φ
∞
(
ω
)
=
lim
⁡
s
→
∞
Φ
(
s
ω
)
s
.
Then for u
∈
B
V
(
Ω
) and setting Φ
(
θ
)
=
Φ
(
|
θ
|
) we have
(21)
∫
Ω
‍
Φ
(
D
u
)
=
∫
Ω
‍
Φ
(
|
∇
u
|
)
d
x
+
Φ
∞
(
1
)
∫
Ω
‍
D
s
u
.
This implies that u
→
∫
Ω
‍
Φ
(
D
u
) is lower semicontinuous for the B
V
(
Ω
) weak* topology.


### 3.2. Existence and Uniqueness of Solution to the Minimization Problem

The linear growth condition makes it natural to consider a solution in the space(22)
U
=
{
u
∈
L
2
(
Ω
)
;
∇
u
∈
L
1
(
Ω
)
2
}
.
However, this space is not reflexive. But sequences bounded in U are also bounded in BV
(
Ω
) and are therefore compact for BV weak* topology. Moreover, due to the fact that L
1-space is separable, weak* topology allows us to obtain compactness results even if the space is not reflexive [34–37]. Hence, by denoting BV
(
Ω
) weak* topology simply as BV
(
Ω
), we therefore seek the solution to the minimization problem (6) in the space BV
(
Ω
)
∩
L
2
(
Ω
).Theorem 7.
The minimization problemI
(
u
) (refer to (6)) has a unique solution u
∈
B
V
(
Ω
)
∩
L
2
(
Ω
).Proof .
From the linear growth property ofΦ
(
s
) we have
(23)
α
|
s
|
≤
Φ
(
|
s
|
)
<
|
s
|
,
for
0
<
α
<
1
,
which implies that
(24)
lim
⁡
s
→
+
∞
Φ
(
s
)
=
+
∞
.
Thus I
(
u
) is coercive. Therefore, let u
n be a minimizing sequence in BV
(
Ω
)
∩
L
2
(
Ω
). Then
(25)
I
(
u
n
)
≤
M
,
where M denotes a generic constant that may differ from line to line. It clearly follows from above that
(26)
∫
Ω
‍
Φ
(
|
∇
u
n
|
)
d
x
≤
I
(
u
n
)
≤
M
,
∫
Ω
‍
|
u
n
|
2
d
x
≤
M
.
This implies that {
u
n
} is bounded in BV
(
Ω
)
∩
L
2
(
Ω
). Hence there exists a subsequence {
u
n
k
} of {
u
n
} and a function u
∈
BV
(
Ω
)
∩
L
2
(
Ω
) such that
(27)
u
n
k
⟶
u
,
strongly
in
L
1
(
Ω
)
,
u
n
k
⇀
u
,
weakly
in
L
2
(
Ω
)
.
And, by Lemma 6 and the weak lower semicontinuity of the L
2-norm, we have that
(28)
I
(
u
)
≤
lim
⁡
inf
⁡
k
→
∞
I
(
u
n
k
)
=
min
⁡
BV
(
Ω
)
∩
L
2
(
Ω
)
I
(
v
)
.
Hence, there exists a solution of the minimization problem. Uniqueness of the solution follows from the strict convexity of the functional.


## 3.1. Preliminaries

Definition 4.
LetΩ be an open subset of R
n. A function u
∈
L
1 has a bounded variation in Ω if
(17)
sup
⁡
{
∫
Ω
‍
u
div
⁡
φ
d
x
:
φ
∈
C
0
1
(
Ω
;
R
n
)
,
|
φ
|
≤
1
}
<
∞
,
where BV
(
Ω
) denotes the space of such functions. Then BV-norm is given by
(18)
∥
∇
u
∥
BV
=
∫
Ω
‍
|
∇
u
|
+
|
u
|
L
1
(
Ω
)
.Definition 5 (see [33]).
Ifu
∈
BV
(
Ω
), then
(19)
D
u
=
∇
u
·
L
n
+
D
s
u
and D
u is a radon measure, where ∇
u is the density of the absolutely continuous part of D
u with respect to the Lebesgue measure, L
n, and D
s
u is the singular part.Lemma 6 (see [34]).
Letϕ
:
R
→
R
+ be convex, even, nondecreasing on R
+ with linear growth at infinity. Also let Φ
∞ be the recession function of Φ defined by
(20)
Φ
∞
(
ω
)
=
lim
⁡
s
→
∞
Φ
(
s
ω
)
s
.
Then for u
∈
B
V
(
Ω
) and setting Φ
(
θ
)
=
Φ
(
|
θ
|
) we have
(21)
∫
Ω
‍
Φ
(
D
u
)
=
∫
Ω
‍
Φ
(
|
∇
u
|
)
d
x
+
Φ
∞
(
1
)
∫
Ω
‍
D
s
u
.
This implies that u
→
∫
Ω
‍
Φ
(
D
u
) is lower semicontinuous for the B
V
(
Ω
) weak* topology.


## 3.2. Existence and Uniqueness of Solution to the Minimization Problem

The linear growth condition makes it natural to consider a solution in the space(22)
U
=
{
u
∈
L
2
(
Ω
)
;
∇
u
∈
L
1
(
Ω
)
2
}
.
However, this space is not reflexive. But sequences bounded in U are also bounded in BV
(
Ω
) and are therefore compact for BV weak* topology. Moreover, due to the fact that L
1-space is separable, weak* topology allows us to obtain compactness results even if the space is not reflexive [34–37]. Hence, by denoting BV
(
Ω
) weak* topology simply as BV
(
Ω
), we therefore seek the solution to the minimization problem (6) in the space BV
(
Ω
)
∩
L
2
(
Ω
).Theorem 7.
The minimization problemI
(
u
) (refer to (6)) has a unique solution u
∈
B
V
(
Ω
)
∩
L
2
(
Ω
).Proof .
From the linear growth property ofΦ
(
s
) we have
(23)
α
|
s
|
≤
Φ
(
|
s
|
)
<
|
s
|
,
for
0
<
α
<
1
,
which implies that
(24)
lim
⁡
s
→
+
∞
Φ
(
s
)
=
+
∞
.
Thus I
(
u
) is coercive. Therefore, let u
n be a minimizing sequence in BV
(
Ω
)
∩
L
2
(
Ω
). Then
(25)
I
(
u
n
)
≤
M
,
where M denotes a generic constant that may differ from line to line. It clearly follows from above that
(26)
∫
Ω
‍
Φ
(
|
∇
u
n
|
)
d
x
≤
I
(
u
n
)
≤
M
,
∫
Ω
‍
|
u
n
|
2
d
x
≤
M
.
This implies that {
u
n
} is bounded in BV
(
Ω
)
∩
L
2
(
Ω
). Hence there exists a subsequence {
u
n
k
} of {
u
n
} and a function u
∈
BV
(
Ω
)
∩
L
2
(
Ω
) such that
(27)
u
n
k
⟶
u
,
strongly
in
L
1
(
Ω
)
,
u
n
k
⇀
u
,
weakly
in
L
2
(
Ω
)
.
And, by Lemma 6 and the weak lower semicontinuity of the L
2-norm, we have that
(28)
I
(
u
)
≤
lim
⁡
inf
⁡
k
→
∞
I
(
u
n
k
)
=
min
⁡
BV
(
Ω
)
∩
L
2
(
Ω
)
I
(
v
)
.
Hence, there exists a solution of the minimization problem. Uniqueness of the solution follows from the strict convexity of the functional.


## 4. Existence and Uniqueness for the Evolution Equation

In this section, we define a weak solution for the evolution equations (12)–(14). That is, if (6) has a minimum point u, then it should formally satisfy the Euler-Lagrange equation (12), subject to the boundary conditions (13)-(14). We define a weak solution that we are seeking for the equation system (12)–(14). We then propose an approximating problem and derive certaina priori estimates to the approximating problem. These will aid the process of proving the existence and uniqueness of the weak solution. Finally, we show the large time behavior (asymptotic stability) of the weak solution.


### 4.1. Definition of Weak Solution

Letv
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
), f
∈
BV
(
Ω
), v
>
0, and ∂
v
/
∂
n
→
=
0, where Ω
T
=
Ω
×
[
0
,
T
] and u is a solution to (12)–(14). Multiplying (12) by (
v
-
u
) and integrating over Ω, we have
(29)
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
+
∫
Ω
‍
K
∇
u
1
+
K
|
∇
u
|
(
∇
v
-
∇
u
)
d
x
=
-
λ
∫
Ω
‍
(
u
-
f
)
(
v
-
u
)
d
x
.
Then applying the convexity condition of Φ, that is, Φ
(
x
)
-
Φ
(
y
)
≥
Φ
′
(
y
)
(
x
-
y
), on both the second term on the left side and on the right side of the above equation gives
(30)
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
+
∫
Ω
Φ
(
|
∇
u
|
)
d
x
‍
+
λ
2
∫
Ω
(
v
-
f
)
2
d
x
‍
≥
∫
Ω
Φ
(
|
∇
u
|
)
‍
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
.
Integrating (30) with respect to t we get
(31)
∫
0
t
‍
∫
Ω
∂
t
u
(
v
-
u
)
d
x
d
t
‍
+
∫
0
t
I
(
v
)
d
t
‍
≥
∫
0
t
‍
I
(
u
)
d
t
.
Now, since v
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
), we may choose v
=
u
+
ξ
ϕ, where ϕ
∈
C
0
∞
(
Ω
). We observe that the left-hand side of (31) has a minimum at ξ
=
0. This shows that u is a solution of (12) in the distributional sense. The facts raised above motivate the following definition for a weak solution to (12)–(14).Definition 8.
A functionu
∈
L
2
(
0
,
T
;
BV
(
Ω
)
∩
L
2
) is called a weak solution of (12)–(14) if ∂
t
u
∈
L
2
(
Ω
T
), u
(
x
,
0
)
=
f
(
x
), on Ω and u satisfies (31) for every v
∈
L
2
(
(
0
,
T
)
;
BV
(
Ω
)
∩
L
2
) and t
∈
[
0
,
T
].


### 4.2. Approximating Energy Functional

LetΩ be an open bounded subset of R
n and let f
∈
BV
(
Ω
)
∩
L
∞. Then let us consider the following approximating energy functional for 1
<
p
≤
2:
(32)
I
p
(
u
)
=
∫
Ω
‍
Φ
p
(
∇
u
)
d
x
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
,
where
(33)
Φ
p
(
s
)
=
1
p
(
|
s
|
p
-
1
K
ln
⁡
(
1
+
K
|
s
|
p
)
)
,
with the following properties: (1)
Φ
p
′
(
s
)
=
K
|
s
|
2
p
-
2
s
/
(
1
+
K
|
s
|
p
) and Φ
p
′
(
s
)
≥
0
∀
s
∈
R
n;
(2)
Φ
p
′′
(
s
)
=
K
|
s
|
2
p
-
2
[
(
2
p
-
1
)
+
K
(
p
-
1
)
|
s
|
p
]
/
[
1
+
K
|
s
|
p
]
2
>
0. Hence Φ
(
s
) is strictly convex in s.Remark 9.
Sincef
∈
BV
(
Ω
) we can have a sequence f
κ
∈
W
1
,
p
(
Ω
) with
(34)
f
κ
⟶
f
as
κ
⟶
0
in
L
2
(
Ω
)
,
∫
Ω
‍
|
∇
f
κ
|
⟶
∫
Ω
‍
|
∇
f
|
.
And as a consequence of Theorem  2.9 in [20] we have
(35)
∫
Ω
‍
|
∇
f
κ
|
≤
C
.
From property (iii) of Φ
(
s
) it is easy to see that
(36)
Φ
(
∇
f
)
≤
|
∇
f
|
.
This implies that
(37)
Φ
(
∇
f
κ
)
≤
|
∇
f
κ
|
≤
C
.
Hence
(38)
∫
Ω
‍
Φ
(
∇
f
κ
)
⟶
∫
Ω
‍
Φ
(
∇
f
)
,
as
κ
⟶
0
.On the strength of (32) and Remark 9 let us consider the approximate evolution problem
(39)
∂
t
u
κ
,
p
=
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
-
λ
(
u
κ
,
p
-
f
κ
)
,
in
Ω
×
[
0
,
T
]
,
(40)
∂
u
κ
,
p
∂
n
→
=
0
;
on
⁡
∂
Ω
×
[
0
,
T
]
,
(41)
u
κ
,
p
(
x
,
0
)
=
f
κ
,
x
∈
Ω
.
Let the following L
∞ bound also hold for the solution of (39)–(41). Then, the following lemma indicates the boundedness of the solution u
κ
,
p of the above approximate evolution problem.Lemma 10.
Supposef
κ
∈
L
∞
(
Ω
)
∩
B
V
(
Ω
) and that u
κ
,
p is a solution to the problem (39)–(41); then
(42)
inf
⁡
f
κ
≤
u
κ
,
p
≤
sup
⁡
f
κ
.Proof.
By a method similar to that of Zhou in [33], let M
=
sup
⁡
f
κ and define (
u
κ
,
p
-
M
)
+ such that
(43)
(
u
κ
,
p
-
M
)
+
:
=
{
(
u
κ
,
p
-
M
)
,
if
u
κ
,
p
-
M
≥
0
,
0
,
otherwise
;
then multiplying (39) and integrating over Ω yields
(44)
∫
Ω
‍
∂
t
u
κ
,
p
(
u
κ
,
p
-
M
)
+
d
x
+
∫
Ω
‍
Φ
p
′
(
∇
u
κ
,
p
)
·
∇
u
κ
,
p
d
x
+
λ
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
(
u
κ
,
p
-
M
)
+
d
x
=
0
.
Observe that the last two integrals in the above equation are nonnegative, based on the definition of (
u
κ
,
p
-
M
)
+ and the fact that Φ
p
′
(
s
)
·
s
≥
0. Therefore, we have
(45)
1
2
∫
Ω
‍
d
d
t
(
u
κ
,
p
-
M
)
+
2
d
x
≤
0
,
which indicates that J
(
t
)
=
(
1
/
2
)
∫
Ω
‍
(
u
κ
,
p
-
M
)
+
2
d
x is a decreasing function in t. Since J
(
0
)
=
0 we have that
(46)
J
(
t
)
=
0
,
∀
t
∈
[
0
,
T
]
.
Hence u
κ
,
p
≤
M
=
sup
⁡
f
κ. Conversely, multiplying (39) by (
u
κ
,
p
-
m
)
-, where m
=
inf
⁡
f
κ, and employing a similar argument, we obtain that u
κ
,
p
≥
-
m, ∀
t
∈
[
0
,
T
]. Therefore, inf
⁡
f
κ
≤
u
κ
,
p
≤
sup
⁡
f
κ, ∀
t
∈
[
0
,
T
].Another estimate (bound) is obtained via the following lemma.Lemma 11.
The approximate evolution problem (39)–(41) has a unique solution u
κ
,
p
∈
L
∞
(
0
,
T
;
BV
(
Ω
)
) with ∂
t
u
κ
,
p
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
) such that
(47)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
I
p
(
v
)
d
t
≥
∫
0
t
‍
I
p
(
u
κ
,
p
)
d
t
,
for any t
∈
[
0
,
T
] and v
∈
L
2
(
0
,
T
;
W
1
,
p
(
Ω
)
) with ∂
v
/
∂
n
→
=
0. Moreover
(48)
∫
0
∞
‍
∫
Ω
‍
|
∂
t
u
κ
,
p
|
2
d
x
d
t
+
sup
⁡
t
∈
[
0
,
T
]
{
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
+
λ
2
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
2
d
x
}
≤
I
p
(
∇
f
κ
)
.Proof.
SinceΦ
(
∇
u
) is a lower semicontinuous and strictly convex function, then, by definitions provided in [34, 38, 39], -
div
⁡
(
Φ
′
(
∇
u
)
), which is a subdifferential of ∫
Ω
‍
Φ
(
∇
u
)
d
x, is maximal monotone. The approximate problem (39)–(41) has a unique solution such that u
κ
,
p
∈
L
∞
(
0
,
T
;
W
1
,
p
(
Ω
)
). Now, to show that u
κ
,
p is the weak solution to the approximating problem (39)–(41), as governed by (47), we multiply (39) by v
-
u
κ
,
p for v
∈
L
2
(
0
,
T
;
W
1
,
p
(
Ω
)
) with ∂
v
/
∂
n
→
=
0 and integrate over Ω and t
∈
[
0
,
T
]. Thus
(49)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
=
∫
0
t
‍
∫
Ω
‍
(
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
(
v
-
u
κ
,
p
)
)
d
x
d
t
-
λ
∫
0
t
‍
∫
Ω
‍
(
(
u
κ
,
p
-
f
κ
)
(
v
-
u
κ
,
p
)
)
d
x
d
t
,
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
p
′
(
∇
u
κ
,
p
)
(
∇
v
-
∇
u
κ
,
p
)
d
x
d
t
=
-
λ
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
(
v
-
u
κ
,
p
)
d
x
d
t
.
Applying convexity condition on both sides of the above equation we obtain
(50)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
p
(
∇
v
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
t
‍
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
d
t
+
λ
2
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
k
)
2
d
x
d
t
,
which implies (47).
Then, multiplying (39) by u
˙
κ
,
p and integrating over Ω and t
∈
[
0
,
T
] we obtain
(51)
∫
0
t
‍
∫
Ω
‍
(
∂
t
u
κ
,
p
)
2
d
x
d
t
=
∫
0
t
‍
∫
Ω
‍
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
u
˙
κ
,
p
d
x
d
t
-
λ
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
u
˙
κ
,
p
d
x
d
t
,
which yields
(52)
∫
0
t
‍
∫
Ω
(
∂
t
u
κ
,
p
)
2
d
x
d
t
‍
+
∫
0
t
‍
∂
t
(
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
)
d
t
+
λ
2
∫
0
t
‍
∂
t
(
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
d
x
)
2
d
t
=
0
.
From the above equation together with (41) we obtain
(53)
∫
0
t
‍
∫
Ω
‍
(
∂
t
u
κ
,
p
)
2
d
x
d
t
+
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
|
0
t
d
x
+
λ
2
∫
Ω
(
u
κ
,
p
-
f
κ
)
2
|
0
t
d
x
=
0
‍
,
which by (40) and (41) produces (48).


### 4.3. Existence and Uniqueness of the Solution of the Evolution Problem

Next, using thea priori estimates obtained above, through the approximate problem, we proceed to prove the existence and uniqueness of the solution of the evolution problem (12)–(14).Theorem 12.
Letf
∈
B
V
(
Ω
)
∩
L
∞
(
Ω
). Then there exists a unique weak solution u
∈
L
∞
(
0
,
∞
;
B
V
(
Ω
)
∩
L
∞
(
Ω
)
), ∂
t
u
∈
L
2
(
Q
∞
), and u
(
x
,
0
)
=
f such that
(54)
∫
0
∞
‍
∫
Ω
‍
|
∂
t
u
|
2
d
x
d
t
+
sup
⁡
t
∈
[
0
,
∞
)
{
I
(
u
)
}
≤
I
(
f
)
,
where t
∈
[
0
,
∞
), Q
∞
:
=
Ω
×
[
0
,
∞
).Proof.
From Lemmas10 and 11 we see that
(55)
∥
∂
t
u
κ
,
p
∥
L
2
(
Q
∞
)
+
∥
u
κ
,
p
∥
L
∞
(
Q
∞
)
+
∥
u
κ
,
p
∥
L
∞
(
0
,
∞
;
W
1
,
p
(
Ω
)
∩
L
∞
(
Ω
)
)
≤
M
,
where M is some constant that may vary from line to line. Let {
u
κ
,
p
} by Lemma 10 be a bounded sequence of solutions of the problem (39)–(41). Then there exists a subsequence denoted by {
u
κ
,
p
i
} of {
u
κ
,
p
} and a function u
κ
∈
L
∞
(
Ω
), with u
˙
κ
∈
L
2
(
Ω
;
[
0
,
T
]
), such that, as p
i
→
1,
(56)
u
κ
,
p
i
⟶
u
κ
strongly
in
L
1
(
Ω
)
for
each
t
∈
[
0
,
∞
)
,
(57)
∂
t
u
κ
,
p
i
⇀
∂
t
u
κ
weakly
in
L
2
(
Q
∞
)
,
(58)
u
κ
,
p
i
⇀
*
u
κ
wea
k
*
in
L
∞
(
Q
∞
)
,
(59)
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
uniformly
in
t
,
(60)
lim
⁡
t
→
0
+
∫
Ω
‍
|
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
|
2
d
x
=
0
.
Indeed, from (55), there is a sequence {
u
κ
,
p
i
} and a function u
κ
∈
L
∞
(
Q
∞
) with u
˙
κ
∈
L
2
(
Q
∞
) such that (57) and (58) hold. Observe also that for any ϕ
∈
L
2
(
Ω
), as i
→
∞,
(61)
∫
Ω
‍
(
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
)
ϕ
(
x
)
d
x
=
∫
0
t
‍
∂
s
(
∫
Ω
‍
u
κ
,
p
i
(
x
,
s
)
ϕ
(
x
)
d
x
)
d
s
⟶
∫
0
t
∂
s
(
∫
Ω
‍
u
κ
(
x
,
s
)
ϕ
(
x
)
d
x
)
d
s
‍
=
∫
Ω
(
u
κ
(
x
,
t
)
-
f
κ
(
x
)
)
ϕ
(
x
)
d
x
‍
,
which indicates that for each t
(62)
u
κ
,
p
i
⇀
u
κ
in
L
2
(
Ω
)
.
By Lemmas 10 and 11, for each t
∈
[
0
,
∞
), {
u
κ
,
p
i
(
x
,
t
)
} is a bounded sequence in W
1,1
(
Ω
). Combining that fact with (62) we obtain that for each t as p
i
→
1
(63)
u
κ
,
p
i
⟶
u
κ
in
L
1
(
Ω
)
.
Moreover, (60) follows from the fact that
(64)
∥
u
κ
,
p
i
(
·
,
t
)
-
u
κ
,
p
i
(
·
,
t
′
)
∥
L
2
(
Ω
)
2
≤
|
t
-
t
′
|
∫
0
t
‍
∫
Ω
(
∂
t
u
κ
,
p
i
)
2
d
x
d
t
‍
.
From (64) t
↦
u
κ
,
p
i
(
·
,
t
)
∈
L
2
(
Ω
) is equicontinuous, and from (55) and (58) we have that
(65)
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
.
It then follows by standard argument that we can have u
κ
,
p
i
→
u
κ in L
2
(
Ω
) uniformly in t, giving us (59). From Lemma 10 and (58) we obtain that u
κ
∈
L
∞
(
0
,
∞
,
BV
(
Ω
)
∩
L
∞
(
Ω
)
) with u
˙
κ
∈
L
2
(
Q
∞
).
Next we show that, for allv
∈
L
2
(
0
,
T
,
W
1
,
p
(
Ω
)
∩
L
2
(
Ω
)
), v
>
0 and ∂
v
/
∂
n
→
=
0 and for each t
∈
[
0
,
∞
)
(66)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
(
∇
v
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
t
‍
∫
Ω
‍
Φ
(
∇
u
κ
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
(
u
κ
-
f
k
)
2
d
x
d
t
‍
.
To show this end, from Lemma11, we obtain
(67)
∫
Ω
‍
∂
t
u
κ
,
p
i
(
v
-
u
κ
,
p
i
)
d
x
+
∫
Ω
‍
Φ
p
i
(
∇
v
)
d
x
+
λ
2
∫
Ω
(
v
-
f
κ
)
2
d
x
‍
≥
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
+
λ
2
∫
Ω
(
u
κ
,
p
i
-
f
κ
)
2
d
x
‍
.
From (56) and (58) we deduce that there exists a subsequence {
u
κ
,
p
i
} such that
(68)
u
κ
,
p
i
⟶
u
κ
strongly
in
L
2
(
Ω
;
[
0
,
T
]
)
.
Using (57), (59), and (68) and we let p
i
→
1; in (67) we obtain
(69)
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
+
∫
Ω
‍
Φ
p
i
(
∇
v
)
d
x
+
λ
2
∫
Ω
(
v
-
f
)
2
d
x
‍
≥
liminf
⁡
i
→
∞
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
+
λ
2
∫
Ω
‍
(
u
κ
-
f
κ
)
2
d
x
.
But we define from the lower semicontinuity theorem that
(70)
∫
Ω
Φ
(
∇
u
κ
)
d
x
‍
≤
liminf
⁡
i
→
∞
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
.
Using (70) in (69) and integrating over s
∈
[
0
,
T
] we obtain
(71)
∫
0
s
‍
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
d
t
+
∫
0
s
‍
∫
Ω
‍
Φ
(
∇
v
)
d
x
d
t
+
λ
2
∫
0
s
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
s
‍
∫
Ω
‍
Φ
κ
(
∇
u
κ
)
d
x
d
t
+
λ
2
∫
0
s
‍
∫
Ω
(
u
κ
-
f
κ
)
2
d
x
d
t
‍
,
which confirms (66).
Now, to complete the proof of the existence of solution to (12)–(14), it remains to pass to the limit as κ
→
0 in (71). Replacing p by p
i in (48), letting i
→
∞
(
p
i
→
1
), and using (57)–(59) and (70) we obtain
(72)
∫
0
∞
‍
∫
Ω
|
∂
t
u
κ
|
2
d
x
d
t
‍
+
sup
⁡
t
∈
[
0
,
∞
)
{
∫
Ω
Φ
(
∇
u
κ
)
d
x
‍
+
λ
2
∫
Ω
(
u
κ
-
f
κ
)
2
d
x
‍
}
≤
∫
Ω
Φ
(
∇
f
κ
)
d
x
‍
.
From Lemma 10 and also from the above inequality we see that u
κ is uniformly bounded in L
∞
(
0
,
∞
,
BV
(
Ω
)
∩
L
∞
(
Ω
)
) and ∂
t
u
κ is uniformly bounded in L
2
(
Q
∞
). This means we can extract a subsequence {
u
κ
i
} of {
u
κ
} such that as i
→
∞
(
κ
i
→
0
) we have
(73)
∂
t
u
κ
i
⇀
∂
t
u
weakly
in
L
2
(
Q
∞
)
,
u
κ
i
⇀
*
u
wea
k
*
in
L
∞
(
Ω
∞
)
,
u
κ
i
⟶
u
strongly
in
L
1
(
Ω
,
[
0
,
T
]
)
,
for
each
t
∈
[
0
,
∞
)
,
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
uniformly
in
t
,
lim
⁡
t
→
0
+
∫
Ω
‍
|
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
|
2
d
x
=
0
.
Now, replacing κ with κ
i in (71), letting i
→
∞
(
κ
i
→
0
), and applying the lower semicontinuity in (70) we obtain
(74)
∫
0
s
‍
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
d
t
+
∫
Ω
‍
I
(
v
)
d
t
≥
∫
Ω
‍
I
(
u
)
d
t
,
for all v
∈
L
2
(
0
,
T
,
W
1
,
p
(
Ω
)
∩
L
2
(
Ω
)
), v
>
0 and ∂
v
/
∂
n
→
=
0 and for each t
∈
[
0
,
∞
). Hence u is a weak solution to (12)–(14). Replacing κ by κ
i in (72), letting i
→
∞
(
κ
i
→
0
), and using (57)–(59) and (70) we obtain (54).
Uniqueness of the Weak Solution. With reference to the definition of solution inequality as given in (31), let u
1 and u
2 be two weak solutions to the problem (12) to (14) such that u
1
(
x
,
0
)
=
u
2
(
x
,
0
)
=
f. Then we have, for two solutions, two inequalities:
(75)
∫
0
s
‍
∫
Ω
‍
∂
t
u
1
(
u
2
-
u
1
)
d
x
d
t
+
∫
0
s
‍
I
(
u
2
)
d
t
≥
∫
0
s
‍
I
(
u
1
)
d
t
,
∫
0
s
‍
∫
Ω
‍
∂
t
u
2
(
u
1
-
u
2
)
d
x
d
t
+
∫
0
s
‍
I
(
u
1
)
d
t
≥
∫
0
s
‍
I
(
u
2
)
d
t
.
Now, adding the two inequalities (75) we obtain a more compact inequality given by
(76)
∫
0
s
‍
∫
Ω
‍
(
∂
t
u
2
-
∂
t
u
1
)
(
u
1
-
u
2
)
d
x
d
t
≥
0
.
This implies that
(77)
∫
0
s
‍
d
d
t
∫
Ω
‍
(
u
1
-
u
2
)
2
d
x
d
t
≤
0
,
which implies ∥
u
1
-
u
2
∥
=
0 for a.e. (
x
,
t
)
∈
Q
∞. This confirms the uniqueness of the weak solution.


### 4.4. Large Time Behavior

Finally, we will assess the asymptotic limit of the weak solutionu
(
·
,
t
) as t
→
∞. The essence of this part is to demonstrate that over time the solution of the evolution problem (12)–(14) ultimately converges to the unique minimizer of the functional (6).Theorem 13.
Ast
→
∞ the weak solution u of the evolution equation (12)–(14) converges weakly in L
2
(
Ω
) to a minimizer u
- of the functional I
(
u
) in (6).Proof.
Letv
∈
BV
(
Ω
)
∩
L
∞
(
Ω
) into (31) to obtain
(78)
-
1
2
∫
Ω
‍
(
v
(
x
)
-
u
)
2
|
0
s
d
x
+
∫
0
s
‍
I
(
v
(
x
)
)
d
t
≥
∫
0
s
‍
I
(
u
)
d
t
.
The equation above simplifies to
(79)
∫
Ω
‍
(
u
(
x
,
s
)
-
f
)
v
(
x
)
d
x
-
1
2
∫
Ω
‍
(
u
2
(
x
,
s
)
-
f
2
)
d
x
+
s
I
(
v
(
x
)
)
≥
∫
0
s
‍
I
(
u
)
d
t
.
By taking u
-
(
x
,
s
)
=
(
1
/
s
)
∫
0
s
‍
u
(
x
,
t
)
d
t, and have for each s, u
-
(
x
,
s
)
∈
BV
(
Ω
)
∩
L
∞
(
Ω
). Given that u
- is uniformly bounded in BV
(
Ω
), we conclude that there exists sequence u
-
(
x
,
s
n
) and its subsequence is still denoted by {
u
-
(
x
,
s
n
)
} such that u
-
(
x
,
s
n
)
→
u
-
(
x
,
s
) in L
1
(
Ω
) and u
-
(
x
,
s
n
)
⇀
u
-
(
x
,
s
) weak* in BV
(
Ω
), as s
n
→
∞, respectively. Hence dividing inequality (79) by s and taking the limit as s
n
→
∞ we obtain
(80)
I
(
v
)
≥
I
(
u
-
)
.
This demonstrates that indeed u
- is a weak solution to (12)–(14) and also the unique minimizer of problem (6).


## 4.1. Definition of Weak Solution

Letv
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
), f
∈
BV
(
Ω
), v
>
0, and ∂
v
/
∂
n
→
=
0, where Ω
T
=
Ω
×
[
0
,
T
] and u is a solution to (12)–(14). Multiplying (12) by (
v
-
u
) and integrating over Ω, we have
(29)
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
+
∫
Ω
‍
K
∇
u
1
+
K
|
∇
u
|
(
∇
v
-
∇
u
)
d
x
=
-
λ
∫
Ω
‍
(
u
-
f
)
(
v
-
u
)
d
x
.
Then applying the convexity condition of Φ, that is, Φ
(
x
)
-
Φ
(
y
)
≥
Φ
′
(
y
)
(
x
-
y
), on both the second term on the left side and on the right side of the above equation gives
(30)
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
+
∫
Ω
Φ
(
|
∇
u
|
)
d
x
‍
+
λ
2
∫
Ω
(
v
-
f
)
2
d
x
‍
≥
∫
Ω
Φ
(
|
∇
u
|
)
‍
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
.
Integrating (30) with respect to t we get
(31)
∫
0
t
‍
∫
Ω
∂
t
u
(
v
-
u
)
d
x
d
t
‍
+
∫
0
t
I
(
v
)
d
t
‍
≥
∫
0
t
‍
I
(
u
)
d
t
.
Now, since v
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
), we may choose v
=
u
+
ξ
ϕ, where ϕ
∈
C
0
∞
(
Ω
). We observe that the left-hand side of (31) has a minimum at ξ
=
0. This shows that u is a solution of (12) in the distributional sense. The facts raised above motivate the following definition for a weak solution to (12)–(14).Definition 8.
A functionu
∈
L
2
(
0
,
T
;
BV
(
Ω
)
∩
L
2
) is called a weak solution of (12)–(14) if ∂
t
u
∈
L
2
(
Ω
T
), u
(
x
,
0
)
=
f
(
x
), on Ω and u satisfies (31) for every v
∈
L
2
(
(
0
,
T
)
;
BV
(
Ω
)
∩
L
2
) and t
∈
[
0
,
T
].


## 4.2. Approximating Energy Functional

LetΩ be an open bounded subset of R
n and let f
∈
BV
(
Ω
)
∩
L
∞. Then let us consider the following approximating energy functional for 1
<
p
≤
2:
(32)
I
p
(
u
)
=
∫
Ω
‍
Φ
p
(
∇
u
)
d
x
+
λ
2
∫
Ω
‍
(
u
-
f
)
2
d
x
,
where
(33)
Φ
p
(
s
)
=
1
p
(
|
s
|
p
-
1
K
ln
⁡
(
1
+
K
|
s
|
p
)
)
,
with the following properties: (1)
Φ
p
′
(
s
)
=
K
|
s
|
2
p
-
2
s
/
(
1
+
K
|
s
|
p
) and Φ
p
′
(
s
)
≥
0
∀
s
∈
R
n;
(2)
Φ
p
′′
(
s
)
=
K
|
s
|
2
p
-
2
[
(
2
p
-
1
)
+
K
(
p
-
1
)
|
s
|
p
]
/
[
1
+
K
|
s
|
p
]
2
>
0. Hence Φ
(
s
) is strictly convex in s.Remark 9.
Sincef
∈
BV
(
Ω
) we can have a sequence f
κ
∈
W
1
,
p
(
Ω
) with
(34)
f
κ
⟶
f
as
κ
⟶
0
in
L
2
(
Ω
)
,
∫
Ω
‍
|
∇
f
κ
|
⟶
∫
Ω
‍
|
∇
f
|
.
And as a consequence of Theorem  2.9 in [20] we have
(35)
∫
Ω
‍
|
∇
f
κ
|
≤
C
.
From property (iii) of Φ
(
s
) it is easy to see that
(36)
Φ
(
∇
f
)
≤
|
∇
f
|
.
This implies that
(37)
Φ
(
∇
f
κ
)
≤
|
∇
f
κ
|
≤
C
.
Hence
(38)
∫
Ω
‍
Φ
(
∇
f
κ
)
⟶
∫
Ω
‍
Φ
(
∇
f
)
,
as
κ
⟶
0
.On the strength of (32) and Remark 9 let us consider the approximate evolution problem
(39)
∂
t
u
κ
,
p
=
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
-
λ
(
u
κ
,
p
-
f
κ
)
,
in
Ω
×
[
0
,
T
]
,
(40)
∂
u
κ
,
p
∂
n
→
=
0
;
on
⁡
∂
Ω
×
[
0
,
T
]
,
(41)
u
κ
,
p
(
x
,
0
)
=
f
κ
,
x
∈
Ω
.
Let the following L
∞ bound also hold for the solution of (39)–(41). Then, the following lemma indicates the boundedness of the solution u
κ
,
p of the above approximate evolution problem.Lemma 10.
Supposef
κ
∈
L
∞
(
Ω
)
∩
B
V
(
Ω
) and that u
κ
,
p is a solution to the problem (39)–(41); then
(42)
inf
⁡
f
κ
≤
u
κ
,
p
≤
sup
⁡
f
κ
.Proof.
By a method similar to that of Zhou in [33], let M
=
sup
⁡
f
κ and define (
u
κ
,
p
-
M
)
+ such that
(43)
(
u
κ
,
p
-
M
)
+
:
=
{
(
u
κ
,
p
-
M
)
,
if
u
κ
,
p
-
M
≥
0
,
0
,
otherwise
;
then multiplying (39) and integrating over Ω yields
(44)
∫
Ω
‍
∂
t
u
κ
,
p
(
u
κ
,
p
-
M
)
+
d
x
+
∫
Ω
‍
Φ
p
′
(
∇
u
κ
,
p
)
·
∇
u
κ
,
p
d
x
+
λ
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
(
u
κ
,
p
-
M
)
+
d
x
=
0
.
Observe that the last two integrals in the above equation are nonnegative, based on the definition of (
u
κ
,
p
-
M
)
+ and the fact that Φ
p
′
(
s
)
·
s
≥
0. Therefore, we have
(45)
1
2
∫
Ω
‍
d
d
t
(
u
κ
,
p
-
M
)
+
2
d
x
≤
0
,
which indicates that J
(
t
)
=
(
1
/
2
)
∫
Ω
‍
(
u
κ
,
p
-
M
)
+
2
d
x is a decreasing function in t. Since J
(
0
)
=
0 we have that
(46)
J
(
t
)
=
0
,
∀
t
∈
[
0
,
T
]
.
Hence u
κ
,
p
≤
M
=
sup
⁡
f
κ. Conversely, multiplying (39) by (
u
κ
,
p
-
m
)
-, where m
=
inf
⁡
f
κ, and employing a similar argument, we obtain that u
κ
,
p
≥
-
m, ∀
t
∈
[
0
,
T
]. Therefore, inf
⁡
f
κ
≤
u
κ
,
p
≤
sup
⁡
f
κ, ∀
t
∈
[
0
,
T
].Another estimate (bound) is obtained via the following lemma.Lemma 11.
The approximate evolution problem (39)–(41) has a unique solution u
κ
,
p
∈
L
∞
(
0
,
T
;
BV
(
Ω
)
) with ∂
t
u
κ
,
p
∈
L
2
(
0
,
T
;
L
2
(
Ω
)
) such that
(47)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
I
p
(
v
)
d
t
≥
∫
0
t
‍
I
p
(
u
κ
,
p
)
d
t
,
for any t
∈
[
0
,
T
] and v
∈
L
2
(
0
,
T
;
W
1
,
p
(
Ω
)
) with ∂
v
/
∂
n
→
=
0. Moreover
(48)
∫
0
∞
‍
∫
Ω
‍
|
∂
t
u
κ
,
p
|
2
d
x
d
t
+
sup
⁡
t
∈
[
0
,
T
]
{
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
+
λ
2
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
2
d
x
}
≤
I
p
(
∇
f
κ
)
.Proof.
SinceΦ
(
∇
u
) is a lower semicontinuous and strictly convex function, then, by definitions provided in [34, 38, 39], -
div
⁡
(
Φ
′
(
∇
u
)
), which is a subdifferential of ∫
Ω
‍
Φ
(
∇
u
)
d
x, is maximal monotone. The approximate problem (39)–(41) has a unique solution such that u
κ
,
p
∈
L
∞
(
0
,
T
;
W
1
,
p
(
Ω
)
). Now, to show that u
κ
,
p is the weak solution to the approximating problem (39)–(41), as governed by (47), we multiply (39) by v
-
u
κ
,
p for v
∈
L
2
(
0
,
T
;
W
1
,
p
(
Ω
)
) with ∂
v
/
∂
n
→
=
0 and integrate over Ω and t
∈
[
0
,
T
]. Thus
(49)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
=
∫
0
t
‍
∫
Ω
‍
(
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
(
v
-
u
κ
,
p
)
)
d
x
d
t
-
λ
∫
0
t
‍
∫
Ω
‍
(
(
u
κ
,
p
-
f
κ
)
(
v
-
u
κ
,
p
)
)
d
x
d
t
,
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
p
′
(
∇
u
κ
,
p
)
(
∇
v
-
∇
u
κ
,
p
)
d
x
d
t
=
-
λ
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
(
v
-
u
κ
,
p
)
d
x
d
t
.
Applying convexity condition on both sides of the above equation we obtain
(50)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
,
p
(
v
-
u
κ
,
p
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
p
(
∇
v
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
t
‍
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
d
t
+
λ
2
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
k
)
2
d
x
d
t
,
which implies (47).
Then, multiplying (39) by u
˙
κ
,
p and integrating over Ω and t
∈
[
0
,
T
] we obtain
(51)
∫
0
t
‍
∫
Ω
‍
(
∂
t
u
κ
,
p
)
2
d
x
d
t
=
∫
0
t
‍
∫
Ω
‍
div
⁡
(
Φ
p
′
(
∇
u
κ
,
p
)
)
u
˙
κ
,
p
d
x
d
t
-
λ
∫
0
t
‍
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
u
˙
κ
,
p
d
x
d
t
,
which yields
(52)
∫
0
t
‍
∫
Ω
(
∂
t
u
κ
,
p
)
2
d
x
d
t
‍
+
∫
0
t
‍
∂
t
(
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
d
x
)
d
t
+
λ
2
∫
0
t
‍
∂
t
(
∫
Ω
‍
(
u
κ
,
p
-
f
κ
)
d
x
)
2
d
t
=
0
.
From the above equation together with (41) we obtain
(53)
∫
0
t
‍
∫
Ω
‍
(
∂
t
u
κ
,
p
)
2
d
x
d
t
+
∫
Ω
‍
Φ
p
(
∇
u
κ
,
p
)
|
0
t
d
x
+
λ
2
∫
Ω
(
u
κ
,
p
-
f
κ
)
2
|
0
t
d
x
=
0
‍
,
which by (40) and (41) produces (48).


## 4.3. Existence and Uniqueness of the Solution of the Evolution Problem

Next, using thea priori estimates obtained above, through the approximate problem, we proceed to prove the existence and uniqueness of the solution of the evolution problem (12)–(14).Theorem 12.
Letf
∈
B
V
(
Ω
)
∩
L
∞
(
Ω
). Then there exists a unique weak solution u
∈
L
∞
(
0
,
∞
;
B
V
(
Ω
)
∩
L
∞
(
Ω
)
), ∂
t
u
∈
L
2
(
Q
∞
), and u
(
x
,
0
)
=
f such that
(54)
∫
0
∞
‍
∫
Ω
‍
|
∂
t
u
|
2
d
x
d
t
+
sup
⁡
t
∈
[
0
,
∞
)
{
I
(
u
)
}
≤
I
(
f
)
,
where t
∈
[
0
,
∞
), Q
∞
:
=
Ω
×
[
0
,
∞
).Proof.
From Lemmas10 and 11 we see that
(55)
∥
∂
t
u
κ
,
p
∥
L
2
(
Q
∞
)
+
∥
u
κ
,
p
∥
L
∞
(
Q
∞
)
+
∥
u
κ
,
p
∥
L
∞
(
0
,
∞
;
W
1
,
p
(
Ω
)
∩
L
∞
(
Ω
)
)
≤
M
,
where M is some constant that may vary from line to line. Let {
u
κ
,
p
} by Lemma 10 be a bounded sequence of solutions of the problem (39)–(41). Then there exists a subsequence denoted by {
u
κ
,
p
i
} of {
u
κ
,
p
} and a function u
κ
∈
L
∞
(
Ω
), with u
˙
κ
∈
L
2
(
Ω
;
[
0
,
T
]
), such that, as p
i
→
1,
(56)
u
κ
,
p
i
⟶
u
κ
strongly
in
L
1
(
Ω
)
for
each
t
∈
[
0
,
∞
)
,
(57)
∂
t
u
κ
,
p
i
⇀
∂
t
u
κ
weakly
in
L
2
(
Q
∞
)
,
(58)
u
κ
,
p
i
⇀
*
u
κ
wea
k
*
in
L
∞
(
Q
∞
)
,
(59)
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
uniformly
in
t
,
(60)
lim
⁡
t
→
0
+
∫
Ω
‍
|
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
|
2
d
x
=
0
.
Indeed, from (55), there is a sequence {
u
κ
,
p
i
} and a function u
κ
∈
L
∞
(
Q
∞
) with u
˙
κ
∈
L
2
(
Q
∞
) such that (57) and (58) hold. Observe also that for any ϕ
∈
L
2
(
Ω
), as i
→
∞,
(61)
∫
Ω
‍
(
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
)
ϕ
(
x
)
d
x
=
∫
0
t
‍
∂
s
(
∫
Ω
‍
u
κ
,
p
i
(
x
,
s
)
ϕ
(
x
)
d
x
)
d
s
⟶
∫
0
t
∂
s
(
∫
Ω
‍
u
κ
(
x
,
s
)
ϕ
(
x
)
d
x
)
d
s
‍
=
∫
Ω
(
u
κ
(
x
,
t
)
-
f
κ
(
x
)
)
ϕ
(
x
)
d
x
‍
,
which indicates that for each t
(62)
u
κ
,
p
i
⇀
u
κ
in
L
2
(
Ω
)
.
By Lemmas 10 and 11, for each t
∈
[
0
,
∞
), {
u
κ
,
p
i
(
x
,
t
)
} is a bounded sequence in W
1,1
(
Ω
). Combining that fact with (62) we obtain that for each t as p
i
→
1
(63)
u
κ
,
p
i
⟶
u
κ
in
L
1
(
Ω
)
.
Moreover, (60) follows from the fact that
(64)
∥
u
κ
,
p
i
(
·
,
t
)
-
u
κ
,
p
i
(
·
,
t
′
)
∥
L
2
(
Ω
)
2
≤
|
t
-
t
′
|
∫
0
t
‍
∫
Ω
(
∂
t
u
κ
,
p
i
)
2
d
x
d
t
‍
.
From (64) t
↦
u
κ
,
p
i
(
·
,
t
)
∈
L
2
(
Ω
) is equicontinuous, and from (55) and (58) we have that
(65)
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
.
It then follows by standard argument that we can have u
κ
,
p
i
→
u
κ in L
2
(
Ω
) uniformly in t, giving us (59). From Lemma 10 and (58) we obtain that u
κ
∈
L
∞
(
0
,
∞
,
BV
(
Ω
)
∩
L
∞
(
Ω
)
) with u
˙
κ
∈
L
2
(
Q
∞
).
Next we show that, for allv
∈
L
2
(
0
,
T
,
W
1
,
p
(
Ω
)
∩
L
2
(
Ω
)
), v
>
0 and ∂
v
/
∂
n
→
=
0 and for each t
∈
[
0
,
∞
)
(66)
∫
0
t
‍
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
d
t
+
∫
0
t
‍
∫
Ω
‍
Φ
(
∇
v
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
t
‍
∫
Ω
‍
Φ
(
∇
u
κ
)
d
x
d
t
-
λ
2
∫
0
t
‍
∫
Ω
(
u
κ
-
f
k
)
2
d
x
d
t
‍
.
To show this end, from Lemma11, we obtain
(67)
∫
Ω
‍
∂
t
u
κ
,
p
i
(
v
-
u
κ
,
p
i
)
d
x
+
∫
Ω
‍
Φ
p
i
(
∇
v
)
d
x
+
λ
2
∫
Ω
(
v
-
f
κ
)
2
d
x
‍
≥
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
+
λ
2
∫
Ω
(
u
κ
,
p
i
-
f
κ
)
2
d
x
‍
.
From (56) and (58) we deduce that there exists a subsequence {
u
κ
,
p
i
} such that
(68)
u
κ
,
p
i
⟶
u
κ
strongly
in
L
2
(
Ω
;
[
0
,
T
]
)
.
Using (57), (59), and (68) and we let p
i
→
1; in (67) we obtain
(69)
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
+
∫
Ω
‍
Φ
p
i
(
∇
v
)
d
x
+
λ
2
∫
Ω
(
v
-
f
)
2
d
x
‍
≥
liminf
⁡
i
→
∞
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
+
λ
2
∫
Ω
‍
(
u
κ
-
f
κ
)
2
d
x
.
But we define from the lower semicontinuity theorem that
(70)
∫
Ω
Φ
(
∇
u
κ
)
d
x
‍
≤
liminf
⁡
i
→
∞
∫
Ω
Φ
p
i
(
∇
u
κ
,
p
i
)
d
x
‍
.
Using (70) in (69) and integrating over s
∈
[
0
,
T
] we obtain
(71)
∫
0
s
‍
∫
Ω
‍
∂
t
u
κ
(
v
-
u
κ
)
d
x
d
t
+
∫
0
s
‍
∫
Ω
‍
Φ
(
∇
v
)
d
x
d
t
+
λ
2
∫
0
s
‍
∫
Ω
‍
(
v
-
f
κ
)
2
d
x
d
t
≥
∫
0
s
‍
∫
Ω
‍
Φ
κ
(
∇
u
κ
)
d
x
d
t
+
λ
2
∫
0
s
‍
∫
Ω
(
u
κ
-
f
κ
)
2
d
x
d
t
‍
,
which confirms (66).
Now, to complete the proof of the existence of solution to (12)–(14), it remains to pass to the limit as κ
→
0 in (71). Replacing p by p
i in (48), letting i
→
∞
(
p
i
→
1
), and using (57)–(59) and (70) we obtain
(72)
∫
0
∞
‍
∫
Ω
|
∂
t
u
κ
|
2
d
x
d
t
‍
+
sup
⁡
t
∈
[
0
,
∞
)
{
∫
Ω
Φ
(
∇
u
κ
)
d
x
‍
+
λ
2
∫
Ω
(
u
κ
-
f
κ
)
2
d
x
‍
}
≤
∫
Ω
Φ
(
∇
f
κ
)
d
x
‍
.
From Lemma 10 and also from the above inequality we see that u
κ is uniformly bounded in L
∞
(
0
,
∞
,
BV
(
Ω
)
∩
L
∞
(
Ω
)
) and ∂
t
u
κ is uniformly bounded in L
2
(
Q
∞
). This means we can extract a subsequence {
u
κ
i
} of {
u
κ
} such that as i
→
∞
(
κ
i
→
0
) we have
(73)
∂
t
u
κ
i
⇀
∂
t
u
weakly
in
L
2
(
Q
∞
)
,
u
κ
i
⇀
*
u
wea
k
*
in
L
∞
(
Ω
∞
)
,
u
κ
i
⟶
u
strongly
in
L
1
(
Ω
,
[
0
,
T
]
)
,
for
each
t
∈
[
0
,
∞
)
,
u
κ
,
p
i
⟶
u
κ
in
L
2
(
Ω
)
uniformly
in
t
,
lim
⁡
t
→
0
+
∫
Ω
‍
|
u
κ
,
p
i
(
x
,
t
)
-
f
κ
(
x
)
|
2
d
x
=
0
.
Now, replacing κ with κ
i in (71), letting i
→
∞
(
κ
i
→
0
), and applying the lower semicontinuity in (70) we obtain
(74)
∫
0
s
‍
∫
Ω
‍
∂
t
u
(
v
-
u
)
d
x
d
t
+
∫
Ω
‍
I
(
v
)
d
t
≥
∫
Ω
‍
I
(
u
)
d
t
,
for all v
∈
L
2
(
0
,
T
,
W
1
,
p
(
Ω
)
∩
L
2
(
Ω
)
), v
>
0 and ∂
v
/
∂
n
→
=
0 and for each t
∈
[
0
,
∞
). Hence u is a weak solution to (12)–(14). Replacing κ by κ
i in (72), letting i
→
∞
(
κ
i
→
0
), and using (57)–(59) and (70) we obtain (54).
Uniqueness of the Weak Solution. With reference to the definition of solution inequality as given in (31), let u
1 and u
2 be two weak solutions to the problem (12) to (14) such that u
1
(
x
,
0
)
=
u
2
(
x
,
0
)
=
f. Then we have, for two solutions, two inequalities:
(75)
∫
0
s
‍
∫
Ω
‍
∂
t
u
1
(
u
2
-
u
1
)
d
x
d
t
+
∫
0
s
‍
I
(
u
2
)
d
t
≥
∫
0
s
‍
I
(
u
1
)
d
t
,
∫
0
s
‍
∫
Ω
‍
∂
t
u
2
(
u
1
-
u
2
)
d
x
d
t
+
∫
0
s
‍
I
(
u
1
)
d
t
≥
∫
0
s
‍
I
(
u
2
)
d
t
.
Now, adding the two inequalities (75) we obtain a more compact inequality given by
(76)
∫
0
s
‍
∫
Ω
‍
(
∂
t
u
2
-
∂
t
u
1
)
(
u
1
-
u
2
)
d
x
d
t
≥
0
.
This implies that
(77)
∫
0
s
‍
d
d
t
∫
Ω
‍
(
u
1
-
u
2
)
2
d
x
d
t
≤
0
,
which implies ∥
u
1
-
u
2
∥
=
0 for a.e. (
x
,
t
)
∈
Q
∞. This confirms the uniqueness of the weak solution.


## 4.4. Large Time Behavior

Finally, we will assess the asymptotic limit of the weak solutionu
(
·
,
t
) as t
→
∞. The essence of this part is to demonstrate that over time the solution of the evolution problem (12)–(14) ultimately converges to the unique minimizer of the functional (6).Theorem 13.
Ast
→
∞ the weak solution u of the evolution equation (12)–(14) converges weakly in L
2
(
Ω
) to a minimizer u
- of the functional I
(
u
) in (6).Proof.
Letv
∈
BV
(
Ω
)
∩
L
∞
(
Ω
) into (31) to obtain
(78)
-
1
2
∫
Ω
‍
(
v
(
x
)
-
u
)
2
|
0
s
d
x
+
∫
0
s
‍
I
(
v
(
x
)
)
d
t
≥
∫
0
s
‍
I
(
u
)
d
t
.
The equation above simplifies to
(79)
∫
Ω
‍
(
u
(
x
,
s
)
-
f
)
v
(
x
)
d
x
-
1
2
∫
Ω
‍
(
u
2
(
x
,
s
)
-
f
2
)
d
x
+
s
I
(
v
(
x
)
)
≥
∫
0
s
‍
I
(
u
)
d
t
.
By taking u
-
(
x
,
s
)
=
(
1
/
s
)
∫
0
s
‍
u
(
x
,
t
)
d
t, and have for each s, u
-
(
x
,
s
)
∈
BV
(
Ω
)
∩
L
∞
(
Ω
). Given that u
- is uniformly bounded in BV
(
Ω
), we conclude that there exists sequence u
-
(
x
,
s
n
) and its subsequence is still denoted by {
u
-
(
x
,
s
n
)
} such that u
-
(
x
,
s
n
)
→
u
-
(
x
,
s
) in L
1
(
Ω
) and u
-
(
x
,
s
n
)
⇀
u
-
(
x
,
s
) weak* in BV
(
Ω
), as s
n
→
∞, respectively. Hence dividing inequality (79) by s and taking the limit as s
n
→
∞ we obtain
(80)
I
(
v
)
≥
I
(
u
-
)
.
This demonstrates that indeed u
- is a weak solution to (12)–(14) and also the unique minimizer of problem (6).


## 5. Numerical Experiments

In this section, we present the performance of our method in denoising images involving a Gaussian white noise. We have then compared our results with the ones obtained by the classical methods of PM method [7], TV method [8], and the more recent D-α-PM method [11].


### 5.1. Numerical Scheme

In the following two subsections, two numerical discrete schemes, the PM scheme (PMS) and the additive operator splitting (AOS) scheme, have been proposed.


#### 5.1.1. PM Scheme

Here, we have proposed a numerical scheme similar to the original PM method, whereby (12)–(14) are discretized as follows:
(81)
C
N
i
,
j
n
=
K
1
+
K
|
∇
N
u
i
,
j
|
,
C
S
i
,
j
n
=
K
1
+
K
|
∇
S
u
i
,
j
|
,
C
E
i
,
j
n
=
K
1
+
K
|
∇
E
u
i
,
j
|
,
C
W
i
,
j
n
=
K
1
+
K
|
∇
W
u
i
,
j
|
,
di
v
i
,
j
n
=
[
C
N
i
,
j
n
∇
N
u
i
,
j
+
C
S
i
,
j
n
∇
S
u
i
,
j
+
C
E
i
,
j
n
∇
E
u
i
,
j
+
C
W
i
,
j
n
∇
W
u
i
,
j
]
,
and λ is dynamically determined according to the following discretization scheme:
(82)
λ
n
=
1
σ
2
|
Ω
|
∑
i
,
j
‍
di
v
i
,
j
n
(
u
i
,
j
-
f
i
,
j
)
,
where
|
Ω
|
=
M
N
is
the
size
of
image
.
Hence from (81) and (82) we have
(83)
u
i
,
j
n
+
1
=
u
i
,
j
n
+
τ
di
v
i
,
j
n
-
λ
n
τ
(
u
i
,
j
-
f
i
,
j
)
,
u
i
,
j
0
=
f
i
,
j
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
M
,
j
n
=
u
M
-
1
,
j
n
,
u
i
,
N
n
=
u
i
,
N
-
1
n
,
where
(84)
∇
N
u
i
,
j
=
u
i
-
1
,
j
-
u
i
,
j
,
∇
S
u
i
,
j
=
u
i
+
1
,
j
-
u
i
,
j
,
∇
E
u
i
,
j
=
u
i
,
j
+
1
-
u
i
,
j
,
∇
W
u
i
,
j
=
u
i
,
j
-
1
-
u
i
,
j
,
for i
=
0,1
,
2
,
…
,
N and j
=
0,1
,
2
,
…
,
M.


#### 5.1.2. AOS Scheme

In this part, using a AOS scheme, the problem (12)–(14) has been discretized as follows:
(85)
λ
0
=
0
,
u
n
+
1
=
1
m
∑
l
=
1
m
‍
[
I
-
m
τ
A
l
(
u
k
)
]
-
1
[
u
n
+
λ
τ
(
f
-
u
n
)
]
,
div
⁡
n
=
(
u
n
+
1
-
u
n
)
τ
,
λ
n
=
1
σ
2
M
N
(
u
-
f
)
div
⁡
n
,
u
i
,
j
0
=
f
i
,
j
=
f
(
i
h
,
j
h
)
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
I
,
i
n
=
u
I
-
1
,
i
n
,
u
i
,
J
n
=
u
i
,
J
-
1
n
,
where A
l
(
u
n
)
=
[
a
i
,
j
(
u
n
)
],
(86)
a
i
,
j
(
u
n
)
:
=
{
C
i
n
+
C
j
n
2
h
2
,
[
j
∈
N
(
i
)
]
,
-
∑
N
∈
N
(
i
)
C
i
n
+
C
N
n
2
h
2
‍
,
(
j
=
i
)
,
0
,
(
else
)
,
C
i
n
:
=
K
1
+
K
|
∇
u
i
,
j
n
|
,
where
(87)
|
∇
u
i
,
j
n
|
=
1
2
∑
p
,
q
∈
N
(
i
)
|
u
p
n
-
u
q
n
|
2
h
‍
,
where N
(
i
) is the set of the two neighbors of pixel i (boundary pixels have only one neighbor).It is observed that AOS schemes with large time steps still reveal average grey value invariance, stability based on extremum principle, Lyapunov functionals, and convergence to a constant steady state [10]. The AOS scheme is less than twice the typical effort needed for the PM scheme, a rather low price for gaining absolute stability [10]. It is worth noting that if we use div
⁡
(
∇
u
/
ε
+
|
∇
u
|
2
) to approximate div
⁡
(
∇
u
/
|
∇
u
|
) (TV kernel), in numerical scheme, the AOS scheme may be unstable because of the small number ε. Our approximation, however, effectively avoids instabilities arising from such a scenario, as there is no need for lifting the denominator, since it cannot not be equal to zero. This has been an additional motivation for considering AOS for our numerical experiments. The possibility of occurrence of a zero denominator in the evolution problem is a phenomenon that makes the numerical implementation of TV problematic.


### 5.2. Comparison with Other Methods

The experiments in this work were performed on a Compaq610 computer, having Intel Core 2 Duo CPU T5870 each 2.00 GHz, physical RAM of 4.00 GB, and Professional Windows 8 64-bit Operating System, on MATLAB R2013b. The image restoration performance was measured in terms of the peak-signal-to-noise ratio (PSNR), mean absolute deviation/error (MAE), structural similarity index measure (SSIM), the measure of similarity of edges (PSNRE), and visual effects. The iteration stopping mechanism was based on the maximal PSNR or the lowest MAE. The PSNR and MAE values were obtained using the formula by Durand et al. [40] and are given by
(88)
PSNR
=
10
log
⁡
10
M
N
|
max
⁡
u
0
-
min
⁡
u
0
|
2
∥
u
-
u
0
∥
L
2
2
d
B
,
MAE
=
∥
u
-
u
0
∥
L
1
M
N
.
Also, taking the edge map EM
(
u
) from the evolution equation, we have
(89)
EM
(
u
)
=
K
(
1
+
K
|
∇
u
|
)
.
And corresponding edge similarity measure is given, according to the formulation by Guo et al. [11], by
(90)
PSNR
E
=
PSNR
(
EM
(
u
)
,
EM
(
u
0
)
)
d
B
.
The SSIM measures have been obtained according to the formula by Wang et al. [41] given by
(91)
SSIM
(
u
,
u
0
)
=
L
(
u
,
u
0
)
·
C
(
u
,
u
0
)
·
R
(
u
,
u
0
)
,
where u
0 denotes the noise-free image, u is the denoised image, M
×
N is the dimension of image, and |

max
⁡
u
0
-
min
⁡
u
0
| yields the gray scale range of the original image. In addition, L
(
u
,
u
0
)
=
(
2
μ
u
μ
u
0
+
k
1
)
/
(
μ
u
2
+
μ
u
0
2
+
k
1
) compares the two images’ mean luminances μ
u and μ
u
0. The maximal value of L
(
u
,
u
0
)
=
1, if μ
u
=
μ
u
0 and C
(
u
,
u
0
)
=
(
2
σ
u
σ
u
0
+
k
2
)
/
(
σ
u
2
+
σ
u
0
2
+
k
2
), measures the closeness of contrast of the two images u and u
0. Contrast is determined in terms of standard deviation, σ. Contrast comparison measure C
(
u
,
u
0
)
=
1 maximally if and only if σ
u
=
σ
u
0, that is, when the images have equal contrast.R
(
u
,
u
0
)
=
(
σ
u
u
0
+
k
3
)
/
(
σ
u
σ
u
0
+
k
3
), where σ
u
u
0 is covariance between u and u
0, is a structure comparison measure which determines the correlation between the images u and u
0. It attains maximal value of 1 if structurally the two images coincide, but its value is equal to zero when there is absolutely no structural coincidence. The quantities k
1, k
2, and k
3 are small positive constants that avert the possibility of having zero denominators.The results of our method were compared to PM method [7], TV method [8], and the D-α-PM method [11]. Tables 1 and 2, respectively, give a summary of the results from the experiments, using synthetic image in Figure 1 and Lena image in Figure 3. The parameters considered here are thresholding parameter K, variance δ, the time step parameter τ, and the convolution parameter σ
1 applied in the D-α-PM method. The fidelity parameter λ was dynamically obtained according to (82) or under the AOS scheme in Section 5.1.2, while the rest of the parameters were chosen to guarantee stability and attainment of optimal results.Table 1
Numerical results for synthetic image(
300
×
300
).
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
PSNRE
SSIM
σ
K
τ
PM
30
12
0.25
232
12.73
34.23
2.84
25.26
0.9829
TV
30
n/a
0.2
203
17.86
32.79
3.76
15.76
0.9748
D-α-PM
30
1
0.25
42
2.40
37.30
3.38
24.76
0.9838
PMS
30
n/a
0.20
262
12.01
35.75
2.15
25.60
0.9870
AOS
30
n/a
3.00
17
2.51
33.30
2.92
25.40
0.9839Table 2
Numerical results for Lena image(
300
×
300
).
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
PSNRE
SSIM
σ
K
τ
PM
30
12
0.25
60
2.43
27.10
7.75
21.70
0.7234
TV
30
n/a
0.2
149
13.87
27.46
7.49
22.63
0.7960
D-α-PM
30
4
2
13
0.90
28.17
7.35
24.78
0.7880
PMS
30
n/a
0.2
122
5.80
27.73
7.45
26.60
0.8060
AOS
30
n/a
2.0
7
1.24
28.06
7.27
27.80
0.8509Synthetic image (300
×
300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ
=
30. (c) Our algorithm by PMS; τ
=
0.2 (262 steps). (d) Our algorithm by AOS scheme; τ
=
3 (17 steps). (e) PM method; K
=
5 and τ
=
0.25 (232 steps). (f) TV method; τ
=
0.2 (203 steps). (g) D-α-PM method; σ
1
=
0.5, σ
=
30, τ
=
0.25, and K
=
1 (42 steps).
(a)
Original image
(b)
Noisy image
(c)
Proposed model I (PMS)
(d)
Proposed model I (AOS)
(e)
PM model
(f)
TV model
(g)
D-α-PM modelFor nontexture images as displayed in Figure1 and the results of our method using PMS scheme, as shown in the fourth row of Table 1, demonstrate better performance than those PM and TV as indicated by the higher PSNR, PSNRE, and SSIM values and lower MAE. And, in spite of higher iterative steps, the corresponding CPU time is lower compared to TV and the traditional Perona-Malik (PM) model. And although D-α-PM model shows better results than PSNR and MAE results, it can be observed that in terms of edge feature recovery measure (PSNRE) and general structural coincidence measure (SSIM) our method performs better. Moreover, implementing our model by the OAS scheme revealed faster execution in terms of both the CPU time and iteration steps (see Table 1). The MAE has lower and even better PSNR results than TV (see Table 1). Looking at visual results of our method using PM scheme (PMS) and AOS, respectively, in Figures 1(c) and 1(d), there is manifestly better visual appeal for our method compared to the PM method (see Figure 3(e)) which shows some speckles, TV method (see Figure 1(f)) which exhibits staircase effects and slight loss in contrast, and the D-α-method (see Figure 3(g)) which shows slightly deformed edges.However, for real image such as the given Lena image in Figure3, the results as shown in Table 2 indicate that our method by AOS scheme gives superior results as shown by the extremely lower iteration steps (7 steps), very short CPU processing time (1.24 sec), better PSNR (28.06), and even lower MAE (7.27) compared to those of the TV method and PM method. Note that our model implemented using AOS performs better than the same model implemented using PM scheme (PMS) for real images. And, in spite of the slightly superior PSNR and MAE values by the D-α-PM method, our method not only gives better edge preservation as evidenced by the higher PSNRE, but also gives closer structural coincidence than the other three models. Further, the visual appeal of our method, whether using PMS (Figure 3(c)) or AOS scheme (Figure 3(d)), excels those of the traditional Perona-Malik (PM) (see Figure 3(e)) method which shows speckles and a bit of blur, TV method (see Figure 3(f)) which introduces staircasing effects on the denoised image, and even the D-α-PM method, which shows blockiness and some slight speckle effects.In addition, from the similarity curves given in Figures2 and 4, it can be observed in the light of the marked areas, for instance, that this model performs better than its comparisons, both for the synthetic image and real image. Note also that, even though in performance metrics, especially in terms of PSNR and MAE, the D-α-PM method seems to perform better, the similarity curves attest that our method generates restored images that more closely match the original image than the results obtained by the D-α-PM method (see Figures 2(e) and 4(e)).Synthetic image (300
×
300). Similarity graphs between the original image and results of our method (AOS and PMS), PM method, TV method, and D-α-PM method, respectively.
(a)
Proposed model I (AOS)
(b)
Proposed model I (PMS)
(c)
PM model
(d)
TV model
(e)
D-α-PM modelLena image (300
×
300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ
=
30. (c) Our algorithm by PMS; τ
=
0.2 (122 steps). (d) Our algorithm by AOS scheme; τ
=
2 (7 steps). (e) PM method; K
=
12 and τ
=
0.25 (60 steps). (f) TV method; τ
=
0.2 (149 steps). (g) D-α-PM method; σ
1
=
0.5, σ
=
30, τ
=
2, and K
=
4 (13 steps).
(a)
Original image
(b)
Noisy image
(c)
Proposed model I (PMS)
(d)
Proposed model I (AOS)
(e)
PM model
(f)
TV model
(g)
D-α-PM modelLena image (300
×
300). Similarity graphs between the original image and results of our method (AOS and PMS), PM method, TV method, and D-α-PM method, respectively.
(a)
Proposed model I (AOS)
(b)
Proposed model I (PMS)
(c)
PM model
(d)
TV model
(e)
D-α-PM model


## 5.1. Numerical Scheme

In the following two subsections, two numerical discrete schemes, the PM scheme (PMS) and the additive operator splitting (AOS) scheme, have been proposed.


### 5.1.1. PM Scheme

Here, we have proposed a numerical scheme similar to the original PM method, whereby (12)–(14) are discretized as follows:
(81)
C
N
i
,
j
n
=
K
1
+
K
|
∇
N
u
i
,
j
|
,
C
S
i
,
j
n
=
K
1
+
K
|
∇
S
u
i
,
j
|
,
C
E
i
,
j
n
=
K
1
+
K
|
∇
E
u
i
,
j
|
,
C
W
i
,
j
n
=
K
1
+
K
|
∇
W
u
i
,
j
|
,
di
v
i
,
j
n
=
[
C
N
i
,
j
n
∇
N
u
i
,
j
+
C
S
i
,
j
n
∇
S
u
i
,
j
+
C
E
i
,
j
n
∇
E
u
i
,
j
+
C
W
i
,
j
n
∇
W
u
i
,
j
]
,
and λ is dynamically determined according to the following discretization scheme:
(82)
λ
n
=
1
σ
2
|
Ω
|
∑
i
,
j
‍
di
v
i
,
j
n
(
u
i
,
j
-
f
i
,
j
)
,
where
|
Ω
|
=
M
N
is
the
size
of
image
.
Hence from (81) and (82) we have
(83)
u
i
,
j
n
+
1
=
u
i
,
j
n
+
τ
di
v
i
,
j
n
-
λ
n
τ
(
u
i
,
j
-
f
i
,
j
)
,
u
i
,
j
0
=
f
i
,
j
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
M
,
j
n
=
u
M
-
1
,
j
n
,
u
i
,
N
n
=
u
i
,
N
-
1
n
,
where
(84)
∇
N
u
i
,
j
=
u
i
-
1
,
j
-
u
i
,
j
,
∇
S
u
i
,
j
=
u
i
+
1
,
j
-
u
i
,
j
,
∇
E
u
i
,
j
=
u
i
,
j
+
1
-
u
i
,
j
,
∇
W
u
i
,
j
=
u
i
,
j
-
1
-
u
i
,
j
,
for i
=
0,1
,
2
,
…
,
N and j
=
0,1
,
2
,
…
,
M.


### 5.1.2. AOS Scheme

In this part, using a AOS scheme, the problem (12)–(14) has been discretized as follows:
(85)
λ
0
=
0
,
u
n
+
1
=
1
m
∑
l
=
1
m
‍
[
I
-
m
τ
A
l
(
u
k
)
]
-
1
[
u
n
+
λ
τ
(
f
-
u
n
)
]
,
div
⁡
n
=
(
u
n
+
1
-
u
n
)
τ
,
λ
n
=
1
σ
2
M
N
(
u
-
f
)
div
⁡
n
,
u
i
,
j
0
=
f
i
,
j
=
f
(
i
h
,
j
h
)
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
I
,
i
n
=
u
I
-
1
,
i
n
,
u
i
,
J
n
=
u
i
,
J
-
1
n
,
where A
l
(
u
n
)
=
[
a
i
,
j
(
u
n
)
],
(86)
a
i
,
j
(
u
n
)
:
=
{
C
i
n
+
C
j
n
2
h
2
,
[
j
∈
N
(
i
)
]
,
-
∑
N
∈
N
(
i
)
C
i
n
+
C
N
n
2
h
2
‍
,
(
j
=
i
)
,
0
,
(
else
)
,
C
i
n
:
=
K
1
+
K
|
∇
u
i
,
j
n
|
,
where
(87)
|
∇
u
i
,
j
n
|
=
1
2
∑
p
,
q
∈
N
(
i
)
|
u
p
n
-
u
q
n
|
2
h
‍
,
where N
(
i
) is the set of the two neighbors of pixel i (boundary pixels have only one neighbor).It is observed that AOS schemes with large time steps still reveal average grey value invariance, stability based on extremum principle, Lyapunov functionals, and convergence to a constant steady state [10]. The AOS scheme is less than twice the typical effort needed for the PM scheme, a rather low price for gaining absolute stability [10]. It is worth noting that if we use div
⁡
(
∇
u
/
ε
+
|
∇
u
|
2
) to approximate div
⁡
(
∇
u
/
|
∇
u
|
) (TV kernel), in numerical scheme, the AOS scheme may be unstable because of the small number ε. Our approximation, however, effectively avoids instabilities arising from such a scenario, as there is no need for lifting the denominator, since it cannot not be equal to zero. This has been an additional motivation for considering AOS for our numerical experiments. The possibility of occurrence of a zero denominator in the evolution problem is a phenomenon that makes the numerical implementation of TV problematic.


## 5.1.1. PM Scheme

Here, we have proposed a numerical scheme similar to the original PM method, whereby (12)–(14) are discretized as follows:
(81)
C
N
i
,
j
n
=
K
1
+
K
|
∇
N
u
i
,
j
|
,
C
S
i
,
j
n
=
K
1
+
K
|
∇
S
u
i
,
j
|
,
C
E
i
,
j
n
=
K
1
+
K
|
∇
E
u
i
,
j
|
,
C
W
i
,
j
n
=
K
1
+
K
|
∇
W
u
i
,
j
|
,
di
v
i
,
j
n
=
[
C
N
i
,
j
n
∇
N
u
i
,
j
+
C
S
i
,
j
n
∇
S
u
i
,
j
+
C
E
i
,
j
n
∇
E
u
i
,
j
+
C
W
i
,
j
n
∇
W
u
i
,
j
]
,
and λ is dynamically determined according to the following discretization scheme:
(82)
λ
n
=
1
σ
2
|
Ω
|
∑
i
,
j
‍
di
v
i
,
j
n
(
u
i
,
j
-
f
i
,
j
)
,
where
|
Ω
|
=
M
N
is
the
size
of
image
.
Hence from (81) and (82) we have
(83)
u
i
,
j
n
+
1
=
u
i
,
j
n
+
τ
di
v
i
,
j
n
-
λ
n
τ
(
u
i
,
j
-
f
i
,
j
)
,
u
i
,
j
0
=
f
i
,
j
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
M
,
j
n
=
u
M
-
1
,
j
n
,
u
i
,
N
n
=
u
i
,
N
-
1
n
,
where
(84)
∇
N
u
i
,
j
=
u
i
-
1
,
j
-
u
i
,
j
,
∇
S
u
i
,
j
=
u
i
+
1
,
j
-
u
i
,
j
,
∇
E
u
i
,
j
=
u
i
,
j
+
1
-
u
i
,
j
,
∇
W
u
i
,
j
=
u
i
,
j
-
1
-
u
i
,
j
,
for i
=
0,1
,
2
,
…
,
N and j
=
0,1
,
2
,
…
,
M.


## 5.1.2. AOS Scheme

In this part, using a AOS scheme, the problem (12)–(14) has been discretized as follows:
(85)
λ
0
=
0
,
u
n
+
1
=
1
m
∑
l
=
1
m
‍
[
I
-
m
τ
A
l
(
u
k
)
]
-
1
[
u
n
+
λ
τ
(
f
-
u
n
)
]
,
div
⁡
n
=
(
u
n
+
1
-
u
n
)
τ
,
λ
n
=
1
σ
2
M
N
(
u
-
f
)
div
⁡
n
,
u
i
,
j
0
=
f
i
,
j
=
f
(
i
h
,
j
h
)
,
u
i
,
0
n
=
u
i
,
1
n
,
u
0
,
j
n
=
u
1
,
j
n
,
u
I
,
i
n
=
u
I
-
1
,
i
n
,
u
i
,
J
n
=
u
i
,
J
-
1
n
,
where A
l
(
u
n
)
=
[
a
i
,
j
(
u
n
)
],
(86)
a
i
,
j
(
u
n
)
:
=
{
C
i
n
+
C
j
n
2
h
2
,
[
j
∈
N
(
i
)
]
,
-
∑
N
∈
N
(
i
)
C
i
n
+
C
N
n
2
h
2
‍
,
(
j
=
i
)
,
0
,
(
else
)
,
C
i
n
:
=
K
1
+
K
|
∇
u
i
,
j
n
|
,
where
(87)
|
∇
u
i
,
j
n
|
=
1
2
∑
p
,
q
∈
N
(
i
)
|
u
p
n
-
u
q
n
|
2
h
‍
,
where N
(
i
) is the set of the two neighbors of pixel i (boundary pixels have only one neighbor).It is observed that AOS schemes with large time steps still reveal average grey value invariance, stability based on extremum principle, Lyapunov functionals, and convergence to a constant steady state [10]. The AOS scheme is less than twice the typical effort needed for the PM scheme, a rather low price for gaining absolute stability [10]. It is worth noting that if we use div
⁡
(
∇
u
/
ε
+
|
∇
u
|
2
) to approximate div
⁡
(
∇
u
/
|
∇
u
|
) (TV kernel), in numerical scheme, the AOS scheme may be unstable because of the small number ε. Our approximation, however, effectively avoids instabilities arising from such a scenario, as there is no need for lifting the denominator, since it cannot not be equal to zero. This has been an additional motivation for considering AOS for our numerical experiments. The possibility of occurrence of a zero denominator in the evolution problem is a phenomenon that makes the numerical implementation of TV problematic.


## 5.2. Comparison with Other Methods

The experiments in this work were performed on a Compaq610 computer, having Intel Core 2 Duo CPU T5870 each 2.00 GHz, physical RAM of 4.00 GB, and Professional Windows 8 64-bit Operating System, on MATLAB R2013b. The image restoration performance was measured in terms of the peak-signal-to-noise ratio (PSNR), mean absolute deviation/error (MAE), structural similarity index measure (SSIM), the measure of similarity of edges (PSNRE), and visual effects. The iteration stopping mechanism was based on the maximal PSNR or the lowest MAE. The PSNR and MAE values were obtained using the formula by Durand et al. [40] and are given by
(88)
PSNR
=
10
log
⁡
10
M
N
|
max
⁡
u
0
-
min
⁡
u
0
|
2
∥
u
-
u
0
∥
L
2
2
d
B
,
MAE
=
∥
u
-
u
0
∥
L
1
M
N
.
Also, taking the edge map EM
(
u
) from the evolution equation, we have
(89)
EM
(
u
)
=
K
(
1
+
K
|
∇
u
|
)
.
And corresponding edge similarity measure is given, according to the formulation by Guo et al. [11], by
(90)
PSNR
E
=
PSNR
(
EM
(
u
)
,
EM
(
u
0
)
)
d
B
.
The SSIM measures have been obtained according to the formula by Wang et al. [41] given by
(91)
SSIM
(
u
,
u
0
)
=
L
(
u
,
u
0
)
·
C
(
u
,
u
0
)
·
R
(
u
,
u
0
)
,
where u
0 denotes the noise-free image, u is the denoised image, M
×
N is the dimension of image, and |

max
⁡
u
0
-
min
⁡
u
0
| yields the gray scale range of the original image. In addition, L
(
u
,
u
0
)
=
(
2
μ
u
μ
u
0
+
k
1
)
/
(
μ
u
2
+
μ
u
0
2
+
k
1
) compares the two images’ mean luminances μ
u and μ
u
0. The maximal value of L
(
u
,
u
0
)
=
1, if μ
u
=
μ
u
0 and C
(
u
,
u
0
)
=
(
2
σ
u
σ
u
0
+
k
2
)
/
(
σ
u
2
+
σ
u
0
2
+
k
2
), measures the closeness of contrast of the two images u and u
0. Contrast is determined in terms of standard deviation, σ. Contrast comparison measure C
(
u
,
u
0
)
=
1 maximally if and only if σ
u
=
σ
u
0, that is, when the images have equal contrast.R
(
u
,
u
0
)
=
(
σ
u
u
0
+
k
3
)
/
(
σ
u
σ
u
0
+
k
3
), where σ
u
u
0 is covariance between u and u
0, is a structure comparison measure which determines the correlation between the images u and u
0. It attains maximal value of 1 if structurally the two images coincide, but its value is equal to zero when there is absolutely no structural coincidence. The quantities k
1, k
2, and k
3 are small positive constants that avert the possibility of having zero denominators.The results of our method were compared to PM method [7], TV method [8], and the D-α-PM method [11]. Tables 1 and 2, respectively, give a summary of the results from the experiments, using synthetic image in Figure 1 and Lena image in Figure 3. The parameters considered here are thresholding parameter K, variance δ, the time step parameter τ, and the convolution parameter σ
1 applied in the D-α-PM method. The fidelity parameter λ was dynamically obtained according to (82) or under the AOS scheme in Section 5.1.2, while the rest of the parameters were chosen to guarantee stability and attainment of optimal results.Table 1
Numerical results for synthetic image(
300
×
300
).
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
PSNRE
SSIM
σ
K
τ
PM
30
12
0.25
232
12.73
34.23
2.84
25.26
0.9829
TV
30
n/a
0.2
203
17.86
32.79
3.76
15.76
0.9748
D-α-PM
30
1
0.25
42
2.40
37.30
3.38
24.76
0.9838
PMS
30
n/a
0.20
262
12.01
35.75
2.15
25.60
0.9870
AOS
30
n/a
3.00
17
2.51
33.30
2.92
25.40
0.9839Table 2
Numerical results for Lena image(
300
×
300
).
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
PSNRE
SSIM
σ
K
τ
PM
30
12
0.25
60
2.43
27.10
7.75
21.70
0.7234
TV
30
n/a
0.2
149
13.87
27.46
7.49
22.63
0.7960
D-α-PM
30
4
2
13
0.90
28.17
7.35
24.78
0.7880
PMS
30
n/a
0.2
122
5.80
27.73
7.45
26.60
0.8060
AOS
30
n/a
2.0
7
1.24
28.06
7.27
27.80
0.8509Synthetic image (300
×
300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ
=
30. (c) Our algorithm by PMS; τ
=
0.2 (262 steps). (d) Our algorithm by AOS scheme; τ
=
3 (17 steps). (e) PM method; K
=
5 and τ
=
0.25 (232 steps). (f) TV method; τ
=
0.2 (203 steps). (g) D-α-PM method; σ
1
=
0.5, σ
=
30, τ
=
0.25, and K
=
1 (42 steps).
(a)
Original image
(b)
Noisy image
(c)
Proposed model I (PMS)
(d)
Proposed model I (AOS)
(e)
PM model
(f)
TV model
(g)
D-α-PM modelFor nontexture images as displayed in Figure1 and the results of our method using PMS scheme, as shown in the fourth row of Table 1, demonstrate better performance than those PM and TV as indicated by the higher PSNR, PSNRE, and SSIM values and lower MAE. And, in spite of higher iterative steps, the corresponding CPU time is lower compared to TV and the traditional Perona-Malik (PM) model. And although D-α-PM model shows better results than PSNR and MAE results, it can be observed that in terms of edge feature recovery measure (PSNRE) and general structural coincidence measure (SSIM) our method performs better. Moreover, implementing our model by the OAS scheme revealed faster execution in terms of both the CPU time and iteration steps (see Table 1). The MAE has lower and even better PSNR results than TV (see Table 1). Looking at visual results of our method using PM scheme (PMS) and AOS, respectively, in Figures 1(c) and 1(d), there is manifestly better visual appeal for our method compared to the PM method (see Figure 3(e)) which shows some speckles, TV method (see Figure 1(f)) which exhibits staircase effects and slight loss in contrast, and the D-α-method (see Figure 3(g)) which shows slightly deformed edges.However, for real image such as the given Lena image in Figure3, the results as shown in Table 2 indicate that our method by AOS scheme gives superior results as shown by the extremely lower iteration steps (7 steps), very short CPU processing time (1.24 sec), better PSNR (28.06), and even lower MAE (7.27) compared to those of the TV method and PM method. Note that our model implemented using AOS performs better than the same model implemented using PM scheme (PMS) for real images. And, in spite of the slightly superior PSNR and MAE values by the D-α-PM method, our method not only gives better edge preservation as evidenced by the higher PSNRE, but also gives closer structural coincidence than the other three models. Further, the visual appeal of our method, whether using PMS (Figure 3(c)) or AOS scheme (Figure 3(d)), excels those of the traditional Perona-Malik (PM) (see Figure 3(e)) method which shows speckles and a bit of blur, TV method (see Figure 3(f)) which introduces staircasing effects on the denoised image, and even the D-α-PM method, which shows blockiness and some slight speckle effects.In addition, from the similarity curves given in Figures2 and 4, it can be observed in the light of the marked areas, for instance, that this model performs better than its comparisons, both for the synthetic image and real image. Note also that, even though in performance metrics, especially in terms of PSNR and MAE, the D-α-PM method seems to perform better, the similarity curves attest that our method generates restored images that more closely match the original image than the results obtained by the D-α-PM method (see Figures 2(e) and 4(e)).Synthetic image (300
×
300). Similarity graphs between the original image and results of our method (AOS and PMS), PM method, TV method, and D-α-PM method, respectively.
(a)
Proposed model I (AOS)
(b)
Proposed model I (PMS)
(c)
PM model
(d)
TV model
(e)
D-α-PM modelLena image (300
×
300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ
=
30. (c) Our algorithm by PMS; τ
=
0.2 (122 steps). (d) Our algorithm by AOS scheme; τ
=
2 (7 steps). (e) PM method; K
=
12 and τ
=
0.25 (60 steps). (f) TV method; τ
=
0.2 (149 steps). (g) D-α-PM method; σ
1
=
0.5, σ
=
30, τ
=
2, and K
=
4 (13 steps).
(a)
Original image
(b)
Noisy image
(c)
Proposed model I (PMS)
(d)
Proposed model I (AOS)
(e)
PM model
(f)
TV model
(g)
D-α-PM modelLena image (300
×
300). Similarity graphs between the original image and results of our method (AOS and PMS), PM method, TV method, and D-α-PM method, respectively.
(a)
Proposed model I (AOS)
(b)
Proposed model I (PMS)
(c)
PM model
(d)
TV model
(e)
D-α-PM model


## 6. Conclusion

In this paper, we have proposed a modified total variation model based on the strictly convex modification for image denoising. The main idea was to offer a better image restoration model that is strictly convex and is, therefore, not subject to backward diffusion which has the potential of introducing blurs and to limit the number of parameters available for manual manipulation. This, probably, explains why, in spite of the fact that our model tends to TV asK tends to +
∞, its practical implementation enjoys better image visual sharpness. This visual sharpness is comparably visible in the results by PM method, but it is still marred behind the oversmoothing and speckle effects in PM images. In fact, the reason why the D-α-PM method tends to give images that have some blur and deformed edges is attributable to the backward diffusion that it introduces in the course of diffusion. Indeed, it is difficult to assure convergence to a solution of minimum energy given the backward-forward motion of the diffusion process induced by the D-α-PM algorithm. At the practical level, the increased number of parameters in the D-α-PM algorithm tends to make it difficult to arrive at a definite permutation of parameter values that would give an optimal result.For the proposed method, we have demonstrated the existence and uniqueness of the solution of the model. Moreover, numerical implementation of our model also does not suffer potential inaccuracies typical of the TV method, since we do not need to add any small perturbation constant in TV method. Our thresholding parameterK depends on the evolution parameter t and therefore does not have to be constrained to smaller values as in the case of PM model. The resultant evolution equation has been discretized and implemented using PM scheme and AOS scheme to demonstrate the performance of our algorithm. From the given experimental results, the PSNR values, MAE values, Iterative steps, the CPU processing time, PSNRE, SSIM, and visual appeal of our denoised images all testify that our method is actually a good balance of the PM, TV, and D-α-PM methods and hence a better image restoration model.However, in real life application, we observe that the success of any denoising algorithm depends on the type of image being considered, the type of noise, the application intended for the results of the restoration, the extent of degradation, and indeed the implementation platform. And although we have only considered additive Gaussian white noise, different kinds of noise (whether Poisson, speckle, salt, or paper, among others) will require different formulations for the fidelity part and may even demand the application of more than just one formulation for effectiveness. The choice of platform and scheme of implementation must also be appropriately made for efficient performance of the formulation.With respect to the use of the restoration results, there are situations where the noise removal, generally, may be counterproductive. This usually occurs when the oscillations due to the noise are of comparable scale to those of the features being targeted for preservation. A case like this may require a combination of formulations [42, 43].For images that are heavily degraded, it might be necessary to do a preconvolution, to blur the noise effects, and then to use an effective formulation such as this one to recover semantically important features such as the edges and contours of the image.



---

*Source: 948392-2014-06-18.xml*