# An Alternative Variational Framework for Image Denoising

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

---

## Abstract

We propose an alternative framework for total variation based image denoising models. The model is based on the minimization of the total variation with a functional coefficient, where, in this case, the functional coefficient is a function of the magnitude of image gradient. We determine the considerations to bear on the choice of the functional coefficient. With the use of an example functional, we demonstrate the effectiveness of a model chosen based on the proposed consideration. In addition, for the illustrative model, we prove the existence and uniqueness of the minimizer of the variational problem. The existence and uniqueness of the solution associated evolution equation are also established. Experimental results are included to demonstrate the effectiveness of the selected model in image restoration over the traditional methods of Perona-Malik (PM), total variation (TV), and the D-α-PM method.


---

## Body

## 1. Introduction

The objective of any image restoration process should not focus only on the removal of noise, but it should also observe, as Perona and Malik [1], Koenderink [2], and Witkin [3] determined, that no new spurious details are created in the restored image; at each scale-space representation, the boundaries/edges are sharp or preserved, and, at all scales, intraregion smoothing is preferred to interregion smoothing.In the light of the above considerations, researchers have observed that it is then logical to obtain or develop edge indicators that would be adapted to the local image structure [4]. Consequently, a number of edge indicators have been proposed and logically grafted into the partial differential equation (PDE) based evolution equations [1, 5]. In addition, energy minimization problems are continually being formulated which focus on producing adaptive partial differential equations [6, 7].Total variation (TV) method, widely considered a powerful technique for smoothing, edge preservation, and general image restoration, was first proposed by Rudin et al.  [8]. The method is based on the strength of the argument that TV norms principally are L1-norms of derivatives and that L1-norms provide the proper basis for image restoration [8, 9]. 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 for “jumps” or discontinuities and thus be able to protect edges.TheL2-norm of the gradient of image allows removing noise; however, it has the adverse effect of penalizing too much the gradients corresponding to edges [10]. A functional based on the L2-norm does not permit discontinuities in the solution, and thus edges cannot be recovered properly [11].The total variation (TV) norm in [8] is a regularization functional of the form
(1)TV(u)=∫Ω‍|∇u|dx=∫Ωux2+uy2dx.
Although TV regularization above allows for edge recovery, it has some demerits. Firstly, the formulation favours solutions which are piecewise constant. This has the effect of causing staircase effects and may even generate false edges on the image [12]. Secondly, the method has the effect of reducing contrast even in regions of the same pixel intensity or in noise-free observed images [13, 14].Various modifications have, therefore, been proposed, in an attempt to address the drawbacks of the TV model, to make it as adaptive as possible to the local image structure. For instance, Strong and Chan in [15] proposed the spatially adaptive regularization functional of the form
(2)J(u)=∫Ω‍α(x)|∇u|dx,
where the control factor α(x) is designed to slow diffusion in the neighbourhood of edges. Then, we have the efforts, according to Blomgren et al. in [13], which produced a denoising functional of the form
(3)J(u)=∫Ω‍|∇u|p(|∇u|),
where p is a nondecreasing function of the magnitude of gradient, p(s)→2 as s→0 and p(s)→1 as s→∞. This model is designed to automatically tap into the benefits of both isotropic diffusion and TV regularization. For other modifications, we refer the reader to, among others, the works in [12, 16, 17].In this paper, however, we propose an alternative framework of variational model for image denoising, where the regularization potential is a product of a gradient based functional coefficient and the norm of gradient of image (potential function for TV); namely,F(u)=∫Ω‍Ψ(|∇u|)·|∇u|dx. That is, the coefficient function is a function of the magnitude of the gradient of the image. We propose the criteria for the choice of the coefficient. Finally, we have selected a model based on the proposed criteria, as an example for further analysis and demonstration of experimental results.The structure of this paper is as follows. In Section2, we present the proposed model (4) and discuss the general criteria for choosing the functional coefficient to the traditional TV potential. Section 2 is concluded by considering, based on the proposed criteria, an example functional for further analysis. In Section 3, we give certain preliminary definitions and lemmas we rely on, variously, in this paper. In addition, we prove the existence and uniqueness of the solution to the minimization problem (11) and the associated evolution equation (27)–(29). In Section 4, we define the weak solution to the evolution problem (27)–(11). Furthermore, we present the formulation of the approximate evolution equation (32)–(36), prove the existence of solutions to the approximate evolution problem, and conclude Section 4 with the existence and uniqueness of the solution to the evolution problem (27)–(29). In Section 5, we give the numerical schemes and experimental results to demonstrate the strength and effectiveness of our method. Additionally, we have presented a brief discussional comparison of our results with those of other methods like PM, original TV, and D-α-PM method by Guo et al. [5]. A brief summary concludes the paper in Section 6.


## 2. Proposed Model

In this section, inspired by the works of Chan and Shen [18], Vese [19], and Chen et al. [12], among others, we propose an alternative framework for total variation based denoising model. The model is based on minimization of a functional, where the regularization potential is a product of the total variation potential and a gradient based functional coefficient. The coefficient, which acts to penalize the norm of gradient and detect any edges, is also taken as function of the norm of the gradient of image. So, we present the general form of the model and determine certain properties of the functional coefficient, especially in terms of linearity, sublinearity, and superlinearity growth at infinity. Ultimately, we have selected a specific example of model, based on the proposed guidelines, for further analysis.


### 2.1. The New Framework for Energy Functional

