# Existence and Uniqueness of Weak Solutions for Novel Anisotropic Nonlinear Diffusion Equations Related to Image Analysis

Authors: Anas Tiarimti Alaoui; Mostafa Jourhmane
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This paper establishes the existence and uniqueness of weak solutions for the initial-boundary value problem of anisotropic nonlinear diffusion partial differential equations related to image processing and analysis. An implicit iterative method combined with a variational approach has been applied to construct approximate solutions for this problem. Then, under some a priori estimates and a monotonicity condition, the existence of unique weak solutions for this problem has been proven. This work has been complemented by a consistent and stable approximation scheme showing its great significance as an image restoration technique.


---

## Body

## 1. Introduction

In the last three decades, nonlinear diffusion equations have inspired numerous research studies in various application ranges. Perona and Malik [1] were the first to introduce such equation in image processing and analysis in the following manner:(1)∂u∂t−∇⋅c∇u∇u=0,inΩ×0,T,c∇u∇u,n=0,on∂Ω×0,T,ux;0=u0x,inΩ,where Ω is an image domain in ℝ2 and c is a positive decreasing function defined on ℝ+.When it comes to processing a digital image, Perona and Malik chose the above model to preserve meaningful features such as edges while reducing irrelevant information such as noise in the homogeneous area. Nevertheless, this model, known as an isotropic nonlinear diffusion equation, handles an image feature with the same amount of blurring in all its directions. For instance, this process cannot successfully eliminate noises at edges [2]. Accordingly, it might be wise to consider the orientation of essential features by using anisotropic diffusion. Weickert [2] introduced this property by defining an orientation descriptor using the structure tensor, which is convenient to identify features such as corners and T-junctions. Besides, digital images present some structural difficulties; that is, they are discrete in space and image intensity values. Accordingly, it would be of great interest to adapt the diffusion to digital images’ structure by considering vertical, horizontal, and diagonal differential operators. Due to these reasons, we modeled and developed anisotropic nonlinear diffusion equations using a novel diffusion tensor.Various tools can be used to examine the existence of solutions for nonlinear partial differential equations (PDEs), such as variational techniques, monotonicity method, fixed-point theorems, iterative methods, and truncation techniques; for more detailed information, we refer to [3–7] and the references therein. These PDEs have been motivated by various applications such as image restoration and reconstruction (see, for example, [3, 4, 8–11]). Moreover, the image processing of the brain allows the localization of epileptogenic foci for the patient. A noninvasive method has been examined numerically as an inverse problem in [12].Under some challenging conditions, the existence and uniqueness of weak solutions for the Perona and Malik model have been investigated in the bounded variation spaceBVΩ [3, 13]. In some other functional frameworks, Wang and Zhou have thoroughly studied in [4] and proved the existence and uniqueness of weak solutions in the Orlicz space LlogLΩ using a new diffusion function cs=s+s+1logs+1/ss+1 for all s≥0.In this paper, we suppose thatΩ is an open-bounded domain of ℝ2 with Lipschitz boundary ∂Ω, and T is a positive number. We denote(2)∂x1u≔ux1≔∇u⋅e1,∂x2u≔ux2≔∇u⋅e2,∂x12u≔ux12≔∇u⋅e1+e2e1+e2,∂x−12u≔ux−12≔∇u⋅−e1+e2−e1+e2,where e1,e2 is the canonical basis of ℝ2. We consider the following anisotropic nonlinear parabolic initial-boundary value problem:(3)∂u∂t−∇⋅D∇u∇u=0,inΩ×0,T,D∇u∇u,n=0,on∂Ω×0,T,ux;0=u0x,inΩ,where D∇u, the diffusion tensor, is a real symmetric positive definite matrix of ℝ2×2 defined as follows:(4)D∇u=gux1+gux12+gux−122gux12−gux−122gux12−gux−122gux2+gux12+gux−122,and g:ℝ+⟶ℝ is a C1 positive decreasing function. Then, we can define ϕ:ℝ+⟶ℝ as a C2 function such that(5)ϕs=∫0srgrdr,s≥0,satisfying(6)ϕ0=ϕ′0=0,ϕs>0,ϕ′s>0,fors∈ℝ+∗,ϕ″s≥0,sϕ″s≤ϕ′s,fors∈ℝ+,0<lims⟶∞ϕsslogs<∞,0<lims⟶∞ϕ′slogs<∞,lims⟶0+ϕ′ss>0,lims⟶∞ϕ′ss=0.To construct an adaptive diffusion tensor, the functiong is approximated numerically by a cubic Hermite spline [14] that interpolates numeric data specified at 0=k0<k1<⋯<km with m∈ℕ∗:(7)gs=pkiP1,kiki+1s+vkiP2,kiki+1s+pki+1P1,ki+1kis+vki+1P2,ki+1kis,s∈ki,ki+1i∈0,1,…,m−1,pkmgkm,1s+vkmgkm,2s,s∈km,∞where p⋅ and v⋅ are the coefficients used to define the position and the velocity vector at a specific point, ki are the threshold parameters, Pj,cd is the family of the basis functions composed of polynomials of degree 3 used on the interval c,d such that(8)P1,cds=s−d22s+d−3cd−c3,P2,cds=s−d2s−cd−c2.And we may consider(9)gkm,1s=kmlogkm+22slogs+1−kmlogkms2,gkm,2s=km2logkm+2slogs+1−kmlogkm+1s2.From the definition ofϕ, we can deduce(10)ϕs=Ci+∑j=03Akiki+1,jj+2sj+2,s∈ki,ki+1,i∈0,1,…,m−1,Akm,2slogs+Akm,1logs+Cm,s∈km,∞,where Ci and Cm are constants determined by the continuity of ϕ at each ki. In this case, the values of the coefficients Akiki+1,j are determined experimentally provided that ϕ satisfying the above conditions on 0,km. Besides, we may introduce some sufficient conditions on km and Akm,. that guarantee the properties of ϕ on km,∞:(11)km≥1,Akm,2>0,Akm,1<kmAkm,2,Akm,1≥−kmlogkm2Akm,2.Anisotropic diffusion model (3) allows strong directional smoothing within the areas where ux1, ux2, ux12, or ux−12 is small and prevents blurring boundaries, contours, or corners that separate neighboring areas, where one or a combination of these differential operators has significant value.Moreover, the matrixD∇u has two eigenvalues λ+/−:(12)λ+/−=12gux1+gux2+gux12+gux−12±gux1−gux22+gux12−gux−122,with θ+/− are the corresponding eigenvectors. We can then expand the first equation of (3) into(13)∂u∂t=∇⋅λ+θ+θ+T∇u+∇⋅λ−θ−θ−T∇u.Accordingly, it is clear from the expression ofλ+/− that λ+≥λ−>0, which means that the diffusion towards θ+ is privileged over θ−. In fact, the difference λ+−λ−2=gux1−gux22+gux12−gux−122 indicates the isotropic diffusion for zero value and anisotropic diffusion for positive values.Henceforth, we will assume that the initial value satisfies(14)u0∈L2Ω,and we will introduce the following Orlicz space:(15)LlogLkmΩ=u:Ω⟶ℝ|∫Ω∩u≥kmulogudx<∞.Next, we define weak solutions for problem (3) on QT≔Ω×0,T with T>0:Definition 1.
A functionu:Q¯T⟶ℝ is a weak solution for problem (3) if the following conditions are satisfied:(i)
u∈C0,T;L2Ω∩L10,T;W1,1Ω with ∂xiu∈LlogLkmΩ for i=1,2.(ii)
For anyφ∈C1Q¯T with φ.,T=0, we have(16)−∫Ωu0xφx,0dx+∫0T∫Ω−uφt+D∇u∇u⋅∇φdxdt=0.
Now, we state our main theorem.Theorem 1.
Under assumption (14), there exists a unique weak solution for initial-boundary value problem (3).Inspired by [4], this paper will investigate the existence and uniqueness of weak solutions for problem (3) according to the following steps:(i)
First, we approximate nonlinear evolution problem (3) by nonlinear elliptic problems using an implicit iterative method (discretization in time-variable only), and then we prove the existence of a unique weak solution for each elliptic problem adopting a variational approach. These solutions constitute approximate solutions for problem (3).(ii)
Next, we show the uniqueness of solutions for initial-boundary value problem (3) using the monotonicity of the vector field D∇∇.:u∈ℝ2⟶D∇u∇u∈ℝ2.(iii)
Finally, passing to limits in some a priori energy estimates and using the monotonicity condition (17), we demonstrate the existence of weak solutions for problem (3).


## 2. Preliminaries

In this section, we state some useful lemmas that will be used later in the proofs.Lemma 1.
For alla≥0 and b≥1, we have ab≤aexpa+blogb.Proof.
Ifb≤expa, then ab≤aexpa≤aexpa+blogb.
Ifexpa<b, then a<logb, which means ab<blogb<aexpa+blogb.Lemma 2.
Supposeϕ:ℝ+⟶ℝ+ is a C2 convex function. Then, for all ξ0,ξ1∈ℝ2, we have(17)Dξ1ξ1−Dξ0ξ0⋅ξ1−ξ0≥0.Proof.
For eacht∈0,1, we put ξt=1−tξ0+tξ1. Then, we have(18)Dξ1ξ1−Dξ0ξ0=∫01dξt⋅e1gξt⋅e1e1+∫01dξt⋅e2gξt⋅e2e2+∫01dξt⋅e12gξt⋅e12e12+∫01dξt⋅e−12gξt⋅e−12e−12.
Sinceϕ″s=gs+sg′s, then we obtain(19)Dξ1ξ1−Dξ0ξ0=∫01ϕ″ξt⋅e1ξ1−ξ0⋅e1e1dt+∫01ϕ″ξt⋅e2ξ1−ξ0⋅e2e2dt+∫01ϕ″ξt⋅e12ξ1−ξ0⋅e12e12dt+∫01ϕ″ξt⋅e−12ξ1−ξ0⋅e−12e−12dt.
We conclude then(20)Dξ1ξ1−Dξ0ξ0⋅ξ1−ξ0=∫01ϕ″ξt⋅e1ξ1−ξ0⋅e12dt+∫01ϕ″ξt⋅e2ξ1−ξ0⋅e22dt+∫01ϕ″ξt⋅e12ξ1−ξ0⋅e122dt+∫01ϕ″ξt⋅e−12ξ1−ξ0⋅e−122dt≥0,which completes the proof.Lemma 3.
Uniform integrability and weak convergence [15].
AssumeΩ⊂ℝ2 is bounded, and let uii=1∞ be a sequence of functions in L1Ω satisfying(21)supiuiL1Ω<∞.
Suppose also(22)liml⟶∞supi∫Ω∩ui≥luidx=0.
Then, there exist a subsequenceuijj=1∞ and u˜∈L1Ω such that(23)uij⇀u˜,weakly inL1Ω.Lemma 4.
AssumeΩ⊂ℝ2 is bounded, and let uii=1∞ be a sequence of functions in L1Ω such that(24)supi∫Ω∩ui≥kmuiloguidx<∞.
Then, there exist a subsequenceuijj=1∞ and u˜∈L1Ω such that(25)uij⇀u˜,weakly inL1Ω,with u˜∈LlogLkmΩ.Proof.
GivenM>0, we may find an l≥km such that Ms≤slogs for all s≥l. Consequently,(26)∫Ωuidx=∫Ω∩ui<kmuidx+∫Ω∩ui≥kmuidx≤kmΩ+1M∫Ω∩ui≥kmuilogui|dx,which implies that(27)supi∫Ωuidx<∞.
On the other hand, there exists a positive constantC such that(28)∫Ω∩ui≥luidx≤1M∫Ω∩ui≥luiloguidx≤1M∫Ω∩ui≥kmuiloguidx≤CM=ɛ,which is true for all i and arbitrary ɛ>0. It follows then that(29)liml⟶∞supi∫Ω∩ui≥luidx=0.
Then, from Lemma3, there exist a subsequence uijj=1∞ of uii=1∞ and a function u˜∈L1Ω such that(30)uij⇀u˜,weakly inL1Ω.
It remains to prove thatu˜∈LlogLkmΩ.
We know that the functionfs=slogs for s≥1 is increasing and convex, and then the function fs is also convex for all s≥1. Therefore, we obtain(31)fu˜≤fuij+f′u˜u˜−uij.
Integrating the above inequality overΩN∩uij≥km with ΩN=Ω∩km≤u˜≤N, we have(32)∫ΩNfu˜dx≤∫Ω∩uij≥kmfuijdx+∫ΩN∩uij≥kmf′u˜u˜−uijdx≤∫Ω∩uij≥kmfuijdx+∫Ω∩uij≥kmf′u˜χkm≤u˜≤Nu˜−uijdx.
Sincef′u˜χkm≤u˜≤N∈L∞Ω and passing to limits as j⟶∞, we get(33)∫ΩNfu˜dx≤liminfj⟶∞∫Ω∩uij≥kmfuijdx<∞.
Then, passing to limits asN⟶∞, we deduce(34)∫Ω∩u˜≥kmu˜logu˜dx<∞.
It follows thenu˜∈LlogLkmΩ. This finishes the proof.