The proposed energy functional is given in the following general form:(4)min⁡u∈BV∩L2(Ω){F(u)=∫Ω‍Ψ(|∇u|)|∇u|dx+λ2∫Ω‍(u-f)2dx},
where Ψ(s) is a function that detects edges and penalizes the norm of gradient and f is the noise image. Observe that if we set Ψ(s)=s, then the kernel of F(u) becomes a variational problem of the L2-norm that is known not to allow for discontinuities, since it leads to the traditional edge obliterating isotropic diffusion [20], while if Ψ(s)=1 the problem becomes the usual TV functional, which allows a diffusion mechanism that is strictly normal to the image gradient [12].The functionalF(u) in (4) should verify the Euler-Lagrange equation
(5)0=-div⁡{(Ψ′(|∇u|)+Ψ(|∇u|)|∇u|)∇u}+λ(u-f).
Next, we establish certain properties of Ψ(s) from (5) that should influence its choice. For the purposes of denoising, to be able to smooth the image in homogeneous regions, that is, regions where s=|∇u|→0,Ψ(s) is expected such that lim⁡s→0⁡(Ψ′(|∇u|)+(Ψ(|∇u|)/|∇u|))=C(constant)>0. On the other hand, to preserve discontinuities (edges) in the image, that is, in areas where s=|∇u|→+∞, Ψ(s) takes a form such that lim⁡s→+∞⁡(Ψ′(|∇u|)+(Ψ(|∇u|)/|∇u|))=0.Now, for a little more precision on the expected behaviour ofΨ(s), we decompose divergence part in terms of the tangent and normal directions to the isophote lines. Equation (5) then becomes
(6)0=[Ψ′(|∇u|)+Ψ(|∇u|)|∇u|]uTT+[2Ψ′(|∇u|)+Ψ′′(|∇u|)|∇u|]uNN-λ(u-f).
Then, to be able to achieve isotropic diffusion within homogeneous regions, from (6), we may impose the condition that
(7)lim⁡s→0[Ψ′(s)+Ψ(s)s]=lim⁡s→0[2Ψ′(s)+Ψ′′(s)·s]=C>0,
where C is a constant. The above condition yields from (6) a diffusion model that diffuses isotropically.In the regions neighboring the edges, the model should dissipate diffusion effects across the edges and favor diffusion along the edges. This implies that ass→∞,  [Ψ′(s)+(Ψ(s)/s)]→C1, where C1 is an arbitrary positive constant, while [2Ψ′(s)+Ψ′′(s)·s]→0. These conditions might be difficult to meet simultaneously for Ψ(s). Therefore, motivated by the works in [10, 19], a weaker compromise is imposed by demanding that both terms approach zero, but at different rates, with diffusion along the direction normal to the isophote lines approaching zero faster than the diffusion along the tangent to isophote lines. This leads to the condition that
(8)lim⁡s→+∞[2Ψ′(s)+Ψ′′(s)·s][Ψ′(s)+(Ψ(s)/s)]=0.Hence, we observe that, forΨ(s) to become an effective coefficient in the functional in (4), it should be such that conditions (7) and (8) are satisfied.For simplicity of notation in subsequent stages, we will denote, from (5),
(9)φ′(∇u):=Ψ′(|∇u|)|∇u|+Ψ(|∇u|),so thatφ(∇u):=Ψ(|∇u|)|∇u|.Furthermore, we are interested in the properties ofΨ(s) with regard to its growth at infinity and its overall impact on the functional F(u):
(10)iflim⁡s→+∞Ψ(s)s={C,thenF(u)is superlinear,0,thenF(u)is eithersuperlinear or linear,+∞,thenF(u)is superlinear.
Note that C>0 is a constant and growth.


#### 2.1.1. Summary on the Characterization of Growth ofΨ(·)

The aim of this section is to discuss the property of the function coefficientΨ(s) that will lead to a functional F(u) of linear growth at infinity, but, of course, subject to conditions (7)-(8). The characterization is based upon the results of process (10).Observe that setting the coefficientΨ(s) in the regularization kernel, as a function of the norm of gradient of the image to be linear growth, leads to a functional F(u) in (4), which is of superlinear growth. Such a functional yields associated diffusion equation which either diffuses images uniformly, without being sensitive to discontinuities or edges, or generally does not yield good denoising results [10, 21].ChoosingΨ(s) that is of superlinear growth leads to F(u) which is invariably of superlinear growth. For instance, if we set Ψ(s)=sp,(p>1) which is of superlinear growth at infinity, then we obtain the function Ψ(s)·s=sp+1. This functional is of superlinear growth and does not give good results in image denoising. This is because, with such a choice of Ψ(s), the derivative of Ψ(s)·s would yield a nondegenerate elliptic differential operator of the second order, which would have an oversmoothing effect under the optimality condition [22].Bildhauer and Fuchs [21] reckon that a regularization function need not necessarily be of power growth. Therefore, another example for Ψ(s) of superlinear growth is Ψ(s)·s=s2log⁡(1+s), where Ψ(s)=slog⁡(s+1) is viewed as some kind of compromise between the cases p>1 and p=1 of coefficient of power growth [21]. This formulation too does not promise any better results.Consequently, we observe that to achieve maximal results, which entails reconstructing a noise image in such a way that the edges are protected while diffusion within homogeneous regions proceeds rather uniformly, it is required thatΨ(s) be of sublinear growth at infinity. Also, the chosen Ψ(s) should be such that the functional F(u) generated is of linear growth. Moreover, Ψ(s) should be such that the conditions in (7)-(8) are satisfied. Functionals which arise from such potential functions tend to give better image denoising results. Some examples of the functions Ψ(s) of sublinear growth are Ψ(s):=s/(1+s),Ψ(s):=log⁡(1+s), and Ψ(s):=1/(1+αs). The function Ψ(s):=s/(1+s) is a good candidate for further consideration.In this paper, therefore, we choose the control functionΨ(s):=s/(1+s), which is function of sublinear growth at infinity. This leads to a nondecreasing potential function for s∈ℝ+, which is strictly convex and of linear growth at infinity given in the form Ψ(s)·s:=s2/(1+s).With this choice, we observe the following further properties ofΨ(s);(i)
lim⁡s→0⁡Ψ(s)=0;(ii)
lim⁡s→∞⁡Ψ(s)=1;(iii)
Ψ(s) is nondecreasing for s≥0. This, as observed from (6), ensures only forward diffusion along the tangent to the isophote lines.


## 2.1. The New Framework for Energy Functional

The proposed energy functional is given in the following general form:(4)min⁡u∈BV∩L2(Ω){F(u)=∫Ω‍Ψ(|∇u|)|∇u|dx+λ2∫Ω‍(u-f)2dx},
where Ψ(s) is a function that detects edges and penalizes the norm of gradient and f is the noise image. Observe that if we set Ψ(s)=s, then the kernel of F(u) becomes a variational problem of the L2-norm that is known not to allow for discontinuities, since it leads to the traditional edge obliterating isotropic diffusion [20], while if Ψ(s)=1 the problem becomes the usual TV functional, which allows a diffusion mechanism that is strictly normal to the image gradient [12].The functionalF(u) in (4) should verify the Euler-Lagrange equation
(5)0=-div⁡{(Ψ′(|∇u|)+Ψ(|∇u|)|∇u|)∇u}+λ(u-f).
Next, we establish certain properties of Ψ(s) from (5) that should influence its choice. For the purposes of denoising, to be able to smooth the image in homogeneous regions, that is, regions where s=|∇u|→0,Ψ(s) is expected such that lim⁡s→0⁡(Ψ′(|∇u|)+(Ψ(|∇u|)/|∇u|))=C(constant)>0. On the other hand, to preserve discontinuities (edges) in the image, that is, in areas where s=|∇u|→+∞, Ψ(s) takes a form such that lim⁡s→+∞⁡(Ψ′(|∇u|)+(Ψ(|∇u|)/|∇u|))=0.Now, for a little more precision on the expected behaviour ofΨ(s), we decompose divergence part in terms of the tangent and normal directions to the isophote lines. Equation (5) then becomes
(6)0=[Ψ′(|∇u|)+Ψ(|∇u|)|∇u|]uTT+[2Ψ′(|∇u|)+Ψ′′(|∇u|)|∇u|]uNN-λ(u-f).
Then, to be able to achieve isotropic diffusion within homogeneous regions, from (6), we may impose the condition that
(7)lim⁡s→0[Ψ′(s)+Ψ(s)s]=lim⁡s→0[2Ψ′(s)+Ψ′′(s)·s]=C>0,
where C is a constant. The above condition yields from (6) a diffusion model that diffuses isotropically.In the regions neighboring the edges, the model should dissipate diffusion effects across the edges and favor diffusion along the edges. This implies that ass→∞,  [Ψ′(s)+(Ψ(s)/s)]→C1, where C1 is an arbitrary positive constant, while [2Ψ′(s)+Ψ′′(s)·s]→0. These conditions might be difficult to meet simultaneously for Ψ(s). Therefore, motivated by the works in [10, 19], a weaker compromise is imposed by demanding that both terms approach zero, but at different rates, with diffusion along the direction normal to the isophote lines approaching zero faster than the diffusion along the tangent to isophote lines. This leads to the condition that
(8)lim⁡s→+∞[2Ψ′(s)+Ψ′′(s)·s][Ψ′(s)+(Ψ(s)/s)]=0.Hence, we observe that, forΨ(s) to become an effective coefficient in the functional in (4), it should be such that conditions (7) and (8) are satisfied.For simplicity of notation in subsequent stages, we will denote, from (5),
(9)φ′(∇u):=Ψ′(|∇u|)|∇u|+Ψ(|∇u|),so thatφ(∇u):=Ψ(|∇u|)|∇u|.Furthermore, we are interested in the properties ofΨ(s) with regard to its growth at infinity and its overall impact on the functional F(u):
(10)iflim⁡s→+∞Ψ(s)s={C,thenF(u)is superlinear,0,thenF(u)is eithersuperlinear or linear,+∞,thenF(u)is superlinear.
Note that C>0 is a constant and growth.


### 2.1.1. Summary on the Characterization of Growth ofΨ(·)

The aim of this section is to discuss the property of the function coefficientΨ(s) that will lead to a functional F(u) of linear growth at infinity, but, of course, subject to conditions (7)-(8). The characterization is based upon the results of process (10).Observe that setting the coefficientΨ(s) in the regularization kernel, as a function of the norm of gradient of the image to be linear growth, leads to a functional F(u) in (4), which is of superlinear growth. Such a functional yields associated diffusion equation which either diffuses images uniformly, without being sensitive to discontinuities or edges, or generally does not yield good denoising results [10, 21].ChoosingΨ(s) that is of superlinear growth leads to F(u) which is invariably of superlinear growth. For instance, if we set Ψ(s)=sp,(p>1) which is of superlinear growth at infinity, then we obtain the function Ψ(s)·s=sp+1. This functional is of superlinear growth and does not give good results in image denoising. This is because, with such a choice of Ψ(s), the derivative of Ψ(s)·s would yield a nondegenerate elliptic differential operator of the second order, which would have an oversmoothing effect under the optimality condition [22].Bildhauer and Fuchs [21] reckon that a regularization function need not necessarily be of power growth. Therefore, another example for Ψ(s) of superlinear growth is Ψ(s)·s=s2log⁡(1+s), where Ψ(s)=slog⁡(s+1) is viewed as some kind of compromise between the cases p>1 and p=1 of coefficient of power growth [21]. This formulation too does not promise any better results.Consequently, we observe that to achieve maximal results, which entails reconstructing a noise image in such a way that the edges are protected while diffusion within homogeneous regions proceeds rather uniformly, it is required thatΨ(s) be of sublinear growth at infinity. Also, the chosen Ψ(s) should be such that the functional F(u) generated is of linear growth. Moreover, Ψ(s) should be such that the conditions in (7)-(8) are satisfied. Functionals which arise from such potential functions tend to give better image denoising results. Some examples of the functions Ψ(s) of sublinear growth are Ψ(s):=s/(1+s),Ψ(s):=log⁡(1+s), and Ψ(s):=1/(1+αs). The function Ψ(s):=s/(1+s) is a good candidate for further consideration.In this paper, therefore, we choose the control functionΨ(s):=s/(1+s), which is function of sublinear growth at infinity. This leads to a nondecreasing potential function for s∈ℝ+, which is strictly convex and of linear growth at infinity given in the form Ψ(s)·s:=s2/(1+s).With this choice, we observe the following further properties ofΨ(s);(i)
lim⁡s→0⁡Ψ(s)=0;(ii)
lim⁡s→∞⁡Ψ(s)=1;(iii)
Ψ(s) is nondecreasing for s≥0. This, as observed from (6), ensures only forward diffusion along the tangent to the isophote lines.


## 2.1.1. Summary on the Characterization of Growth ofΨ(·)

The aim of this section is to discuss the property of the function coefficientΨ(s) that will lead to a functional F(u) of linear growth at infinity, but, of course, subject to conditions (7)-(8). The characterization is based upon the results of process (10).Observe that setting the coefficientΨ(s) in the regularization kernel, as a function of the norm of gradient of the image to be linear growth, leads to a functional F(u) in (4), which is of superlinear growth. Such a functional yields associated diffusion equation which either diffuses images uniformly, without being sensitive to discontinuities or edges, or generally does not yield good denoising results [10, 21].ChoosingΨ(s) that is of superlinear growth leads to F(u) which is invariably of superlinear growth. For instance, if we set Ψ(s)=sp,(p>1) which is of superlinear growth at infinity, then we obtain the function Ψ(s)·s=sp+1. This functional is of superlinear growth and does not give good results in image denoising. This is because, with such a choice of Ψ(s), the derivative of Ψ(s)·s would yield a nondegenerate elliptic differential operator of the second order, which would have an oversmoothing effect under the optimality condition [22].Bildhauer and Fuchs [21] reckon that a regularization function need not necessarily be of power growth. Therefore, another example for Ψ(s) of superlinear growth is Ψ(s)·s=s2log⁡(1+s), where Ψ(s)=slog⁡(s+1) is viewed as some kind of compromise between the cases p>1 and p=1 of coefficient of power growth [21]. This formulation too does not promise any better results.Consequently, we observe that to achieve maximal results, which entails reconstructing a noise image in such a way that the edges are protected while diffusion within homogeneous regions proceeds rather uniformly, it is required thatΨ(s) be of sublinear growth at infinity. Also, the chosen Ψ(s) should be such that the functional F(u) generated is of linear growth. Moreover, Ψ(s) should be such that the conditions in (7)-(8) are satisfied. Functionals which arise from such potential functions tend to give better image denoising results. Some examples of the functions Ψ(s) of sublinear growth are Ψ(s):=s/(1+s),Ψ(s):=log⁡(1+s), and Ψ(s):=1/(1+αs). The function Ψ(s):=s/(1+s) is a good candidate for further consideration.In this paper, therefore, we choose the control functionΨ(s):=s/(1+s), which is function of sublinear growth at infinity. This leads to a nondecreasing potential function for s∈ℝ+, which is strictly convex and of linear growth at infinity given in the form Ψ(s)·s:=s2/(1+s).With this choice, we observe the following further properties ofΨ(s);(i)
lim⁡s→0⁡Ψ(s)=0;(ii)
lim⁡s→∞⁡Ψ(s)=1;(iii)
Ψ(s) is nondecreasing for s≥0. This, as observed from (6), ensures only forward diffusion along the tangent to the isophote lines.


## 3. The Minimization Problem

Our model based on the presentation above is then given by(11)min⁡u∈BV∩L2(Ω){F(u)=∫Ωφ(∇u)dx+λ2∫Ω(u-f)2dx},
where
(12)φ(∇u):=|∇u|21+|∇u|.
Hence, in this section, we study the existence and uniqueness of the minimization problem (11). The following preliminary information prefaces our reasoning here and in the subsequent sections of this paper.


### 3.1. Preliminaries

Definition 1.
LetΩ be an open subset of ℝn. A function u∈L1(Ω) has a bounded variation in Ω, denoted by BV(Ω), if
(13)∫Ω|Du|=sup⁡{∫Ω‍udiv⁡φdx:φ∈C01(Ω;ℝn),|φ|≤1}<∞,
where BV(Ω) denotes the space of such functions. Then, BV-norm is given by
(14)∥u∥BV=∫Ω‍|∇u|+∥u∥L1(Ω).Definition 2 (see [23]).
Ifu∈BV(Ω), then
(15)Du=∇udx+Dsu.
And writing Du and its total variation on Ω,|Du|(Ω), one has (see [7])
(16)|Du|(Ω)=∫Ω‍|∇u|dx+∫Ω‍Dsu.Du is a Radon measure, where ∇u is the density of the absolutely continuous part of Du with respect to the Lebesgue measure, ℒn, and Dsu is the singular part.Lemma 3 (see [24]).
Letφ:ℝn→ℝ. The function φ is convex if and only if
(17)φ(x)≥φ′(y)(x-y)+φ(y),∀x,y∈ℝn,
and if φ∈C2(ℝ2), then φ is convex if and only if ∇2φ≥0.Lemma 4 (see [10]).
LetΦ:ℝ→ℝ+ be convex, even, and nondecreasing on ℝ+ with linear growth at infinity. Also, let Φ∞ be the recession function of Φ defined by
(18)Φ∞(ω)=lim⁡s→∞Φ(sω)s.
Then, for u∈BV(Ω) and setting Φ(θ)=Φ(|θ|), one has
(19)∫Ω‍Φ(Du)=∫Ω‍Φ(|∇u|)dx+Φ∞(1)∫Ω‍Dsu.
This implies that u→∫Ω‍Φ(Du) is lower semicontinuous for the BV(Ω) topology.


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

Since the objective here is to preserve edges, which are viewed as discontinuities, the natural function space to seek the solution is theBV(Ω) space. In image analysis, the BV(Ω) space allows for discontinuities across the edges of images. Furthermore, due to nonreflexivity of the L1(Ω) space, which would have been a more natural space within which to seek the solution, it is noted that the solution to the minimization problem might not exist there. Therefore, the BV-w*(Ω), which denotes BV-weak* topology, provides the most reasonable alternative space for the existence of the solution. In what follows, we will, for convenience of notation, refer to BV-weak*(Ω) simply as BV(Ω). This space allows obtaining compactness because of the separability of L1(Ω) space even though it is not reflexive [10]. Hence, we state and prove the following existence and uniqueness theorem for our problem.Theorem 5.
Givenφ(∇u) as assigned above, there exists u∈BV(Ω)∩L2(Ω) which is a unique solution to the minimization problem (11).Proof.
Sinceφ(∇u) as specifically defined above for (11) is of linear growth and given the fact thatlim⁡s→∞⁡φ(s)=+∞, it implies that it is also coercive. Consequently, there exists a sequence {un}∈BV(Ω)∩L2(Ω) such that F(un)≤C. Hence, we have ∫Ω‍φ(∇un)<C, and, from the inequality of the second component of F(u), namely, ∫Ω‍(un-f)2dx≤C, if we define wn=(1/|Ω|)∫Ω‍undx, vn=un-wn, we observe that ∫Ω‍vndx=0. We then deduce from [10] that |∫Ω‍undx|≤C and ∥un∥L2(Ω)≤C.
This indicates that{un} is bounded in L2(Ω) and L1(Ω). And from F(un)≤C it is also clear that un is abounded in BV(Ω). Thus, there exists a subsequence {unj}∈BV(Ω) of {un} such that
(20)unj⟶u,weakly inL2(Ω),unj⟶u,strongly inL1(Ω).
And on the strength of the lemma on convex functions of measures (see Lemma 4) and the weak lower semicontinuity of the L2-norm for the second component of F(u), we deduce that F(u) is lower semicontinuous. In fact, the recession function Φ∞(ω) in the definition of convex function of measures is finite for our functional. That is, Φ∞(ω)=1. Hence, we obtain that
(21)∫Ω‍φ(∇u)dx+∫Ω‍(u-f)2dx≤liminf⁡j→∞∫Ω‍φ(∇unj)dx+liminf⁡j→∞∫Ω‍(unj-f)2dx
or
(22)∫Ω‍φ(∇u)dx+∫Ω‍(u-f)2dx≤liminf⁡j→∞{∫Ω‍φ(∇unj)dx+∫Ω‍(unj-f)2dx}.
We, thus, have
(23)F(u)≤liminf⁡j→∞F(unj)=min⁡u∈BV(Ω)∩L2(Ω)F(u).
Thus, u∈BV(Ω)∩L2(Ω) is a minimizer of the problem (11). Uniqueness of u is drawn from the strict convexity of φ(∇u) and convexity of the second component of F(u), which implies the overall strict convexity of the functional F(u). Additionally, a strictly convex functional admits at most one minimum. This, then, implies that u is a unique minimizer of (11).


### 3.3. The Associated Evolution Equation

From the energy minimization problem as assigned above, namely,(24)min⁡u∈BV∩L2(Ω){F(u)=∫Ω‍|∇u|21+|∇u|dx+λ2∫Ω‍(u-f)2dx},
we have the associated Euler-Lagrange equation given by
(25)0=-div⁡(2+|∇u|(1+|∇u|)2∇u)+λ(u-f),x∈Ω,
with the Neumann boundary condition
(26)∂u∂n→=0,x∈∂Ω.
We compute for u in the Euler-Lagrange equation (25) by putting it into a dynamical system, where the time t is used as an evolution parameter. Hence, the evolution equation associated with the minimization problem (4) is given by
(27)∂u∂t=div⁡(2+|∇u|(1+|∇u|)2∇u)-λ(u-f),(x,t)∈QT,(28)u(x,0)=f(x),x∈Ω,(29)∂u∂n→=0,(x,t)∈∂QT,
where f(x) is the noise image, QT:=[0,T]×Ω, and ∂QT:=∂Ω×[0,T].In order to see whether the potentialφ(s):=Ψ(s)·s, as defined above, respects the general principle of image reconstruction, where it is required that reconstructed image be formed by homogeneous regions separated by sharp edges, we decompose the divergence term of (25) using the local image structures like tangent and normal directions to the isophote lines. Writing it in its nonconservative form, analogous to (6), yields
(30)0=-(2+|∇u|(1+|∇u|)2)uTT-(2(1+|∇u|)3)uNN+λ(u-f).
Notice from (30) that as |∇u|→0, signaling homogeneous regions, the potential φ(s) behaves like a linear isotropic diffusion encouraging uniform smoothing in both the uNN and uTT directions. However, as |∇u|→∞, corresponding to the neighbourhood of the edges, diffusion rate along the normal (or uNN) direction is diminished, while the diffusion along the tangent (or uTT) direction is preferred, thereby preserving the edges. The model therefore is well-behaved, since it reasonably satisfies the principle cited above.


## 3.1. Preliminaries

Definition 1.
LetΩ be an open subset of ℝn. A function u∈L1(Ω) has a bounded variation in Ω, denoted by BV(Ω), if
(13)∫Ω|Du|=sup⁡{∫Ω‍udiv⁡φdx:φ∈C01(Ω;ℝn),|φ|≤1}<∞,
where BV(Ω) denotes the space of such functions. Then, BV-norm is given by
(14)∥u∥BV=∫Ω‍|∇u|+∥u∥L1(Ω).Definition 2 (see [23]).
Ifu∈BV(Ω), then
(15)Du=∇udx+Dsu.
And writing Du and its total variation on Ω,|Du|(Ω), one has (see [7])
(16)|Du|(Ω)=∫Ω‍|∇u|dx+∫Ω‍Dsu.Du is a Radon measure, where ∇u is the density of the absolutely continuous part of Du with respect to the Lebesgue measure, ℒn, and Dsu is the singular part.Lemma 3 (see [24]).
Letφ:ℝn→ℝ. The function φ is convex if and only if
(17)φ(x)≥φ′(y)(x-y)+φ(y),∀x,y∈ℝn,
and if φ∈C2(ℝ2), then φ is convex if and only if ∇2φ≥0.Lemma 4 (see [10]).
LetΦ:ℝ→ℝ+ be convex, even, and nondecreasing on ℝ+ with linear growth at infinity. Also, let Φ∞ be the recession function of Φ defined by
(18)Φ∞(ω)=lim⁡s→∞Φ(sω)s.
Then, for u∈BV(Ω) and setting Φ(θ)=Φ(|θ|), one has
(19)∫Ω‍Φ(Du)=∫Ω‍Φ(|∇u|)dx+Φ∞(1)∫Ω‍Dsu.
This implies that u→∫Ω‍Φ(Du) is lower semicontinuous for the BV(Ω) topology.


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

Since the objective here is to preserve edges, which are viewed as discontinuities, the natural function space to seek the solution is theBV(Ω) space. In image analysis, the BV(Ω) space allows for discontinuities across the edges of images. Furthermore, due to nonreflexivity of the L1(Ω) space, which would have been a more natural space within which to seek the solution, it is noted that the solution to the minimization problem might not exist there. Therefore, the BV-w*(Ω), which denotes BV-weak* topology, provides the most reasonable alternative space for the existence of the solution. In what follows, we will, for convenience of notation, refer to BV-weak*(Ω) simply as BV(Ω). This space allows obtaining compactness because of the separability of L1(Ω) space even though it is not reflexive [10]. Hence, we state and prove the following existence and uniqueness theorem for our problem.Theorem 5.
Givenφ(∇u) as assigned above, there exists u∈BV(Ω)∩L2(Ω) which is a unique solution to the minimization problem (11).Proof.
Sinceφ(∇u) as specifically defined above for (11) is of linear growth and given the fact thatlim⁡s→∞⁡φ(s)=+∞, it implies that it is also coercive. Consequently, there exists a sequence {un}∈BV(Ω)∩L2(Ω) such that F(un)≤C. Hence, we have ∫Ω‍φ(∇un)<C, and, from the inequality of the second component of F(u), namely, ∫Ω‍(un-f)2dx≤C, if we define wn=(1/|Ω|)∫Ω‍undx, vn=un-wn, we observe that ∫Ω‍vndx=0. We then deduce from [10] that |∫Ω‍undx|≤C and ∥un∥L2(Ω)≤C.
This indicates that{un} is bounded in L2(Ω) and L1(Ω). And from F(un)≤C it is also clear that un is abounded in BV(Ω). Thus, there exists a subsequence {unj}∈BV(Ω) of {un} such that
(20)unj⟶u,weakly inL2(Ω),unj⟶u,strongly inL1(Ω).
And on the strength of the lemma on convex functions of measures (see Lemma 4) and the weak lower semicontinuity of the L2-norm for the second component of F(u), we deduce that F(u) is lower semicontinuous. In fact, the recession function Φ∞(ω) in the definition of convex function of measures is finite for our functional. That is, Φ∞(ω)=1. Hence, we obtain that
(21)∫Ω‍φ(∇u)dx+∫Ω‍(u-f)2dx≤liminf⁡j→∞∫Ω‍φ(∇unj)dx+liminf⁡j→∞∫Ω‍(unj-f)2dx
or
(22)∫Ω‍φ(∇u)dx+∫Ω‍(u-f)2dx≤liminf⁡j→∞{∫Ω‍φ(∇unj)dx+∫Ω‍(unj-f)2dx}.
We, thus, have
(23)F(u)≤liminf⁡j→∞F(unj)=min⁡u∈BV(Ω)∩L2(Ω)F(u).
Thus, u∈BV(Ω)∩L2(Ω) is a minimizer of the problem (11). Uniqueness of u is drawn from the strict convexity of φ(∇u) and convexity of the second component of F(u), which implies the overall strict convexity of the functional F(u). Additionally, a strictly convex functional admits at most one minimum. This, then, implies that u is a unique minimizer of (11).


## 3.3. The Associated Evolution Equation

From the energy minimization problem as assigned above, namely,(24)min⁡u∈BV∩L2(Ω){F(u)=∫Ω‍|∇u|21+|∇u|dx+λ2∫Ω‍(u-f)2dx},
we have the associated Euler-Lagrange equation given by
(25)0=-div⁡(2+|∇u|(1+|∇u|)2∇u)+λ(u-f),x∈Ω,
with the Neumann boundary condition
(26)∂u∂n→=0,x∈∂Ω.
We compute for u in the Euler-Lagrange equation (25) by putting it into a dynamical system, where the time t is used as an evolution parameter. Hence, the evolution equation associated with the minimization problem (4) is given by
(27)∂u∂t=div⁡(2+|∇u|(1+|∇u|)2∇u)-λ(u-f),(x,t)∈QT,(28)u(x,0)=f(x),x∈Ω,(29)∂u∂n→=0,(x,t)∈∂QT,
where f(x) is the noise image, QT:=[0,T]×Ω, and ∂QT:=∂Ω×[0,T].In order to see whether the potentialφ(s):=Ψ(s)·s, as defined above, respects the general principle of image reconstruction, where it is required that reconstructed image be formed by homogeneous regions separated by sharp edges, we decompose the divergence term of (25) using the local image structures like tangent and normal directions to the isophote lines. Writing it in its nonconservative form, analogous to (6), yields
(30)0=-(2+|∇u|(1+|∇u|)2)uTT-(2(1+|∇u|)3)uNN+λ(u-f).
Notice from (30) that as |∇u|→0, signaling homogeneous regions, the potential φ(s) behaves like a linear isotropic diffusion encouraging uniform smoothing in both the uNN and uTT directions. However, as |∇u|→∞, corresponding to the neighbourhood of the edges, diffusion rate along the normal (or uNN) direction is diminished, while the diffusion along the tangent (or uTT) direction is preferred, thereby preserving the edges. The model therefore is well-behaved, since it reasonably satisfies the principle cited above.


## 4. Weak Solution to the Flow Associated with the Minimization Problem

In this section, we present the definition of weak solution to the evolution problem (27)–(29), propose an approximating evolution equation, establish existence result for the solution of the approximate evolution equation, and, then, by logical mathematical manipulation and passing to the limits, present a proof of the existence and uniqueness of the solution of the evolution problem (27)–(29).We refer to the works in [25–28] as the motivation for our definition of the weak solution to the evolution problem (27)–(29).Definition 6.
A measurable functionu:QT→ℝ, is called an entropy solution of (27)–(29) if u∈L∞([0,T];BV(Ω)∩L2(Ω)),  ∂u/∂t∈L2(QT), u(x,0)=f in Ω in the trace sense, and if there exists z∈L1(QT;ℝN),   ∥z∥L∞(QT;ℝN)≤1, ∂u/∂t=div⁡(z)-λ(u-f) in 𝒟′(QT) such that
(31)∫Ω‍(u(t)-υ(t))∂u∂tdx≤∫Ω‍z(t)·∇vdx-∥Du(t)∥-λ∫Ω‍(u-f)(u(t)-υ)dx
for every υ∈L∞(0,T;W1,1(Ω)), a.e. for t∈[0,T].


### 4.1. Existence of the Solution of the Approximate Evolution Problem (32)–(34)

Before we study the existence and uniqueness of the evolution problem (27)–(29) above, let us consider the following approximate evolution problem: for 1<p≤2 and fp∈W1,p(Ω), we construct the approximation
(32)∂up∂t=div⁡(φp′(∇up))-λ(up-fp),(x,t)∈QT,(33)∂up∂n→=0,(x,t)∈∂QT,(34)up(x,0)=fp,x∈Ω,
where
(35)φp(∇up)=|∇up|p+11+|∇up|,
which implies that
(36)φp′(∇up)=((p+1)+p|∇up|)|∇up|p-1∇up(1+|∇up|)2.Theorem 7.
The approximate evolution equation (32)–(34) with (35) and (36) has a weak solution up∈L∞(0,T;W1,p(Ω)), ∂up/∂t∈L2(QT) such that
(37)∫‍∫QT‍[ϕ∂up∂t+φp′(∇up)·∇ϕ]dxdt=-λ∫‍∫QT‍(up-fp)ϕdxdt,(38)lim⁡t→0+∥up(x,t)-fp∥L2(Ω)=0,
where 1<p≤2, for any t∈[0,T], for every ϕ∈L2(0,T;W1,p(Ω)),  ∂ϕ/∂n→=0, a.e. on [0,T], and that ∥up∥L∞(0,T;W1,p(Ω))+∥up∥L∞(0,T;L2(Ω))+∥∂up/∂t∥L2(QT)≤C.Proof.
We apply Rothe’s method in [29] to construct an approximating solution sequence upn for the approximate evolution problem (32)–(34). We divide the interval [0,T] into n equal parts, where h=T/n. For any k:1≤k≤n, for any integer n>0 and a function u(x,t), we have
(39)upn,k(x)=up(x,kh)
for k=0,1,…,n.
We then consider the difference approximating equation of (32) as follows:
(40)upn,k-upn,k-1h=div⁡(φp′(∇upn,k))-λ(upn,k-fp).
Denoting vp:=upn,k, then the above equation becomes
(41)div⁡(φ′(∇vp))-(1h+λ)vp+(λfp+un,k-1h)=0.
The idea here is to prove that if the value of upn,k-1 is known and upn,0=fp, then (41) admits a weak solution vp:=upn,k.
From (41), we may back-project to the general functional J defined in H1(Ω), given by
(42)J(vp)=∫Ω‍(|∇vp|p+11+|∇vp|)dx+12∫Ω‍(1h+λ)vp2dx-∫Ω‍(λfp+un,k-1h)vpdx.
It can be shown that J(·) above is convex and lower semicontinuous in W1,p(Ω). Hence, there exists a minimizing sequence {upn} for J(vp) such that J(upn)=inf⁡vp∈W1,p⁡J(vp). For simplicity and without loss of generality, we will let λ=1. Then, for any test function ϕ(x)∈C0∞(Ω) with (40) and integration by parts, we have
(43)1h∫Ω‍(upn,k-upn,k-1)ϕ(x)dx+∫Ω‍φp′(∇upn,k)·∇ϕ(x)dx+∫Ω‍(upn,k-fp)ϕ(x)dx=0.
To obtain the approximate solution to the whole domain QT, we denote
(44)upn(x,t)=∑k=1n‍χn,k(t)upn,k;upn(x,0)=fp(x),vpn(x,t)=∑k=1n‍χn,k(t)wwwwwww×(upn,k-1(x)+λn,k(t)(upn,k(x)-upn,k-1(x))),
where χn,k(t) is the indicator function of t∈[(k-1)h,kh],
(45)λn,k(t)={th-(k-1),t∈[(k-1)h,kh),0,otherwise.
Equation (43) implies that
(46)∫‍∫QT‍[ϕ(x,t)∂vpn∂t+φp′(∇upn)·∇ϕ(x,t)]dxdt+∫‍∫QT‍(upn-fp)ϕ(x,t)dxdt=0
for ϕ∈C0∞(QT).
In the steps that follow, we obtain some estimates forupn(x,t) and vpn(x,t). For this purpose, let us choose the test function ϕ(x) in (43) such that ϕ(x)=upn,k-upn,k-1. We then obtain
(47)1h∫Ω‍(upn,k-upn,k-1)2dx+∫Ω‍φ′(∇upn,k)·(∇upn,k-∇upn,k-1)dx+∫Ω‍(upn,k-fp)(upn,k-upn,k-1)dx=0.
Applying convexity leads to
(48)1h∫Ω‍(upn,k-upn,k-1)2dx+12∫Ω‍(un,k-fp)2dx+∫Ω‍φp(∇upn,k)dx≤∫Ω‍φp(∇upn,k-1)dx+12∫Ω‍(upn,k-1-fp)2dx.
Summing (48) from k to r for 1≤k≤r≤n yields
(49)∫Ω‍φp(∇upn,r)dx≤∫Ω‍φp(∇fp)dx,
which implies that
(50)sup⁡0<t<T∫Ω‍φp(∇upn)dx≤∫Ω‍φp(∇fp)dx=C1(C1≡constant).
In addition, let us consider that
(51)C1≥∫Ω‍φp(∇upn)dx=∫Ω‍|∇upn|p+11+|∇upn|dx=∫|∇upn|≤k‍|∇upn|p+11+|∇upn|dx+∫|∇upn|>k‍|∇upn|p+11+|∇upn|dx≥∫|∇upn|≤1‍|∇u|p+1dx+∫|∇upn|>1‍12|∇upn|p,
for some appropriately defined constant k. Observe from the above inequality that
(52)∫|∇upn|≤1‍|∇upn|p+1dx≤C1,(53)12∫|∇upn|>1‍|∇upn|pdx≤C1.
Application of Hölder inequality yields the following inequality:
(54)∫|∇upn|≤1‍|∇upn|pdx≤(∫|∇upn|≤1‍|∇upn|p+1dx)p/(p+1)(∫|∇upn|≤1‍1dx)1/(p+1).≤C1
It is then clear from (53) and (54) that
(55)∫Ω‍|∇upn|pdx≤C1.
Hence, we can conclude from (50)–(54) that
(56)sup⁡0≤t≤T∫Ω‍|∇upn|pdx≤C1.
Now, summing (48) from k to n leads to
(57)1h∑k=1n‍∫Ω‍(upn,k-upn,k-1)2dx+12∫Ω‍(upn-fp)2dx≤∫Ω‍φp(∇fp)dx.
Equation (57) implies that
(58)1h∑k=1n‍∫Ω‍(upn,k-upn,k-1)2dx≤∫Ω‍φp(∇fp)dx=C1.
Moreover, by definition of vpn(x,t), we see that
(59)∂vpn∂t=1h∑k=1n‍χn,k(t)(upn,k-upn,k-1),
which by (58) implies that
(60)∥∂vpn∂t∥L2(QT)2=1h2∑k=1n‍h∥upn,k-upn,k-1∥2≤C1.
It can be deduced from (50) that
(61)sup⁡0<t<T∫Ω‍φp(∇vpn)dx≤C1.
Now, observe from (57) that
(62)12∫Ω‍(upn-fp)2dx≤C1.
Using Minkowski’s inequality on the strength of (62), we deduce that
(63)∥upn∥L2(Ω)-∥fp∥L2(Ω)≤∥upn-fp∥L2(Ω)≤C1,
which leads to
(64)∥upn∥L2(Ω)≤C1+∥fp∥L2(Ω)=C2.
From the above estimate, it can be deduced that
(65)sup⁡0<t<T∥vpn∥L2(Ω)≤C2.
Denoting φp′(∇upn):=Bupn, then from (56), (60), (61), (64), and (65), it can be seen that the sequences {upn},  {vpn}, {∂vpn/∂t}, and {Bupn} are bounded. Hence, there exist subsequences of {upn}, {vpn}, {∂vpn/∂t}, and {Bupn}, respectively, denoted by the same sequences such that, as n→∞,
(66)upn⇀up,weakly inL∞(0,T;L2(Ω)),upn⇀vp,weakly inL∞(0,T;L2(Ω)),upp⇀up,weakly inL∞(0,T;W1,p(Ω)),∂vpn∂t⇀∂vp∂t,weakly inL2(QT),Bupn⇀Z,weak inL∞(0,T;Lq(Ω))N,forq=pp-1,
for some up,vp, and Z.
Next, we show thatvp=up. Observe from (44)-(45) that
(67)vpn-upn=∑k=1n‍χn,k(1-λn,k)(upn,k-1-upn,k),
which by (58) leads to
(68)∥vpn-upn∥L2(QT)2≤∑k=1n‍h∥upn,k-upn,k-1∥L2(Ω)2≤C1h2.
This implies that ∥vpn-upn∥L2(QT)2≤0 as h→0. Since from the subsequences above it can be seen that upn converges to up and vpn converges to vp, both in L2(QT), it follows from (68) that vp=up.
Furthermore, from (56), (60), (64), and (65) and letting n→∞, it can be deduced that
(69)∥∂vp∂t∥L2(QT)+∥up∥L∞(0,T;L2(Ω))+∥up∥L∞(0,T;W1,p(Ω))≤C3,whereC3=(2C1+C2).
Thus, the above convergence results imply that, as n→∞ and for any ϕ∈C0∞(QT), (46) yields
(70)∫‍∫QT‍∂up∂tϕdxdt+∫‍∫QT‍Z·∇ϕdxdt+∫‍∫QT‍(up-fp)ϕdxdt=0.
Now, it remains to show thatZ=Bup. Here, we follow the works in [25, 30, 31].
Recalling thatφp′(∇upn,k):=Bupn,k, then for any w∈Lp(0,T;W1,p(Ω)) and for k from 1 to n, by the monotonicity condition, we obtain the inequality
(71)∫Ω‍(Bupn,k-Bw(t))(∇upn,k-∇w(t))dx≥0.
Now, lettingϕ=upn,k in (43), we obtain
(72)1h∫Ω‍(upn,k-upn,k-1)upn,kdx+∫Ω‍Bupn,k·∇upn,kdx+∫Ω‍(upn,k-fp)upn,kdx=0.
Applying Young’s inequality on the first term of (72) together with the inequality (71) and integrating over((k-1)h,kh), we obtain
(73)12∫Ω‍[|upn,k|2-|upn,k-1|2]dx+∫(k-1)hkh‍∫Ω‍Bupn,k·∇wdxdt+∫(k-1)hkh‍∫Ω‍Bw(∇upn,k-∇w)dxdt+h∫Ω‍(upn,k-fp)upn,kdx≤0.
Summing up (73) for k from 1 to n, we obtain
(74)12∫Ω‍[|upn(T)|2-|fp|2]dx+∫∫QT‍Bupn·∇wdxdt+∫‍∫QT‍Bw(∇upn-∇w)dxdt+∫∫QT‍(upn-fp)upndxdt≤0.
Recalling the convergence sets in (66) and letting n→∞, (74) yields
(75)12∫Ω‍[|up(T)|2-|fp|2]dx+∫‍∫QT‍Z·∇wdxdt+∫∫QT‍Bw(∇up-∇w)dxdt+∫∫QT‍(up-f)updxdt≤0.
Equation (75) may be rewritten in the form
(76)∫‍∫QT‍up∂up∂tdxdt+∫∫QT‍Z·∇wdxdt+∫∫QT‍Bw(∇up-∇w)dxdt+∫∫QT‍(up-fp)updxdt≤0.
Now, setting ϕ=up in (70), we obtain
(77)∫∫QT‍up∂up∂tdxdt+∫‍∫QT‍Z·∇updxdt+∫‍∫QT‍(up-fp)updxdt=0.
Then, substituting equation (77) into (76) leads to
(78)∫∫QT‍(Z-Bw)(∇up-∇w)dxdt≥0.
Since w is arbitrary, we may set w=up-κs where κ>0 and ∇s∈L∞(0,T;W1,p(Ω)). We then have
(79)∫∫QT‍(Z-B(up-κs))∇sdxdt≥0.
Sending κ→0, we obtain
(80)∫∫QT‍(Z-Bup)∇sdxdt≥0,∀s∈L∞(0,T;W1,p(Ω)).
Observe that, in fact, equality holds if we set s=-s in the inequality above. We thus deduce that Bup=Z. Hence, φp′(∇up)=B(up)=Z a.e. in QT. Therefore, (70) together with (80) leads to the identity (37). Thus, up is a weak solution to (32)–(34), taken together with (35) and (36).
To prove the second part, that is, the relation (38), we let in (70) ϕ=up(x,t) and then ϕ=up(x,t1), for0≥t1≥t2≥T. Consider
(81)∫Ω‍(up2(x,t2)-up2(x,t1))dx=-2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t)dxdt-2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t)dxdt,∫Ω‍(up(x,t1)up(x,t2)-up2(x,t1))dx=-∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t1)dxdt-∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t1)dxdt.Observe that from the above equation we have(82)∫Ω‍|up(x,t2)-up(x,t1)|2dx=∫Ω‍(up2(x,t2)-up2(x,t1))dx+2∫Ω‍(up2(x,t1)-up(x,t1)up(x,t2))dx=-2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t)dxdt-2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t)dxdt+2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t1)dxdt+2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t1)dxdt.
We deduce from the above equation that
(83)lim⁡t→0+∥up(x,t)-fp∥L2(Ω)=0,sinceup(x,0)=fp.