## 3. Approximate Solutions

In this section, we will discretize the time-variable interval0,T to get approximate solutions for problem (3). We denote h=T/N with N∈ℕ∗, and we designate by un an approximate solution at time nh. We define gradually from n=1,2,…,N the following elliptic problems:(35)un−un−1h−∇⋅D∇un∇un=0,inΩ,D∇un∇un,n=0,on∂Ω.To solve these equations step by step, we only need to prove the existence and uniqueness of weak solutions of the following elliptic problems:(36)u−u0h−∇⋅D∇u∇u=0,inΩ,D∇u∇u,n=0,on∂Ω,where h>0 and u0∈L2Ω.Definition 2.
A functionu∈L2Ω∩W1,1Ω with ∂xiu∈LlogLkmΩ for i=1,2 is called a weak solution for problem (36); if for any φ∈C1Ω¯, we have(37)∫Ωu−u0hφdx+∫ΩD∇u∇u⋅∇φdx=0.
And whenφ is a constant function, we obtain(38)∫Ωudx=∫Ωu0dx.
In order to prove the existence and uniqueness of weak solutions for problem (36), we consider the variational problem(39)minEu|u∈U,where(40)U=u∈L2Ω∩W1,1Ω|∂xiu∈LlogLkmΩwithi=1,2,∫Ωudx=∫Ωu0dx,and when u∈U, the functional E is defined as(41)Eu=∫Ωϕux1+ϕux2+ϕux12+ϕux−12+12h∫Ωu−u02dx.
It is easy to prove that (36) is the Euler–Lagrange equations of the functional E [16].Theorem 2.
Problem (36) has a unique weak solution.Proof.
Since(42)0≤infu∈UEu≤E0=12h∫Ωu02dx,then we can construct a minimizing sequence uqq=1∞ in U such that Euq<E0+1 and(43)limq⟶∞Euq=infu∈UEu.
Besides,(44)∫Ωuq2dx=∫Ωuq−u0+u02dx≤2∫Ωuq−u02+2∫Ωu02dx≤4hEuq+E0≤4h2E0+1.
It follows then(45)supquqL2Ω<∞.
On the other hand, givenɛ0>0, we may find l0=km such that(46)∫Ω∩∂xiuq≥km∂xiuqlog∂xiuqdx≤C∫Ω∩∂xiuq≥kmϕ∂xiuqdx≤CEuq<CE0+1with C=ɛ0+1/Akm,2>0 and i=1,2. It follows then that for i=1,2,(47)supq∫Ω∩∂xiuq≥km∂xiuqlog∂xiuqdx<∞.
Therefore, thanks to Lemma4 and the weak compactness of L2Ω, we can find a subsequence uqjj=1∞ of uqq=1∞ and a function u1∈L2Ω∩W1,1Ω such that(48)uqj⇀u1,weakly inL2Ω,and for i=1,2,(49)∂xiuqj⇀∂xiu1,weakly inL1Ω,∂xiu1∈LlogLkmΩ.
Therefore, we have(50)∫Ωu1dx=limj⟶∞∫Ωuqjdx=∫Ωu0dx,∫Ωu1−u02dx≤liminfj⟶∞∫Ωuqj−u02dx,and following the reasoning in the proof of Lemma 4, it is easy to show that for any a∈x1,x2,x12,x−12 and for a fixed ϵ>0, there exists l≥km such that(51)∫Ω∩∂au1≥l∂au1log∂au1dx≤liminfj⟶∞∫Ω∩∂auqj≥l∂auqjlog∂auqjdx,(52)∫Ω∩∂au1≥lϕ∂au1dx≤ɛ+Akm,2ɛ+1Akm,2liminfj⟶∞∫Ω∩∂auqj≥lϕ∂auqjdx.
Similarly, sinceϕ is increasing and convex in 0,l, then we can prove that(53)∫Ω∩∂au1<lϕ∂au1dx≤liminfj⟶∞∫Ω∩∂auqj<lϕ∂auqjdx.
Therefore, we obtain from (52) and (53) that(54)∫Ωϕ∂au1dx=∫Ω∩∂au1<lϕ∂au1dx+∫Ω∩∂au1≥lϕ∂au1dx≤ɛ+Akm,2ɛ+1Akm,2liminfj⟶∞∫Ωϕ∂auqjdx.
Thus, by lettingɛ⟶0, we get(55)∫Ωϕ∂au1dx≤liminfj⟶∞∫Ωϕ∂auqjdx,for any a∈x1,x2,x12,x−12. It follows then that(56)Eu1≤liminfj⟶∞Euqj=infu∈UEu,which signifies that u1∈U is a minimizer of the energy functional Eu, i.e.,(57)Eu1=minu∈UEu.
Furthermore, for allφ∈C1Ω¯ and t∈ℝ, we have u1+tφ−φΩ∈U with φΩ=1/Ω∫Ωφdx. Then, ρ0≤ρt where(58)ρt=Eu1+tφ−φΩ.
Hence, we haveρ′0=0, which means(59)∫Ωu1−u0hφ−φΩdx+∫ΩD∇u1∇u1⋅∇φdx=0.
Because of (50), we get(60)∫Ωu1−u0hφdx+∫ΩD∇u1∇u1⋅∇φdx=0.
We conclude then thatu1 is a weak solution for problem (36).
Now, assume that there is another weak solutionu^ of (36). Then, for every φ∈C1Ω¯, we have(61)∫Ωu^−u0hφdx+∫ΩD∇u^∇u^⋅∇φdx=0,which leads to(62)∫Ωu^−u1hφdx+∫ΩD∇u^∇u^−D∇u1∇u1⋅∇φdx=0.
Then, if we chooseφ=u^−u1 as a test function in (62), we get(63)∫Ωu^−u12hdx+∫ΩD∇u^∇u^−D∇u1∇u1⋅∇u^−∇u1dx=0.
Thanks to Lemma2, we deduce that(64)∫Ωu^−u12hdx=0.
Therefore,u^=u1 a.e. in Ω.
In conclusion, we have shown that there exists a unique weak solutionun∈U satisfying (35) for every n∈1,2,…,N. Consequently, we define an approximate solution uh for problem (3) as(65)uhx,t=u0x,t=0,u1x,t∈0,h,………,……,ujx,t∈j−1h,jh,………,……,uNx,t∈N−1h,T,for every h=T/N.


## 4. Existence and Uniqueness of Weak Solutions

Proof.
of Theorem1. In the beginning, we establish the uniqueness of solutions for problem (3). For this purpose, we suppose there exist two weak solutions u and v for problem (3). Then, we obtain the following:(66)∂u−v∂t−∇⋅D∇u∇u−D∇v∇v=0,inQT,D∇u∇u−D∇v∇v,n=0,on∂Ω×0,T,u−vx;0=0,inΩ.
By multiplying the first equation of the above problem byu−v and integrating over Ω and 0,t, we get(67)12∫Ωu−v2tdx+∫0t∫ΩD∇u∇u−D∇v∇v⋅∇u−vdxdτ=0,for every t∈0,T. Since the second term of the above equation is nonnegative (thanks to Lemma 2), it follows then u=v a.e. in QT.
Let us now find our weak solution for problem (3). We intend to send h to zero and show that a subsequence of our solutions uh of the approximate problems (35) converges to a weak solution for problem (3). To this end, we need to find some a priori estimates.
It follows from (35) that for every φ∈C1Ω¯,(68)∫Ωun−un−1hφdx+∫ΩD∇un∇un⋅∇φdx=0.
Then, by takingun as a test function in (68) and using unun−1≤un2+un−12/2, we get(69)12∫Ωun2dx+h∫ΩD∇un∇un⋅∇undx≤12∫Ωun−12dx.
For eacht∈0,T, we can find j∈1,…,N such that t∈j−1h,jh. Then, by adding all the inequalities (69) from n=1 to n=j, we get(70)12∫Ωuj2dx+h∑n=1j∫ΩD∇un∇un⋅∇undx≤12∫Ωu02dx.
Then, by definition ofuh, we obtain for t∈j−1h,jh that(71)12∫Ωuh2x,tdx+∫0jh∫ΩD∇uh∇uh⋅∇uhdxdτ≤12∫Ωu02dx.
SinceD∇uh is a symmetric positive definite matrix, we have also(72)12∫Ωuh2x,tdx+∫0t∫ΩD∇uh∇uh⋅∇uhdxdτ≤12∫Ωu02dx.
Therefore, after taking the supremum over0,T, we deduce that(73)sup0≤t≤T∫Ωuh2x,tdx+∫0T∫ΩD∇uh∇uh⋅∇uhdxdτ≤2∫Ωu02dx.
Recalling that0≤ϕs≤sϕ′s for all s≥0, then we can derive the following:(74)D∇uh∇uh⋅∇uh=∇uh⋅e1ϕ′∇uh⋅e1+∇uh⋅e2ϕ′∇uh⋅e2+∇uh⋅e12ϕ′∇uh⋅e12+∇uh⋅e−12ϕ′∇uh⋅e−12≥ϕ∇uh⋅e1+ϕ∇uh⋅e2+ϕ∇uh⋅e12+ϕ∇uh⋅e−12.
Besides, as in (46), for ∂x1uh,∂x2uh≥km, we may find a positive constant C such that(75)∂x1uhlog∂x1uh+∂x2uhlog∂x2uh≤Cϕ∇uh⋅e1+ϕ∇uh⋅e2≤CD∇uh∇uh⋅∇uh.
Thus, we conclude(76)sup0≤t≤T∫Ωuh2x,tdx<∞,∫0T∫Ω∩∂x1uh≥km∂x1uhlog∂x1uhdxdτ<∞,∫0T∫Ω∩∂x2uh≥km∂x2uhlog∂x2uhdxdτ<∞.
By Lemma4, we can find a subsequence of uh (for simplicity, we also denote it by uh) such that [17](77)uh⇀u,weakly_∗inL∞0,T;L2Ω,uh⇀u,weakly inL10,T;W1,1Ω,with(78)sup0≤t≤T∫Ωu2x,tdx<∞,∫0T∫Ω∩∂x1u≥km∂x1ulog∂x1udxdτ<∞,∫0T∫Ω∩∂x2u≥km∂x2ulog∂x2udxdτ<∞.
So, it remains to prove thatu is just a weak solution for problem (3). Let us now denote ξh=D∇uh∇uh. We will show that ξh is bounded in L2QT2, so we may find a subsequence of ξh that converges weakly in L2QT2 to a particular vector-valued function. Then, we will prove that this vector-valued function is equal almost everywhere to D∇u∇u in QT through monotonicity condition (17).
From the expression ofD∇uh, we can derive the following:(79)ξh=∂x1uh∂x1uhϕ′∂x1uhe1+∂x2uh∂x2uhϕ′∂x2uhe2+∂x12uh∂x12uhϕ′∂x12uhe12+∂x−12uh∂x−12uhϕ′∂x−12uhe−12|≤4ϕ′∂x1uh+∂x2uh.
Givenɛ1,ɛ2>0, we may find l1=l2=km such that(80)ϕ′s≤Mlogs,s≤ɛ2slogs,for all s≥km with M=ɛ1+Akm,2. Thus, we can distinguish two cases:(i)
If∂x1uh+∂x2uh<km then ξh2≤4ϕ′km2.(ii)
If∂x1uh+∂x2uh≥km then(81)ξh≤4Mlog∂x1uh+∂x2uh,ξh2≤4M21+ɛ2log2∂x1uhlog∂x1uh+∂x2uhlog∂x2uh,ξhexpξh4M≤4M1+ɛ2log2∂x1uhlog∂x1uh+∂x2uhlog∂x2uh.
Then,ξh is bounded in L2QT2, which means that we can find a subsequence of ξh (denote it also by ξh) and a function ξ∈L2QT2 such that(82)ξh⇀ξ,weakly inL2QT2.
Sinces↦sexpss≥0 is increasing and convex, then as in the proof of Lemma 4, we deduce that(83)∫0T∫Ω∩∂x1u+∂x2u≥kmξexpξ4Mdxdt≤liminfh⟶0∫0T∫Ω∩∂x1uh+∂x2uh≥kmξhexpξh4Mdxdt<∞.
Then, by using Lemma1, we get(84)∫0T∫Ωξ⋅∇udxdt≤∫0T∫Ωξ∇udxdt≤∫0T∫Ωξ∂x1u+∂x2udxdt≤km∫0T∫Ω∩∂x1u+∂x2u<kmξdxdt+∫0T∫Ω∩∂x1u+∂x2u≥kmξexpξ4Mdxdt+4M1+ɛ2∫0T∫Ω∩∂x1u≥km∂x1ulog∂x1udxdt+∫0T∫Ω∩∂x2u≥km∂x2ulog∂x2udxdt<∞.
It follows thenξ⋅∇u∈L1QT. Next, we will show that ξ=D∇u∇u a.e. in QT.
For eachφ∈C1QT with φ.,T=0, we take φx,nh as a test function in (35):(85)∫Ωunx−un−1xhφx,nhdx+∫ΩD∇un∇un⋅∇φx,nhdx=0,with n∈1,2,…,N. By summing n from 1 to N, we obtain(86)−1h∫Ωu0xφx,0dx+∑n=0N−1∫Ωunxφx,nh−φx,n+1hhdx+∑n=1N∫ΩD∇un∇un⋅∇φx,nhdx=0.
From the definition ofuh (65), we have(87)∑n=0N−1∫Ωunxφx,nh−φx,n+1hhdx=−∑n=0N−1∫nhn+1h∫Ωuhx,tφtx,thdxdt=−1h∫0T∫Ωuhx,tφtx,tdxdt.
Therefore,(88)−∫Ωu0xφx,0dx−∫0T∫Ωuhx,tφtx,tdxdt+∫0T∫ΩD∇uh∇uh⋅∇φdxdt+∑n=1N∫n−1hnh∫ΩD∇un∇un⋅∇φx,nh−∇φx,tdxdt=0.
Lettingh tend to zero, we get(89)∫Ωu0xφx,0dx+∫0T∫Ωuφtdxdt=∫0T∫Ωξ⋅∇φdxdt.
On the other hand, we letv∈L1QT with(90)∫0T∫Ω∩∂xiv≥km∂xivlog∂xivdxdt<∞,for i=1,2. We sum up inequalities (69):(91)12∫Ωuh2Tdx+∫0T∫ΩD∇uh∇uh⋅∇uhdxdt≤12∫Ωu02dx.
We have from Lemma2 that(92)∫0T∫ΩD∇uh∇uh−D∇v∇v⋅∇uh−∇vdxdt≥0.
Then, we obtain(93)12∫Ωuh2Tdx+∫0T∫ΩD∇uh∇uh⋅∇vdxdt+∫0T∫ΩD∇v∇v⋅∇uhdxdt−∫0T∫ΩD∇v∇v⋅∇vdxdt≤12∫Ωu02dx.
Lettingh⟶0 and noting that(94)∫Ωu2Tdx≤liminfh⟶0∫Ωuh2Tdx,we obtain(95)12∫Ωu2Tdx+∫0T∫Ωξ⋅∇vdxdt+∫0T∫ΩD∇v∇v⋅∇udxdt−∫0T∫ΩD∇v∇v⋅∇vdxdt≤12∫Ωu02dx.
By usingφ=u in (89), we get(96)12∫Ωu2Tdx+12∫Ωu02dx=∫0T∫Ωξ⋅∇udxdt.
Combining (95) with (96), we have(97)∫0T∫Ωξ−D∇v∇v⋅∇v−∇udxdt≤−∫Ωu2Tdx.
Now, settingv=u+λw for any λ>0, w∈W1,2QT, we derive from the above inequality that(98)∫0T∫Ωξ−D∇u+λw∇u+λw⋅∇wdxdt≤0.
By lettingλ⟶0 and using Lebesgue’s dominated convergence theorem, we obtain(99)∫0T∫Ωξ−D∇u∇u⋅ψdxdt=0,for every ψ∈L2Ω2. It follows then(100)ξ=D∇u∇u,a.e. inQT.
Therefore, we conclude from (89) that(101)−∫Ωu0xφx,0dx+∫0T∫Ω−uφt+D∇u∇u⋅∇φdxdt=0,for any φ∈C1Q¯T with φ.,T=0. Finally, we need to prove that u∈C0,T,L2Ω. If we choose φ∈C0∞QT in (89), we obtain(102)∫0T∫Ωuφtdxdt=∫0T∫Ωξ⋅∇φdxdt.
Sinceξ∈L2QT2, we conclude that ut∈L10,T;H−1Ω where H−1Ω is the dual space of W01,2Ω. Since(103)u=∫0tuτdτ+u0,u0∈L2Ω↪H−1Ω.
It follows then thatu∈C0,T;H−1Ω. Besides, for every h>0, let vhx,t=ux,t+h be the weak solution for problem (3) satisfying vhx;0=ux,h. Then, whx,t=ux,t+h−ux,t satisfies(104)∂wh∂t−∇⋅D∇vh∇vh−D∇u∇u=0,inΩ×0,T,D∇vh∇vh−D∇u∇u,n=0,on∂Ω×0,T,whx;0=ux,h−u0x,inΩ.
For eacht0∈0,T, we may choose wh as a test function in the first equation for problem (104) over 0,t0:(105)12∫Ωwh2x,t0dx+∫0t0∫ΩD∇vh∇vh−D∇u∇u⋅∇vh−∇udxdt≤12∫Ωwh2x,0dx.
Because of Lemma2, we deduce(106)∫Ωux,t0+h−ux,t02dx≤∫Ωux,h−u0x2dx.
Now, in order to prove thatu∈C0,T,L2Ω, we need to prove(107)limsuph⟶0+∫Ωux,h−u0x2dx=0.
We suppose that (107) is not true. Then, there exist a positive number δ and a sequence hi with hi⟶0 as i⟶∞ such that(108)limhi⟶0+∫Ωux,hi−u0x2dx≥δ.
From estimate (72), we have(109)∫Ωux,hi2dx≤∫Ωu0x2dx.
Then, from (108), we get(110)liminfhi⟶0+∫Ωu0x2dx−∫Ωu0xux,hidx≥δ2.
From (109), we conclude that ux,hi is a bounded sequence in L2Ω. Then, we may find a subsequence (denote it also by ux,hi) such that there exists a function u˜0∈L2Ω such that(111)ux,hi⇀u˜0,weakly inL2Ω.
Sinceu∈C0,T;H−1Ω, it follows that(112)ux,hi⇀u0,weakly inH−1Ω.
Therefore, we must haveu˜0=u0, and since u∈C0,T;H−1Ω, it follows that(113)ux,hi⇀u0,weakly inL2Ω,which is contradictory with (110). Therefore, we conclude that (107) is true and u∈C0,T,L2Ω. This completes the proof of Theorem 1.