### 4.2. Existence and Uniqueness of the Solution to the Evolution Problem (27)–(29)

Theorem 8.
Letf∈BV(Ω)∩L2(Ω). Then, there exists a unique weak solution u∈L∞(0,T;BV(Ω)∩L2(Ω)), ∂u/∂t∈L2(QT), and u(x,0)=f in the trace sense.Proof.
By Theorem7, there exists up, which is a weak solution to the approximate problem (32)–(34) and a constant C such that
(84)∥up∥L∞(0,T;W1,p(Ω))+∥up∥L∞(0,T;L2(Ω))+∥∂up∂t∥L2(QT)≤C.
From (84), we deduce that there exists a subsequence of {up} denoted by {up} itself and a function u∈L∞(0,T;BV(Ω)∩L2(Ω))  with ∂u/∂t∈L2(QT) such that, as p→1+,
(85)up⟶uinL1(Ω),given∥Du∥≤liminf⁡p→1+∥Dup∥Lp(Ω),a.e.∈(0.T).
From (84), we may also deduce that
(86)∂up∂t⇀∂u∂tweakly inL2(QT).
Since BV(Ω)↪Lp(Ω), (1≤p≤N/(N-1)) continuously [7, 27, 31, 32] and also considering (84), we deduce that
(87)up⟶u,strongly inL2(QT).
Applying the method in [25], we next show that φp′(∇up) is weakly compact in L1(QT;ℝN). Using Jensen’s inequality and Hölder’s inequality, we may have that
(88)|∫‍∫QT‍φp′(∇up)dxdt|≤∫∫QT‍|φp′(∇up)|dxdt≤p∫∫QT‍|∇up|p-1dxdt≤p(∫∫QT‍|∇up|pdxdt)(p-1)/p(1)1/p≤p(C)(p-1)/pℒN(QT)1/p.
Hence, φp′(∇up) is bounded and equi-integrable in L1(QT;ℝN) and is therefore weakly compact in L1(QT;ℝN). Thus, we may deduce that
(89)φp′(∇up)⇀zweakly inL1(QT;ℝN).
Hence, we obtain
(90)∫‍∫QT‍∂u∂tϕdxdt+∫‍∫QT‍z·∇ϕdxdt+∫‍∫QT‍(u-f)ϕdxdt=0
for every ϕ∈C0∞(QT) and ∂u/∂t=div⁡(z)-λ(u-f) in 𝒟′(QT).
Now, it remains to show that∥z∥L∞(QT;ℝN)≤1.
For anyr>0 and setting Bp,r={(x,t)∈QT:|∇up|>r}, we have
(91)|Bp,r|≤Crp,for anyp>1,r>0.
Then, as above, there exist a function gr∈L1(QT:ℝN) and the indicator function χBp,r of Bp,r such that
(92)φp′(∇up)χBp,r⟶grweakly inL1(QT;ℝN)asp⟶1+.
And, for any ϕ∈L∞(QT;ℝN), it can be shown that
(93)|∫‍∫QT‍φp′(∇up)·ϕχBp,rdxdt|≤Cr.
Letting p→1+, we obtain
(94)∫‍∫QT‍|gr|dxdt≤Cr,for anyr>0.
Now, taking into consideration equations (91), (93) and (94), we deduce that
(95)|∫‍∫QT‍φp′(∇up)·ϕχQT/Bp,r|≤rp-1,for anyp>1,
letting p→1+, we get that φp′(∇up)·ϕχQT/Bp,r converges weakly to some function fr∈L1(QT;ℝN) with ∥fr∥L∞(QT;ℝN)≤1. And since for any r>0 we may write z=fr+gr, with gr satisfying condition (94), it follows that ∥z∥L∞(QT;ℝN)≤1.
Next, we verify the solution definition inequality (31). By setting ϕ=(up-υn)η(t), (37) leads to
(96)∫‍∫QT‍(up-υn)η(t)∂up∂tdxdt=-∫‍∫QT‍φp′(∇up)·∇((up-υn)η(t))dxdt-λ∫‍∫QT‍(up-fp)(up-υn)η(t)dxdt.
Taking the limit as p→1+ gives the inequality
(97)∫‍∫QT‍(u(t)-υn(t))η(t)∂u∂tdxdt≤∫‍∫QT‍z(t)·∇υnη(t)dxdt-∫0T‍∥∇u∥η(t)dt-λ∫‍∫QT‍(u-f)(u-υn)η(t)dxdt.
Taking the limit as n→∞ for every υ∈L∞(0,T;W1,1(Ω)) and recognizing the arbitrary nature of η(t), we obtain
(98)∫Ω‍(u(t)-υ(t))∂u∂tdx≤∫Ω‍z(t)·∇υdx-∥Du∥-λ∫Ω‍(u-f)(u-υ)dx,
for every υ∈L∞(0,T;W1,1(Ω)), a.e. in [0,T]. That concludes the proof of the existence of the entropy solution of problem (27)–(29).
Uniqueness of Weak Solution. Let u1 and u2 be two entropy solutions of the problem (27)–(29), such that, respectively, their initial data are u1(x,0)=f10 and u2(x,0)=f20. Then, for every υ∈L∞(0,T;W1,1(Ω)), a.e. on [0,T], there exists z1,z2∈L∞(QT;ℝN) such that
(99)∫Ω‍(u1-υ)∂u1∂tdx≤∫Ω‍z1·∇υdx-∥∇u1∥-λ∫Ω‍(u1-f)(u1-υ)dx,(100)∫Ω‍(u2-υ)∂u2∂tdx≤∫Ω‍z2·∇υdx-∥∇u2∥-λ∫Ω‍(u2-f)(u2-υ)dx.
In addition, let u1n and u2n be approximations, respectively, for u1 and u2, such that
(101)lim⁡n→∞(∥∇u1n∥L1(Ω)-∥Du1∥)=0,lim⁡n→∞(∥∇u2n∥L1(Ω)-∥Du2∥)=0,lim⁡n→∞∥u1n-u1∥=0,lim⁡n→∞∥u2n-u2∥=0,
a.e. on [0,T]. By setting υ=u2n in (99) and υ=u1n in (100), adding the two equations, applying Lemma 3, and rearranging the result, we obtain
(102)∫Ω‍(u1-u2)(∂u1∂t-∂u2∂t)dx+∫Ω‍(u1-u1n)∂u2∂tdx+∫Ω‍(u2-u2n)∂u1∂tdx≤∫Ω‍z1·∇u2ndx-∥∇u2∥+∫Ω‍z2·∇u1ndx-∥∇u1∥-λ2∫Ω‍((u1n-f)2-(u1-f)2)dx-λ2∫Ω‍((u2n-f)2-(u2-f)2)dx.
Then, integrating the inequality above from 0 to t and taking the limit as n→∞, yield
(103)∫Ω‍(u1-u2)2dx≤(f10-f20)2
This establishes the uniqueness of the entropy solution.