## 5. Numerical Implementation and Experimental Results

### 5.1. Consistent and Stable Symmetric Finite Difference Approximation

In this section, we provide a consistent and stable discretization scheme using symmetric finite difference approximation: at timetn=nδt, n≥0, and the mesh points xi=iδ, yj=jδ0≤i≤N+1and0≤j≤M+1, and we denote by ui,jn the finite difference approximation of uxi,yj;tn. The time-space derivatives are discretized in the following manner:(114)ux1xi,yj;tn=uxi+1/2,yj;tn−uxi−1/2,yj;tnδ+Oδ2,ux2xi,yj;tn=uxi,yj+1/2;tn−uxi,yj−1/2;tnδ+Oδ2,ux12xi,yj;tn=uxi+1/2,yj+1/2;tn−uxi−1/2,yj−1/2;tn2δ+Oδ2,ux−12xi,yj;tn=uxi−1/2,yj+1/2;tn−uxi+1/2,yj−1/2;tn2δ+Oδ2,utxi,yj;tn=uxi,yj;tn+1−uxi,yj;tnδt+Oδt.By assumeδ=1 and denote(115)gNi,jn=gΔNui,jn,gEi,jn=gΔEui,jn,gSi,jn=gΔSui,jn,gWi,jn=gΔWui,jn,gNEi,jn=gΔNEui,jn2,gSEi,jn=gΔSEui,jn2,gSWi,jn=gΔSWui,jn2,gNWi,jn=gΔNWui,jn2,,withΔNui,jn=ui,j+1n−ui,jn,ΔEui,jn=ui+1,jn−ui,jn,ΔSui,jn=ui,j−1n−ui,jn,ΔWui,jn=ui−1,jn−ui,jn,ΔNEui,jn=ui+1,j+1n−ui,jn,ΔSEui,jn=ui+1,j−1n−ui,jn,ΔSWui,jn=ui−1,j−1n−ui,jn,ΔNWui,jn=ui−1,j+1n−ui,jn.Then, we may approximate problem (3) using the above scheme to obtain the following nonlinear diffusion filter:(116)ui,jn+1=ui,jn+δtgNΔNu+gEΔEu+gSΔSu+gWΔWu+gNEΔNEu+gSEΔSEu+gSWΔSWu+gNWΔNWu2i,jn,for 1≤i≤N, 1≤j≤1≤i≤M, and n≥0, with the initial condition ui,j0 and the discrete Neumann boundary condition:(117)u0,jn=u1,jn,uN+1,jn=uN,jn,for1≤j≤M,ui,0n=ui,1n,ui,M+1n=ui,Mn,for1≤i≤N,u0,0n=u1,1n,uN+1,0n=uN,1n,u0,M+1n=u1,Mn,uN+1,M+1n=uN,Mn.A unique sequenceunn∈ℕ is produced when using filter (116) on a particular initial image u0 [2]. Besides, due to the continuity of the function g, the sequence un depends continuously on u0 for every finite n. Furthermore, equation (116) satisfies the following maximum-minimum principle, which describes a stability condition for the discrete scheme.Theorem 3.
Discrete extremum principle [1, 2].
For an iteration step satisfying(118)0<δt<16g0,scheme (116) satisfies(119)mini,jui,j0≤ui,jn≤maxi,jui,j0,for all 1≤i≤N, 1≤j≤M, and n∈ℕ.


### 5.2. Experimental Results

This section will show the performance of proposed diffusion filter (116) in the image denoising process, under the boundary and initial conditions (117) while respecting the requirements concerning ϕ (Section 1), and δt (118). We will use the Peak Signal-to-Noise Ratio (PSNR that is a positive value) [18] and the Structural SIMilarity Index (SSIM that lies in 0,1) [19] to evaluate the quality of the restored images. The best results for the denoising process are equivalent to the higher value of these metrics.For comparative purposes, we will examine the proposed diffusion function with another one that has the same properties using the same filter (116). Therefore, we will use the following diffusion functions:(i)
The proposed diffusion function (m=1 for instance):(120)gs=p0P1,0ks+v0P2,0ks+pkP1,k0s+vkP2,k0s,s∈0,k,pkgk,1s+vkgk,2s,s∈k,∞.(ii)
The Wang and Zhou diffusion function (WZ) [4]:(121)gs=1s+1+logs+1s.Additionally, we will consider real test images Figure1 and evaluate our model’s performance on these images, which will be corrupted with different levels of Gaussian white noises with zero mean and variance σ2.Figure 1
Brain MRI scans: Patient30: sagittal T1 of a 30-year-old female patient [20]. Patient50: coronal T2 of a 50-year-old female patient [21]. Patient55: axial T2 of a 55-year-old patient [22]. (a) Patient30, (b) Patient50, and (c) Patient55.
(a)(b)(c)Table1 shows the quantitative results on real images, corrupted with various Gaussian noises, filtered by discrete model (116) using proposed diffusion function (120) and the one proposed by WZ (121). These results are obtained using the optimal parameters determined experimentally, as in Table 2 for each diffusion function.Table 1
PSNR and SSIM values of the images in Figure 1 affected by different values of Gaussian noise σ2 and their corresponding iteration number for both functions.
NoisyWZ [4]Proposedσ2PSNRSSIMPSNRSSIMIterPSNRSSIMIterPatient300.00523.43110.342633.37570.91143033.83330.9372130.01020.68390.243231.33670.87664131.57370.9145210.01519.11900.196030.15650.85214830.20080.8980240.02017.97630.166029.09050.82985629.00070.8825260.10012.07560.059622.68210.653511721.89530.738956Patient500.00523.66860.442431.11850.87082431.29340.8769110.01020.75610.327829.08930.81863529.15270.8258190.01519.13130.270727.91370.78204127.89810.7920260.02017.98380.233427.00110.75084726.97140.7625270.10011.99850.091121.44690.55859320.96740.593262Patient550.00524.01790.386731.31900.90212631.46680.9258190.01021.23030.289229.13100.86003629.09970.8887280.01519.59900.240327.66400.82214427.56350.8561350.02018.42920.209626.77170.79384926.54680.8305420.10012.22340.088220.69480.56429820.17550.623799Table 2
The best possible parameters for different diffusion functions.
WZ [4]Proposedσ2δtδtkp0pkv0vkPatient300.0050.083310.147014.614111.131910.66151−0.00011−0.104410.0100.083310.149915.001911.108910.45921−0.00011−0.043510.0150.083310.150515.202211.106010.45671−0.00021−0.038510.0200.083310.147015.864111.132810.46991−0.00021−0.039610.1000.083310.150015.080811.101110.45791−0.00091−0.03951Patient500.0050.083310.149413.180111.092410.56271−0.00021−0.011510.0100.083210.162311.899311.001710.56991−0.00001−0.002010.0150.083310.165912.069510.990110.56981−0.00011−0.114710.0200.083310.161011.798911.000010.56981−0.00011−0.010510.1000.083310.153013.870010.987010.55511−0.00031−0.13411Patient550.0050.083210.146813.500811.128910.58971−0.00021−0.165710.0100.083310.147013.314311.126210.58281−0.00131−0.173110.0150.083310.151413.099911.090210.58881−0.01401−0.189910.0200.083310.148813.396911.087210.53891−0.00011−0.158010.1000.083310.149912.700911.089010.56861−0.00021−0.20951It can be seen from Table1 and Figure 2 that the proposed model shows remarkable results against the WZ model. From a visual comparison, Figure 2 shows that the restored images using the proposed diffusion function have considerable noise removal and preserve the image essential features better than the restored images by the WZ diffusion function. Besides, compared with the WZ diffusion function, the results from Table 1 prove that the suggested approach has higher values in SSIM, whereas the WZ model shows significant results in PSNR while σ2-value increases.Figure 2
Visual comparison on different real images corrupted by Gaussian noise withσ2=0.1 and restored images using both functions. (a) Original-Patient30. (b) Original-Patient50. (c) Original-Patient55. (d) Noisy-Patient30. (e) Noisy-Patient50. (f) Noisy-Patient50. (g) WZ-Patient30. (h) WZ-Patient50. (i) WZ-Patient55. (j) Proposed-Patient30. (k) Proposed-Patient50. (l) Proposed-Patient55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 5.1. Consistent and Stable Symmetric Finite Difference Approximation