## 4.1. Existence of the Solution of the Approximate Evolution Problem (32)–(34)

Before we study the existence and uniqueness of the evolution problem (27)–(29) above, let us consider the following approximate evolution problem: for 1<p≤2 and fp∈W1,p(Ω), we construct the approximation
(32)∂up∂t=div⁡(φp′(∇up))-λ(up-fp),(x,t)∈QT,(33)∂up∂n→=0,(x,t)∈∂QT,(34)up(x,0)=fp,x∈Ω,
where
(35)φp(∇up)=|∇up|p+11+|∇up|,
which implies that
(36)φp′(∇up)=((p+1)+p|∇up|)|∇up|p-1∇up(1+|∇up|)2.Theorem 7.
The approximate evolution equation (32)–(34) with (35) and (36) has a weak solution up∈L∞(0,T;W1,p(Ω)), ∂up/∂t∈L2(QT) such that
(37)∫‍∫QT‍[ϕ∂up∂t+φp′(∇up)·∇ϕ]dxdt=-λ∫‍∫QT‍(up-fp)ϕdxdt,(38)lim⁡t→0+∥up(x,t)-fp∥L2(Ω)=0,
where 1<p≤2, for any t∈[0,T], for every ϕ∈L2(0,T;W1,p(Ω)),  ∂ϕ/∂n→=0, a.e. on [0,T], and that ∥up∥L∞(0,T;W1,p(Ω))+∥up∥L∞(0,T;L2(Ω))+∥∂up/∂t∥L2(QT)≤C.Proof.
We apply Rothe’s method in [29] to construct an approximating solution sequence upn for the approximate evolution problem (32)–(34). We divide the interval [0,T] into n equal parts, where h=T/n. For any k:1≤k≤n, for any integer n>0 and a function u(x,t), we have
(39)upn,k(x)=up(x,kh)
for k=0,1,…,n.
We then consider the difference approximating equation of (32) as follows:
(40)upn,k-upn,k-1h=div⁡(φp′(∇upn,k))-λ(upn,k-fp).
Denoting vp:=upn,k, then the above equation becomes
(41)div⁡(φ′(∇vp))-(1h+λ)vp+(λfp+un,k-1h)=0.
The idea here is to prove that if the value of upn,k-1 is known and upn,0=fp, then (41) admits a weak solution vp:=upn,k.
From (41), we may back-project to the general functional J defined in H1(Ω), given by
(42)J(vp)=∫Ω‍(|∇vp|p+11+|∇vp|)dx+12∫Ω‍(1h+λ)vp2dx-∫Ω‍(λfp+un,k-1h)vpdx.
It can be shown that J(·) above is convex and lower semicontinuous in W1,p(Ω). Hence, there exists a minimizing sequence {upn} for J(vp) such that J(upn)=inf⁡vp∈W1,p⁡J(vp). For simplicity and without loss of generality, we will let λ=1. Then, for any test function ϕ(x)∈C0∞(Ω) with (40) and integration by parts, we have
(43)1h∫Ω‍(upn,k-upn,k-1)ϕ(x)dx+∫Ω‍φp′(∇upn,k)·∇ϕ(x)dx+∫Ω‍(upn,k-fp)ϕ(x)dx=0.
To obtain the approximate solution to the whole domain QT, we denote
(44)upn(x,t)=∑k=1n‍χn,k(t)upn,k;upn(x,0)=fp(x),vpn(x,t)=∑k=1n‍χn,k(t)wwwwwww×(upn,k-1(x)+λn,k(t)(upn,k(x)-upn,k-1(x))),
where χn,k(t) is the indicator function of t∈[(k-1)h,kh],
(45)λn,k(t)={th-(k-1),t∈[(k-1)h,kh),0,otherwise.
Equation (43) implies that
(46)∫‍∫QT‍[ϕ(x,t)∂vpn∂t+φp′(∇upn)·∇ϕ(x,t)]dxdt+∫‍∫QT‍(upn-fp)ϕ(x,t)dxdt=0
for ϕ∈C0∞(QT).
In the steps that follow, we obtain some estimates forupn(x,t) and vpn(x,t). For this purpose, let us choose the test function ϕ(x) in (43) such that ϕ(x)=upn,k-upn,k-1. We then obtain
(47)1h∫Ω‍(upn,k-upn,k-1)2dx+∫Ω‍φ′(∇upn,k)·(∇upn,k-∇upn,k-1)dx+∫Ω‍(upn,k-fp)(upn,k-upn,k-1)dx=0.
Applying convexity leads to
(48)1h∫Ω‍(upn,k-upn,k-1)2dx+12∫Ω‍(un,k-fp)2dx+∫Ω‍φp(∇upn,k)dx≤∫Ω‍φp(∇upn,k-1)dx+12∫Ω‍(upn,k-1-fp)2dx.
Summing (48) from k to r for 1≤k≤r≤n yields
(49)∫Ω‍φp(∇upn,r)dx≤∫Ω‍φp(∇fp)dx,
which implies that
(50)sup⁡0<t<T∫Ω‍φp(∇upn)dx≤∫Ω‍φp(∇fp)dx=C1(C1≡constant).
In addition, let us consider that
(51)C1≥∫Ω‍φp(∇upn)dx=∫Ω‍|∇upn|p+11+|∇upn|dx=∫|∇upn|≤k‍|∇upn|p+11+|∇upn|dx+∫|∇upn|>k‍|∇upn|p+11+|∇upn|dx≥∫|∇upn|≤1‍|∇u|p+1dx+∫|∇upn|>1‍12|∇upn|p,
for some appropriately defined constant k. Observe from the above inequality that
(52)∫|∇upn|≤1‍|∇upn|p+1dx≤C1,(53)12∫|∇upn|>1‍|∇upn|pdx≤C1.
Application of Hölder inequality yields the following inequality:
(54)∫|∇upn|≤1‍|∇upn|pdx≤(∫|∇upn|≤1‍|∇upn|p+1dx)p/(p+1)(∫|∇upn|≤1‍1dx)1/(p+1).≤C1
It is then clear from (53) and (54) that
(55)∫Ω‍|∇upn|pdx≤C1.
Hence, we can conclude from (50)–(54) that
(56)sup⁡0≤t≤T∫Ω‍|∇upn|pdx≤C1.
Now, summing (48) from k to n leads to
(57)1h∑k=1n‍∫Ω‍(upn,k-upn,k-1)2dx+12∫Ω‍(upn-fp)2dx≤∫Ω‍φp(∇fp)dx.
Equation (57) implies that
(58)1h∑k=1n‍∫Ω‍(upn,k-upn,k-1)2dx≤∫Ω‍φp(∇fp)dx=C1.
Moreover, by definition of vpn(x,t), we see that
(59)∂vpn∂t=1h∑k=1n‍χn,k(t)(upn,k-upn,k-1),
which by (58) implies that
(60)∥∂vpn∂t∥L2(QT)2=1h2∑k=1n‍h∥upn,k-upn,k-1∥2≤C1.
It can be deduced from (50) that
(61)sup⁡0<t<T∫Ω‍φp(∇vpn)dx≤C1.
Now, observe from (57) that
(62)12∫Ω‍(upn-fp)2dx≤C1.
Using Minkowski’s inequality on the strength of (62), we deduce that
(63)∥upn∥L2(Ω)-∥fp∥L2(Ω)≤∥upn-fp∥L2(Ω)≤C1,
which leads to
(64)∥upn∥L2(Ω)≤C1+∥fp∥L2(Ω)=C2.
From the above estimate, it can be deduced that
(65)sup⁡0<t<T∥vpn∥L2(Ω)≤C2.
Denoting φp′(∇upn):=Bupn, then from (56), (60), (61), (64), and (65), it can be seen that the sequences {upn},  {vpn}, {∂vpn/∂t}, and {Bupn} are bounded. Hence, there exist subsequences of {upn}, {vpn}, {∂vpn/∂t}, and {Bupn}, respectively, denoted by the same sequences such that, as n→∞,
(66)upn⇀up,weakly inL∞(0,T;L2(Ω)),upn⇀vp,weakly inL∞(0,T;L2(Ω)),upp⇀up,weakly inL∞(0,T;W1,p(Ω)),∂vpn∂t⇀∂vp∂t,weakly inL2(QT),Bupn⇀Z,weak inL∞(0,T;Lq(Ω))N,forq=pp-1,
for some up,vp, and Z.
Next, we show thatvp=up. Observe from (44)-(45) that
(67)vpn-upn=∑k=1n‍χn,k(1-λn,k)(upn,k-1-upn,k),
which by (58) leads to
(68)∥vpn-upn∥L2(QT)2≤∑k=1n‍h∥upn,k-upn,k-1∥L2(Ω)2≤C1h2.
This implies that ∥vpn-upn∥L2(QT)2≤0 as h→0. Since from the subsequences above it can be seen that upn converges to up and vpn converges to vp, both in L2(QT), it follows from (68) that vp=up.
Furthermore, from (56), (60), (64), and (65) and letting n→∞, it can be deduced that
(69)∥∂vp∂t∥L2(QT)+∥up∥L∞(0,T;L2(Ω))+∥up∥L∞(0,T;W1,p(Ω))≤C3,whereC3=(2C1+C2).
Thus, the above convergence results imply that, as n→∞ and for any ϕ∈C0∞(QT), (46) yields
(70)∫‍∫QT‍∂up∂tϕdxdt+∫‍∫QT‍Z·∇ϕdxdt+∫‍∫QT‍(up-fp)ϕdxdt=0.
Now, it remains to show thatZ=Bup. Here, we follow the works in [25, 30, 31].
Recalling thatφp′(∇upn,k):=Bupn,k, then for any w∈Lp(0,T;W1,p(Ω)) and for k from 1 to n, by the monotonicity condition, we obtain the inequality
(71)∫Ω‍(Bupn,k-Bw(t))(∇upn,k-∇w(t))dx≥0.
Now, lettingϕ=upn,k in (43), we obtain
(72)1h∫Ω‍(upn,k-upn,k-1)upn,kdx+∫Ω‍Bupn,k·∇upn,kdx+∫Ω‍(upn,k-fp)upn,kdx=0.
Applying Young’s inequality on the first term of (72) together with the inequality (71) and integrating over((k-1)h,kh), we obtain
(73)12∫Ω‍[|upn,k|2-|upn,k-1|2]dx+∫(k-1)hkh‍∫Ω‍Bupn,k·∇wdxdt+∫(k-1)hkh‍∫Ω‍Bw(∇upn,k-∇w)dxdt+h∫Ω‍(upn,k-fp)upn,kdx≤0.
Summing up (73) for k from 1 to n, we obtain
(74)12∫Ω‍[|upn(T)|2-|fp|2]dx+∫∫QT‍Bupn·∇wdxdt+∫‍∫QT‍Bw(∇upn-∇w)dxdt+∫∫QT‍(upn-fp)upndxdt≤0.
Recalling the convergence sets in (66) and letting n→∞, (74) yields
(75)12∫Ω‍[|up(T)|2-|fp|2]dx+∫‍∫QT‍Z·∇wdxdt+∫∫QT‍Bw(∇up-∇w)dxdt+∫∫QT‍(up-f)updxdt≤0.
Equation (75) may be rewritten in the form
(76)∫‍∫QT‍up∂up∂tdxdt+∫∫QT‍Z·∇wdxdt+∫∫QT‍Bw(∇up-∇w)dxdt+∫∫QT‍(up-fp)updxdt≤0.
Now, setting ϕ=up in (70), we obtain
(77)∫∫QT‍up∂up∂tdxdt+∫‍∫QT‍Z·∇updxdt+∫‍∫QT‍(up-fp)updxdt=0.
Then, substituting equation (77) into (76) leads to
(78)∫∫QT‍(Z-Bw)(∇up-∇w)dxdt≥0.
Since w is arbitrary, we may set w=up-κs where κ>0 and ∇s∈L∞(0,T;W1,p(Ω)). We then have
(79)∫∫QT‍(Z-B(up-κs))∇sdxdt≥0.
Sending κ→0, we obtain
(80)∫∫QT‍(Z-Bup)∇sdxdt≥0,∀s∈L∞(0,T;W1,p(Ω)).
Observe that, in fact, equality holds if we set s=-s in the inequality above. We thus deduce that Bup=Z. Hence, φp′(∇up)=B(up)=Z a.e. in QT. Therefore, (70) together with (80) leads to the identity (37). Thus, up is a weak solution to (32)–(34), taken together with (35) and (36).
To prove the second part, that is, the relation (38), we let in (70) ϕ=up(x,t) and then ϕ=up(x,t1), for0≥t1≥t2≥T. Consider
(81)∫Ω‍(up2(x,t2)-up2(x,t1))dx=-2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t)dxdt-2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t)dxdt,∫Ω‍(up(x,t1)up(x,t2)-up2(x,t1))dx=-∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t1)dxdt-∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t1)dxdt.Observe that from the above equation we have(82)∫Ω‍|up(x,t2)-up(x,t1)|2dx=∫Ω‍(up2(x,t2)-up2(x,t1))dx+2∫Ω‍(up2(x,t1)-up(x,t1)up(x,t2))dx=-2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t)dxdt-2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t)dxdt+2∫t1t2‍∫Ω‍φp′(∇up(x,t))∇up(x,t1)dxdt+2∫t1t2‍∫Ω‍(up(x,t)-fp)up(x,t1)dxdt.
We deduce from the above equation that
(83)lim⁡t→0+∥up(x,t)-fp∥L2(Ω)=0,sinceup(x,0)=fp.