In this section, we provide a consistent and stable discretization scheme using symmetric finite difference approximation: at timetn=nδt, n≥0, and the mesh points xi=iδ, yj=jδ0≤i≤N+1and0≤j≤M+1, and we denote by ui,jn the finite difference approximation of uxi,yj;tn. The time-space derivatives are discretized in the following manner:(114)ux1xi,yj;tn=uxi+1/2,yj;tn−uxi−1/2,yj;tnδ+Oδ2,ux2xi,yj;tn=uxi,yj+1/2;tn−uxi,yj−1/2;tnδ+Oδ2,ux12xi,yj;tn=uxi+1/2,yj+1/2;tn−uxi−1/2,yj−1/2;tn2δ+Oδ2,ux−12xi,yj;tn=uxi−1/2,yj+1/2;tn−uxi+1/2,yj−1/2;tn2δ+Oδ2,utxi,yj;tn=uxi,yj;tn+1−uxi,yj;tnδt+Oδt.By assumeδ=1 and denote(115)gNi,jn=gΔNui,jn,gEi,jn=gΔEui,jn,gSi,jn=gΔSui,jn,gWi,jn=gΔWui,jn,gNEi,jn=gΔNEui,jn2,gSEi,jn=gΔSEui,jn2,gSWi,jn=gΔSWui,jn2,gNWi,jn=gΔNWui,jn2,,withΔNui,jn=ui,j+1n−ui,jn,ΔEui,jn=ui+1,jn−ui,jn,ΔSui,jn=ui,j−1n−ui,jn,ΔWui,jn=ui−1,jn−ui,jn,ΔNEui,jn=ui+1,j+1n−ui,jn,ΔSEui,jn=ui+1,j−1n−ui,jn,ΔSWui,jn=ui−1,j−1n−ui,jn,ΔNWui,jn=ui−1,j+1n−ui,jn.Then, we may approximate problem (3) using the above scheme to obtain the following nonlinear diffusion filter:(116)ui,jn+1=ui,jn+δtgNΔNu+gEΔEu+gSΔSu+gWΔWu+gNEΔNEu+gSEΔSEu+gSWΔSWu+gNWΔNWu2i,jn,for 1≤i≤N, 1≤j≤1≤i≤M, and n≥0, with the initial condition ui,j0 and the discrete Neumann boundary condition:(117)u0,jn=u1,jn,uN+1,jn=uN,jn,for1≤j≤M,ui,0n=ui,1n,ui,M+1n=ui,Mn,for1≤i≤N,u0,0n=u1,1n,uN+1,0n=uN,1n,u0,M+1n=u1,Mn,uN+1,M+1n=uN,Mn.A unique sequenceunn∈ℕ is produced when using filter (116) on a particular initial image u0 [2]. Besides, due to the continuity of the function g, the sequence un depends continuously on u0 for every finite n. Furthermore, equation (116) satisfies the following maximum-minimum principle, which describes a stability condition for the discrete scheme.Theorem 3.
Discrete extremum principle [1, 2].
For an iteration step satisfying(118)0<δt<16g0,scheme (116) satisfies(119)mini,jui,j0≤ui,jn≤maxi,jui,j0,for all 1≤i≤N, 1≤j≤M, and n∈ℕ.