## 4.2. Existence and Uniqueness of the Solution to the Evolution Problem (27)–(29)

Theorem 8.
Letf∈BV(Ω)∩L2(Ω). Then, there exists a unique weak solution u∈L∞(0,T;BV(Ω)∩L2(Ω)), ∂u/∂t∈L2(QT), and u(x,0)=f in the trace sense.Proof.
By Theorem7, there exists up, which is a weak solution to the approximate problem (32)–(34) and a constant C such that
(84)∥up∥L∞(0,T;W1,p(Ω))+∥up∥L∞(0,T;L2(Ω))+∥∂up∂t∥L2(QT)≤C.
From (84), we deduce that there exists a subsequence of {up} denoted by {up} itself and a function u∈L∞(0,T;BV(Ω)∩L2(Ω))  with ∂u/∂t∈L2(QT) such that, as p→1+,
(85)up⟶uinL1(Ω),given∥Du∥≤liminf⁡p→1+∥Dup∥Lp(Ω),a.e.∈(0.T).
From (84), we may also deduce that
(86)∂up∂t⇀∂u∂tweakly inL2(QT).
Since BV(Ω)↪Lp(Ω), (1≤p≤N/(N-1)) continuously [7, 27, 31, 32] and also considering (84), we deduce that
(87)up⟶u,strongly inL2(QT).
Applying the method in [25], we next show that φp′(∇up) is weakly compact in L1(QT;ℝN). Using Jensen’s inequality and Hölder’s inequality, we may have that
(88)|∫‍∫QT‍φp′(∇up)dxdt|≤∫∫QT‍|φp′(∇up)|dxdt≤p∫∫QT‍|∇up|p-1dxdt≤p(∫∫QT‍|∇up|pdxdt)(p-1)/p(1)1/p≤p(C)(p-1)/pℒN(QT)1/p.
Hence, φp′(∇up) is bounded and equi-integrable in L1(QT;ℝN) and is therefore weakly compact in L1(QT;ℝN). Thus, we may deduce that
(89)φp′(∇up)⇀zweakly inL1(QT;ℝN).
Hence, we obtain
(90)∫‍∫QT‍∂u∂tϕdxdt+∫‍∫QT‍z·∇ϕdxdt+∫‍∫QT‍(u-f)ϕdxdt=0
for every ϕ∈C0∞(QT) and ∂u/∂t=div⁡(z)-λ(u-f) in 𝒟′(QT).
Now, it remains to show that∥z∥L∞(QT;ℝN)≤1.
For anyr>0 and setting Bp,r={(x,t)∈QT:|∇up|>r}, we have
(91)|Bp,r|≤Crp,for anyp>1,r>0.
Then, as above, there exist a function gr∈L1(QT:ℝN) and the indicator function χBp,r of Bp,r such that
(92)φp′(∇up)χBp,r⟶grweakly inL1(QT;ℝN)asp⟶1+.
And, for any ϕ∈L∞(QT;ℝN), it can be shown that
(93)|∫‍∫QT‍φp′(∇up)·ϕχBp,rdxdt|≤Cr.
Letting p→1+, we obtain
(94)∫‍∫QT‍|gr|dxdt≤Cr,for anyr>0.
Now, taking into consideration equations (91), (93) and (94), we deduce that
(95)|∫‍∫QT‍φp′(∇up)·ϕχQT/Bp,r|≤rp-1,for anyp>1,
letting p→1+, we get that φp′(∇up)·ϕχQT/Bp,r converges weakly to some function fr∈L1(QT;ℝN) with ∥fr∥L∞(QT;ℝN)≤1. And since for any r>0 we may write z=fr+gr, with gr satisfying condition (94), it follows that ∥z∥L∞(QT;ℝN)≤1.
Next, we verify the solution definition inequality (31). By setting ϕ=(up-υn)η(t), (37) leads to
(96)∫‍∫QT‍(up-υn)η(t)∂up∂tdxdt=-∫‍∫QT‍φp′(∇up)·∇((up-υn)η(t))dxdt-λ∫‍∫QT‍(up-fp)(up-υn)η(t)dxdt.
Taking the limit as p→1+ gives the inequality
(97)∫‍∫QT‍(u(t)-υn(t))η(t)∂u∂tdxdt≤∫‍∫QT‍z(t)·∇υnη(t)dxdt-∫0T‍∥∇u∥η(t)dt-λ∫‍∫QT‍(u-f)(u-υn)η(t)dxdt.
Taking the limit as n→∞ for every υ∈L∞(0,T;W1,1(Ω)) and recognizing the arbitrary nature of η(t), we obtain
(98)∫Ω‍(u(t)-υ(t))∂u∂tdx≤∫Ω‍z(t)·∇υdx-∥Du∥-λ∫Ω‍(u-f)(u-υ)dx,
for every υ∈L∞(0,T;W1,1(Ω)), a.e. in [0,T]. That concludes the proof of the existence of the entropy solution of problem (27)–(29).
Uniqueness of Weak Solution. Let u1 and u2 be two entropy solutions of the problem (27)–(29), such that, respectively, their initial data are u1(x,0)=f10 and u2(x,0)=f20. Then, for every υ∈L∞(0,T;W1,1(Ω)), a.e. on [0,T], there exists z1,z2∈L∞(QT;ℝN) such that
(99)∫Ω‍(u1-υ)∂u1∂tdx≤∫Ω‍z1·∇υdx-∥∇u1∥-λ∫Ω‍(u1-f)(u1-υ)dx,(100)∫Ω‍(u2-υ)∂u2∂tdx≤∫Ω‍z2·∇υdx-∥∇u2∥-λ∫Ω‍(u2-f)(u2-υ)dx.
In addition, let u1n and u2n be approximations, respectively, for u1 and u2, such that
(101)lim⁡n→∞(∥∇u1n∥L1(Ω)-∥Du1∥)=0,lim⁡n→∞(∥∇u2n∥L1(Ω)-∥Du2∥)=0,lim⁡n→∞∥u1n-u1∥=0,lim⁡n→∞∥u2n-u2∥=0,
a.e. on [0,T]. By setting υ=u2n in (99) and υ=u1n in (100), adding the two equations, applying Lemma 3, and rearranging the result, we obtain
(102)∫Ω‍(u1-u2)(∂u1∂t-∂u2∂t)dx+∫Ω‍(u1-u1n)∂u2∂tdx+∫Ω‍(u2-u2n)∂u1∂tdx≤∫Ω‍z1·∇u2ndx-∥∇u2∥+∫Ω‍z2·∇u1ndx-∥∇u1∥-λ2∫Ω‍((u1n-f)2-(u1-f)2)dx-λ2∫Ω‍((u2n-f)2-(u2-f)2)dx.
Then, integrating the inequality above from 0 to t and taking the limit as n→∞, yield
(103)∫Ω‍(u1-u2)2dx≤(f10-f20)2
This establishes the uniqueness of the entropy solution.