## 5.2. Experimental Results

This section will show the performance of proposed diffusion filter (116) in the image denoising process, under the boundary and initial conditions (117) while respecting the requirements concerning ϕ (Section 1), and δt (118). We will use the Peak Signal-to-Noise Ratio (PSNR that is a positive value) [18] and the Structural SIMilarity Index (SSIM that lies in 0,1) [19] to evaluate the quality of the restored images. The best results for the denoising process are equivalent to the higher value of these metrics.For comparative purposes, we will examine the proposed diffusion function with another one that has the same properties using the same filter (116). Therefore, we will use the following diffusion functions:(i)
The proposed diffusion function (m=1 for instance):(120)gs=p0P1,0ks+v0P2,0ks+pkP1,k0s+vkP2,k0s,s∈0,k,pkgk,1s+vkgk,2s,s∈k,∞.(ii)
The Wang and Zhou diffusion function (WZ) [4]:(121)gs=1s+1+logs+1s.Additionally, we will consider real test images Figure1 and evaluate our model’s performance on these images, which will be corrupted with different levels of Gaussian white noises with zero mean and variance σ2.Figure 1
Brain MRI scans: Patient30: sagittal T1 of a 30-year-old female patient [20]. Patient50: coronal T2 of a 50-year-old female patient [21]. Patient55: axial T2 of a 55-year-old patient [22]. (a) Patient30, (b) Patient50, and (c) Patient55.
(a)(b)(c)Table1 shows the quantitative results on real images, corrupted with various Gaussian noises, filtered by discrete model (116) using proposed diffusion function (120) and the one proposed by WZ (121). These results are obtained using the optimal parameters determined experimentally, as in Table 2 for each diffusion function.Table 1
PSNR and SSIM values of the images in Figure 1 affected by different values of Gaussian noise σ2 and their corresponding iteration number for both functions.
NoisyWZ [4]Proposedσ2PSNRSSIMPSNRSSIMIterPSNRSSIMIterPatient300.00523.43110.342633.37570.91143033.83330.9372130.01020.68390.243231.33670.87664131.57370.9145210.01519.11900.196030.15650.85214830.20080.8980240.02017.97630.166029.09050.82985629.00070.8825260.10012.07560.059622.68210.653511721.89530.738956Patient500.00523.66860.442431.11850.87082431.29340.8769110.01020.75610.327829.08930.81863529.15270.8258190.01519.13130.270727.91370.78204127.89810.7920260.02017.98380.233427.00110.75084726.97140.7625270.10011.99850.091121.44690.55859320.96740.593262Patient550.00524.01790.386731.31900.90212631.46680.9258190.01021.23030.289229.13100.86003629.09970.8887280.01519.59900.240327.66400.82214427.56350.8561350.02018.42920.209626.77170.79384926.54680.8305420.10012.22340.088220.69480.56429820.17550.623799Table 2
The best possible parameters for different diffusion functions.
WZ [4]Proposedσ2δtδtkp0pkv0vkPatient300.0050.083310.147014.614111.131910.66151−0.00011−0.104410.0100.083310.149915.001911.108910.45921−0.00011−0.043510.0150.083310.150515.202211.106010.45671−0.00021−0.038510.0200.083310.147015.864111.132810.46991−0.00021−0.039610.1000.083310.150015.080811.101110.45791−0.00091−0.03951Patient500.0050.083310.149413.180111.092410.56271−0.00021−0.011510.0100.083210.162311.899311.001710.56991−0.00001−0.002010.0150.083310.165912.069510.990110.56981−0.00011−0.114710.0200.083310.161011.798911.000010.56981−0.00011−0.010510.1000.083310.153013.870010.987010.55511−0.00031−0.13411Patient550.0050.083210.146813.500811.128910.58971−0.00021−0.165710.0100.083310.147013.314311.126210.58281−0.00131−0.173110.0150.083310.151413.099911.090210.58881−0.01401−0.189910.0200.083310.148813.396911.087210.53891−0.00011−0.158010.1000.083310.149912.700911.089010.56861−0.00021−0.20951It can be seen from Table1 and Figure 2 that the proposed model shows remarkable results against the WZ model. From a visual comparison, Figure 2 shows that the restored images using the proposed diffusion function have considerable noise removal and preserve the image essential features better than the restored images by the WZ diffusion function. Besides, compared with the WZ diffusion function, the results from Table 1 prove that the suggested approach has higher values in SSIM, whereas the WZ model shows significant results in PSNR while σ2-value increases.Figure 2
Visual comparison on different real images corrupted by Gaussian noise withσ2=0.1 and restored images using both functions. (a) Original-Patient30. (b) Original-Patient50. (c) Original-Patient55. (d) Noisy-Patient30. (e) Noisy-Patient50. (f) Noisy-Patient50. (g) WZ-Patient30. (h) WZ-Patient50. (i) WZ-Patient55. (j) Proposed-Patient30. (k) Proposed-Patient50. (l) Proposed-Patient55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 6. Conclusion

This paper principally investigates the class of anisotropic diffusion partial differential equations related to image processing and analysis. The existence and uniqueness of weak solutions for this problem have been proven under sufficient conditions satisfied byϕ. A consistent and stable numerical approximation has been applied, and a discrete nonlinear filter has been tested and revealed its efficiency in the image restoration field.



---

*Source: 5553126-2021-03-08.xml*