## 5. Numerical Experiments

In this section, we show the performance of the proposed formulation in denoising images involving a Gaussian white noise. The results are then compared with those obtained by the classical methods of Perona and Malik (PM) [1], Rudin et al. (TV) [8], and Guo et al. [5].The following numerical scheme has been proposed for the implementation of the model.


### 5.1. Additive Operator Splitting (AOS) Scheme

Here, we have implemented the evolution problem (27)–(29) using AOS scheme proposed by Weickert et al. in [33]. Thus, the equations are discretized as follows:
(104)λ0=0,un+1=1m∑l=1m‍[I-mτAl(uk)]-1[un+λτ(f-un)],div⁡n=(un+1-un)τ,λn=1σ2MN(u-f)div⁡n,ui,j0=fi,j=f(ih,jh),ui,0n=ui,1n,u0,jn=u1,jn,uI,in=uI-1,in,ui,Jn=ui,J-1n,
where Al(un)=[ai,j(un)],
(105)ai,j(un):={Cin+Cjn2h2,[j∈𝒩(i)],-∑N∈𝒩(i)‍Cin+CNn2h2,(j=i),0,else,Cin:=2+K|∇ui,jn|(1+|∇ui,jn|)2,
where
(106)|∇ui,jn|=12∑p,q∈𝒩(i)‍|upn-uqn|2h,
where 𝒩(i) is the set of the two neighbors of pixel i (boundary pixels have only one neighbor) and K is introduced, for convenience, as a tuning parameter in the implementation of the proposed formulation.


### 5.2. Discussion of the Results

The experiments in this work have been performed on a Compaq610 computer, having Intel(R) Core (TM)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 has been 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 is based on the maximal PSNR.At the end of iteration process, the PSNR, MAE, SSIM, and PSNRE values are recorded. PSNR and MAE values as discussed in [34] are, respectively, given by the following formulas:
(107)PSNR(u,u0)=10log⁡10MN|max⁡u0-min⁡u0|2∥u-u0∥L22dB,MAE(u,u0)=∥u-u0∥L1MN,
where u0 denotes the noise-free image, u is the denoised image, M×N is the dimension of image, and |max⁡u0-min⁡u0| yields the gray scale range of the original image.SSIM, designed by Wang et al. [35], is a quality metric used to measure the similarity between any two images. It is widely considered to work in a manner analogous to the human visual system. As opposed to the PSNR, SNR, MAE, and MSE, which are error based measures, SSIM models image distortion as a combination of loss of correlation, luminance degradation, and contrast distortion. Given any two images u and u0, SSIM measure is given by the formula
(108)SSIM(u,u0)=L(u,u0)·C(u,u0)·R(u,u0).L(u,u0)=(2μuμu0+k1)/(μu2+μu02+k1) compares the two images’ mean luminance μu and μu0. The maximal value of L(u,u0)=1, if μu=μu0, C(u,u0)=(2σuσu0+k2)/(σu2+σu02+k2), measures the closeness of contrast of the two images u and u0. Contrast is determined in terms of standard deviation, σ. Contrast comparison measure C(u,u0)=1 maximally if and only if σu=σu0, that is, when the images have equal contrast.R
(
u
,
u
0
)
=
(
σ
u
u
0
+
k
3
)
/
(
σ
u
σ
u
0
+
k
3
), where σuu0 is covariance between u and u0, is a structure comparison measure which determines the correlation between the images u and u0. 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 k1, k2, and k3 are small positive constants that avert the possibility of having zero denominators.Taking the edge mapEM(u) from (27) with K the tuning parameter, we have
(109)EM(u)=2+K|∇u|(1+|∇u|)2.
Hence, the corresponding edge similarity measure is given by
(110)PSNRE=PSNR(EM(u),EM(u0))dB[5].All these measures have been deployed not only to judge the extent of signal recovery, but also to see the extent of structural coincidence between the original image and the reconstruction image. It should be noted, however, that there is no hard and fast rule for selecting between the error based measures and the structural similarity index measure [36]. Therefore, any or both can be used at the same time to measure the quality of image recovery.We have compared results of the construction by our method to those of Perona-Malik method [1], total variation method of Rudin et al. [8], and the D-α-PM algorithm proposed by [5]. Adjustable parameters in the experiments were time step τ, thresholding\tuning parameter K, and others depending on the respective model. However, fidelity parameter λ was dynamically obtained in the course and process of iteration.In Tables1 and 2, we have, respectively, presented the results of the reconstruction for the synthetic image, Figure 1, and those of the textured image, Lena, Figure 2. The comparisons are based upon PSNR, MAE, SSIM, and PSNRE values and visual effects.Table 1
Numerical results for synthetic image(300×300) experiments.
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
SSIM
PSNRE
σ
K
τ
PM
30
4
0.25
356
1.56
36.57
3.78
0.9817
25.06
TV
30
0.25
205
1.66
34.56
4.73
0.9725
14.75
D-α-PM
30
1
0.25
40
2.33
37.2
3.18
0.9825
24.7
Our method
30
0.01
5
65
3.98
38.52
3.16
0.9905
24.90Table 2
Numerical results for Lena image(300×300) experiments.
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
SSIM
PSNRE
σ
K
τ
PM
30
5
0.25
113
0.45
25.94
12.88
0.6939
20.69
TV
30
0.25
129
1.15
27.75
10.45
0.7698
20.70
D-α-PM
30
4
2
10
0.50
28.15
7.21
0.7780
24.57
Our method
30
0.1
5
15
1.07
27.84
10.20
0.7805
26.34Synthetic image (300×300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ=30. (c) PM algorithm, K=4, τ=0.25 (356 steps). (d) TV algorithm, τ=0.25 (205 steps). (e) Our method, K=0.01,  τ=5 (65 steps). (f) D-α-PM algorithm, k=1,  τ=0.25, (40 steps).
(a)
(b)
(c)
(d)
(e)
(f)Lena image (300×300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ=30. (c) PM algorithm, K=5, τ=0.25 (113 steps). (d) TV method, τ=0.25 (129 steps). (e) Our method, K=0.1,  τ=5 (15 steps). (f) D-α-PM algorithm, K=4,  τ=2 (10 steps).
(a)
(b)
(c)
(d)
(e)
(f)For the nontextured synthetic image, Figure1, it is observed that, in spite of higher CPU time, our method enjoys superior signal and structural recovery, at PSNR=38.52, MAE=3.16, and SSIM=0.9905 (see Table 1). In addition, as evidenced by Figure 1(e), reconstructions due the proposed model do not show the relative blurriness speckles evident in Figure 1(c) (by PM method); it does not manifest the relative blockiness visible on Figure 1(d) (by TV method) neither does it show heavily jagged edges apparent in Figure 1(f) due to the D-α-PM algorithm. Comparing Figure 1(e) produced by our method and Figure 1(a)of the original image, it can be observed that the differences are extremely marginal. The edges also seem better preserved and sharper.Furthermore, although in terms of numerical metrics our method compares fairly well with those obtained by the D-α-PM, it should be noted that the increased number of parameters in the D-α that have to be manually tweaked or manipulated to obtain optimal results aggregately hinders the effectiveness of the method.For the textured image of Lena in Figure2, the metrics evidence that our method still gives better signal recovery and structural coincidence compared to its comparisons (see Table 2). By visual observation, Figure 2(c) produced by PM method shows speckles and blur, while Figure 2(d) shows heavy presence of staircase effect, and Figure 2(f) shows heavier ramp-like features. However, Figure 2(e) produced by our method not only shows no speckles and blur but also manifests significant reduction of staircase effect.In respect of the edge similarity measure PSNRE, for the nontexture image (see Figure 1), we observe a phenomenon where the PM method only marginally beats the proposed method. However, the method far outperforms the TV method. For the texture image (see Figure 2), experiments show that the proposed method generally outperforms the PM, TV, and D-α-PM methods in reconstruction of image from noise (see Tables 1 and 2).In addition, experiments demonstrate that the tuning parameterK introduced in the implementation of the proposed model is texture sensitive. In the nontextured images, K tends to be very small, while as texture increases the value of K also tends to increase. Generally observed is that the AOS implementation algorithm shows gray value invariance, stability based on extremum principle, Lyapunov functionals, and convergence even for larger values of the scale parameter τ [33]. This does not apply in the case of TV and PM, since they become unstable for τ>0.25. The D-α-PM survives this eventuality only when implemented by the AOS algorithm.Generally, therefore, the proposed formulation performs better than TV method and PM method and indeed fairly better than the D-α-PM algorithm not only in the quantitative metric measures but also in terms of visual effect. The D-α-PM method competes fairly well with the proposed method. However, the increased number of parameters in the D-α-PM may disadvantage its practical ability to produce optimal results. Moreover, our formulation seems to handle texture images better than its comparison. Finally, given the nonmonotone nature of the flux function of the D-α-PM model proposed by Guo et al. [5], it is an extremely ill-posed problem, generating both backward and forward diffusion at various times in the evolution process. Therefore, analytically speaking, its ability to arrive at steady state and produce data with minimum energy is highly in mathematical doubt. Notwithstanding the foregoing, the discrete or numerical implementation of the D-α-PM method still does give an approximated solution.


## 5.1. Additive Operator Splitting (AOS) Scheme

Here, we have implemented the evolution problem (27)–(29) using AOS scheme proposed by Weickert et al. in [33]. Thus, the equations are discretized as follows:
(104)λ0=0,un+1=1m∑l=1m‍[I-mτAl(uk)]-1[un+λτ(f-un)],div⁡n=(un+1-un)τ,λn=1σ2MN(u-f)div⁡n,ui,j0=fi,j=f(ih,jh),ui,0n=ui,1n,u0,jn=u1,jn,uI,in=uI-1,in,ui,Jn=ui,J-1n,
where Al(un)=[ai,j(un)],
(105)ai,j(un):={Cin+Cjn2h2,[j∈𝒩(i)],-∑N∈𝒩(i)‍Cin+CNn2h2,(j=i),0,else,Cin:=2+K|∇ui,jn|(1+|∇ui,jn|)2,
where
(106)|∇ui,jn|=12∑p,q∈𝒩(i)‍|upn-uqn|2h,
where 𝒩(i) is the set of the two neighbors of pixel i (boundary pixels have only one neighbor) and K is introduced, for convenience, as a tuning parameter in the implementation of the proposed formulation.


## 5.2. Discussion of the Results

The experiments in this work have been performed on a Compaq610 computer, having Intel(R) Core (TM)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 has been 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 is based on the maximal PSNR.At the end of iteration process, the PSNR, MAE, SSIM, and PSNRE values are recorded. PSNR and MAE values as discussed in [34] are, respectively, given by the following formulas:
(107)PSNR(u,u0)=10log⁡10MN|max⁡u0-min⁡u0|2∥u-u0∥L22dB,MAE(u,u0)=∥u-u0∥L1MN,
where u0 denotes the noise-free image, u is the denoised image, M×N is the dimension of image, and |max⁡u0-min⁡u0| yields the gray scale range of the original image.SSIM, designed by Wang et al. [35], is a quality metric used to measure the similarity between any two images. It is widely considered to work in a manner analogous to the human visual system. As opposed to the PSNR, SNR, MAE, and MSE, which are error based measures, SSIM models image distortion as a combination of loss of correlation, luminance degradation, and contrast distortion. Given any two images u and u0, SSIM measure is given by the formula
(108)SSIM(u,u0)=L(u,u0)·C(u,u0)·R(u,u0).L(u,u0)=(2μuμu0+k1)/(μu2+μu02+k1) compares the two images’ mean luminance μu and μu0. The maximal value of L(u,u0)=1, if μu=μu0, C(u,u0)=(2σuσu0+k2)/(σu2+σu02+k2), measures the closeness of contrast of the two images u and u0. Contrast is determined in terms of standard deviation, σ. Contrast comparison measure C(u,u0)=1 maximally if and only if σu=σu0, that is, when the images have equal contrast.R
(
u
,
u
0
)
=
(
σ
u
u
0
+
k
3
)
/
(
σ
u
σ
u
0
+
k
3
), where σuu0 is covariance between u and u0, is a structure comparison measure which determines the correlation between the images u and u0. 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 k1, k2, and k3 are small positive constants that avert the possibility of having zero denominators.Taking the edge mapEM(u) from (27) with K the tuning parameter, we have
(109)EM(u)=2+K|∇u|(1+|∇u|)2.
Hence, the corresponding edge similarity measure is given by
(110)PSNRE=PSNR(EM(u),EM(u0))dB[5].All these measures have been deployed not only to judge the extent of signal recovery, but also to see the extent of structural coincidence between the original image and the reconstruction image. It should be noted, however, that there is no hard and fast rule for selecting between the error based measures and the structural similarity index measure [36]. Therefore, any or both can be used at the same time to measure the quality of image recovery.We have compared results of the construction by our method to those of Perona-Malik method [1], total variation method of Rudin et al. [8], and the D-α-PM algorithm proposed by [5]. Adjustable parameters in the experiments were time step τ, thresholding\tuning parameter K, and others depending on the respective model. However, fidelity parameter λ was dynamically obtained in the course and process of iteration.In Tables1 and 2, we have, respectively, presented the results of the reconstruction for the synthetic image, Figure 1, and those of the textured image, Lena, Figure 2. The comparisons are based upon PSNR, MAE, SSIM, and PSNRE values and visual effects.Table 1
Numerical results for synthetic image(300×300) experiments.
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
SSIM
PSNRE
σ
K
τ
PM
30
4
0.25
356
1.56
36.57
3.78
0.9817
25.06
TV
30
0.25
205
1.66
34.56
4.73
0.9725
14.75
D-α-PM
30
1
0.25
40
2.33
37.2
3.18
0.9825
24.7
Our method
30
0.01
5
65
3.98
38.52
3.16
0.9905
24.90Table 2
Numerical results for Lena image(300×300) experiments.
Algorithm
Parameters
Number of steps
CPU time (sec)
PSNR
MAE
SSIM
PSNRE
σ
K
τ
PM
30
5
0.25
113
0.45
25.94
12.88
0.6939
20.69
TV
30
0.25
129
1.15
27.75
10.45
0.7698
20.70
D-α-PM
30
4
2
10
0.50
28.15
7.21
0.7780
24.57
Our method
30
0.1
5
15
1.07
27.84
10.20
0.7805
26.34Synthetic image (300×300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ=30. (c) PM algorithm, K=4, τ=0.25 (356 steps). (d) TV algorithm, τ=0.25 (205 steps). (e) Our method, K=0.01,  τ=5 (65 steps). (f) D-α-PM algorithm, k=1,  τ=0.25, (40 steps).
(a)
(b)
(c)
(d)
(e)
(f)Lena image (300×300). (a) Original image. (b) Noisy image corrupted by Gaussian noise for σ=30. (c) PM algorithm, K=5, τ=0.25 (113 steps). (d) TV method, τ=0.25 (129 steps). (e) Our method, K=0.1,  τ=5 (15 steps). (f) D-α-PM algorithm, K=4,  τ=2 (10 steps).
(a)
(b)
(c)
(d)
(e)
(f)For the nontextured synthetic image, Figure1, it is observed that, in spite of higher CPU time, our method enjoys superior signal and structural recovery, at PSNR=38.52, MAE=3.16, and SSIM=0.9905 (see Table 1). In addition, as evidenced by Figure 1(e), reconstructions due the proposed model do not show the relative blurriness speckles evident in Figure 1(c) (by PM method); it does not manifest the relative blockiness visible on Figure 1(d) (by TV method) neither does it show heavily jagged edges apparent in Figure 1(f) due to the D-α-PM algorithm. Comparing Figure 1(e) produced by our method and Figure 1(a)of the original image, it can be observed that the differences are extremely marginal. The edges also seem better preserved and sharper.Furthermore, although in terms of numerical metrics our method compares fairly well with those obtained by the D-α-PM, it should be noted that the increased number of parameters in the D-α that have to be manually tweaked or manipulated to obtain optimal results aggregately hinders the effectiveness of the method.For the textured image of Lena in Figure2, the metrics evidence that our method still gives better signal recovery and structural coincidence compared to its comparisons (see Table 2). By visual observation, Figure 2(c) produced by PM method shows speckles and blur, while Figure 2(d) shows heavy presence of staircase effect, and Figure 2(f) shows heavier ramp-like features. However, Figure 2(e) produced by our method not only shows no speckles and blur but also manifests significant reduction of staircase effect.In respect of the edge similarity measure PSNRE, for the nontexture image (see Figure 1), we observe a phenomenon where the PM method only marginally beats the proposed method. However, the method far outperforms the TV method. For the texture image (see Figure 2), experiments show that the proposed method generally outperforms the PM, TV, and D-α-PM methods in reconstruction of image from noise (see Tables 1 and 2).In addition, experiments demonstrate that the tuning parameterK introduced in the implementation of the proposed model is texture sensitive. In the nontextured images, K tends to be very small, while as texture increases the value of K also tends to increase. Generally observed is that the AOS implementation algorithm shows gray value invariance, stability based on extremum principle, Lyapunov functionals, and convergence even for larger values of the scale parameter τ [33]. This does not apply in the case of TV and PM, since they become unstable for τ>0.25. The D-α-PM survives this eventuality only when implemented by the AOS algorithm.Generally, therefore, the proposed formulation performs better than TV method and PM method and indeed fairly better than the D-α-PM algorithm not only in the quantitative metric measures but also in terms of visual effect. The D-α-PM method competes fairly well with the proposed method. However, the increased number of parameters in the D-α-PM may disadvantage its practical ability to produce optimal results. Moreover, our formulation seems to handle texture images better than its comparison. Finally, given the nonmonotone nature of the flux function of the D-α-PM model proposed by Guo et al. [5], it is an extremely ill-posed problem, generating both backward and forward diffusion at various times in the evolution process. Therefore, analytically speaking, its ability to arrive at steady state and produce data with minimum energy is highly in mathematical doubt. Notwithstanding the foregoing, the discrete or numerical implementation of the D-α-PM method still does give an approximated solution.


## 6. Conclusion

In this paper, we have proposed an alternative framework for total variation based image denoising models. The model is based on minimization of total variation with a functional coefficient. We set out to determine factors to consider when choosing the functional coefficient for the minimization potential, and we have established that, for a well-behaved minimization functional, it is desirable that the functional coefficient be of sublinear growth at infinity. The chosen functional coefficient, apart from its sublinear growth at infinity, must be one that results in a potential which is nondecreasing, of linear growth at infinity, and convex. Convexity can then be used to obtain lower semicontinuity, as it is also instrumental in determination of existence of a minimizer to the functional.We picked on an example to illustrate the effectiveness of the consideration. We have shown for the example that existence and uniqueness of minimizer can be established and that the evolution equation has a unique entropy solution. In addition, numerical implementation of the evolution equation on images demonstrates better denoising results than those offered by traditional competition due to PM method, TV method, and even the D-α-PM in [5].The successful application or otherwise of any regularization formulation in practical denoising always depends on the type of image under consideration, the type of noise involved, the application intended for the restored image, and indeed the platform upon which the formulation is implemented. Consequently, some of the constraints that might affect the practical performance of a regularization formulation modeled along the lines proposed here, including the sample formulation used, include the nature of the noise and even the hardware in which the implementation has been carried out. For more efficient performance, we suggest that platform with higher processing capabilities be deployed. The extent of degradation and the level of preservation required are also other factors that determine the performance of the formulation. For instance, in case of heavy degradation, it is widely recommended that a preconvolution of the noise image with Gaussian be done before the implementation of the suggested formulation to recover edges and other semantically important features. In addition, there may also be cases where noise removal may be counterproductive. This may occur when the oscillation due to noise is of the same scale as that of the features or patterns for which preservation is intended [37]. To surmount such a challenge, one may require a combination of formulations.However, generally, we have observed that models derived along the lines of the approach proposed here give formulations which not only generate better results but also submit rather progressively to mathematical analysis.



---

*Source: 939131-2014-05-05.xml*