# Existence and Uniqueness of Weak Solutions for a New Class of Convex Optimization Problems Related to Image Analysis

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

---

## Abstract

This paper proposes a new anisotropic diffusion model in image restoration that is understood from a variational optimization of an energy functional. Initially, a family of new diffusion functions based on cubic Hermite spline is provided for optimal image denoising. After that, the existence and uniqueness of weak solutions for the corresponding Euler–Lagrange equation are proven in an appropriate functional space, and a consistent and stable numerical model is also shown. We complement this work by illustrating some experiments on different actual brain Magnetic Resonance Imaging (MRI) scans, showing the proposed model’s impressive results.


---

## Body

## 1. Introduction

During the last few decades, the digital image has been significantly utilized as a noninvasive medical diagnosis or treatment tool. For instance, a different method has been developed in the electroencephalogram (EEG) (see [1]). Since a digital image contains relevant information that should be carefully extracted from it, many engineers and mathematicians have established numerous techniques and theories [2–4] to collect quantitative or qualitative data. At that point, by analyzing these reliable data and understanding its contents, one can make the right decisions at the diagnosis time.Along with the meaningful details, digital images carry false information as noise and other undesirable artifacts due to various reasons. This erroneous information can occur during image formation, transmission, or recording processes. Therefore, one needs good models to identify and eliminate these noises while preserving relevant information and structures.Image denoising is a necessary preprocessing step for many applications that rely on image quality, such as image segmentation and pattern recognition. Hence, how to restore a degraded image to its original form? Many approaches have been developed to process images and showed great interest in using the variational approach [5–15], which attempts to establish mathematical models strictly related to the physical world by defining a diffusion equation from a particular optimization problem, which generally has the following form:(1)minuFu≔∫Ωφ∇u+ρ2u−u02dΩ,where Ω⊂ℝ2 is an open bounded domain, u0 and u are the observed and the reconstructed images defined as functions of Ω⊂ℝ2⟶ℝ that associate each pixel x∈Ω to the gray level ux or u0x, φ:ℝ+⟶ℝ is a nonnegative increasing function with φ0=0, and ρ∈ℝ+ is a weighting parameter that enables to adjust the influence of the data term in the regularizing term.One common approach to solve the problem (1) is by seeking the steady-state solution of the following PDE, which corresponds to the Euler–Lagrange equation of the energy functional Fu:(2)∂u∂t=∇⋅φ′∇u∇u∇u−ρu−u0.Perona and Malik (Perona–Malik) [16] were the first to propose such a model in image processing, which corresponds to nonconvex energy functional. Unfortunately, this property leads to multiple solutions with staircase effects observed in practice [6, 17]. In [8], Charbonnier et al. suggested a strictly convex energy functional to circumvent the Perona–Malik model’s ill-posedness. Besides, many variational approaches have also been developed during the last thirty years. The most famous one was the total variation denoising suggested by Rudin et al. in their seminal paper [5]. Several improved nonlinear diffusion models for image denoising derived from Perona–Malik or TV models have been proposed in the last twenty years; for more detailed information, we refer to [2, 10–15, 18] and the references therein.Nevertheless, because of the above model’s isotropy, the diffusion equation (2) is controlled by the function s⟶φ′s/s and produces the same amount of blurring in all its directions. This means that the process cannot successfully eliminate noises at the edges. Thus, it would be important to consider smoothing along edges by adopting anisotropic diffusion [19].Furthermore, the type (1) model uses the gradient magnitude as a local descriptor operator for the image edge detector. However, digital images present some difficulties in their discrete structure; they are discrete in space and discrete in intensity value. Consequently, one may need to adapt to the digital image structure by considering differential operators that respond to vertical, horizontal, and diagonal edges and using a consistent approximation. Therefore, motivated by the above reasoning and inspired by Weickert’s anisotropic idea [19], we conceived new convex optimization problems that lead to anisotropic diffusion equations with a novel matrix diffusion tensor.Therefore, this paper provides a new variational PDE model based on directional edge detectors in Section2. The construction of a new diffusion function using cubic Hermite spline and the existence and uniqueness of weak solutions are illustrated in Section 3. Next, Section 4 implements a consistent explicit symmetric finite difference approximation that involves the 3×3 neighborhoods at each point and shows the experiments on different real images. We terminate this work with a conclusion in Section 5.


## 2. Novel Anisotropic Diffusion Model for Image Restoration

Let us now consider the rectangular domainΩ≔0,a×0,b⊂ℝ2 with Γ≔∂Ω being its boundary. This approach’s general idea is to search for a solution among the minima of some functional energy. The proposed functional energy has a regularizing term that depends on the combination of a function φ and local structures of the image expressed by the four directional derivatives:(3)∂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. Hence, our problem consists of solving the following energy minimization problem:(4)minu∈UEu,such that(5)Eu=∫Ωφux1+φux2+φux12+φux−12+ρ2u−u02dx=∫ΩLx,u,ux1,ux2dx,under the following assumptions:(6)H1φ:ℝ+⟶ℝisaC2−function/φ0=φ′0=0,φs>0,φ′s>0,s>0,φ″s>0,s≥0,lims⟶+∞φss=+∞,ρ∈ℝ+∗.The Euler–Lagrange equation ofEu (5) is obtained from(7)∂L∂u−∂x1∂L∂ux1−∂x2∂L∂ux2=0,which leads to(8)ρu−u0−∂x1φ′ux1ux1ux1+φ′ux122ux12ux12−φ′ux−122ux−12ux−12+∂x2φ′ux2ux2ux2+φ′ux122ux12ux12+φ′ux−122ux−12ux−12=0.This is equivalent to(9)ρu−u0−∂x1φ′ux1ux1ux1+∂x2φ′ux2ux2ux2+∂x12φ′ux12ux12ux12+∂x−12φ′ux−12ux−12ux−12=0.Letg:ℝ+⟶ℝ be a positive function such that gs=φ′s/s. We may reformulate equation (9) and present it as follows.Proposition 1.
The Euler–Lagrange equation (9) is equivalent to(10)ρu−u0−∇⋅D∇u∇u=0,where D∇u is a real symmetric positive definite matrix of ℝ2×2 defined as follows:(11)D∇u=gux1+gux12+gux−122gux12−gux−122gux12−gux−122gux2+gux12+gux−122.Proof.
First, it is clear that(12)∂x1gux1ux1+∂x2gux2ux2=∇⋅gux100gux2ux1ux2.
On the other hand, we have(13)∂x12gux12ux12+∂x−12gux−12ux−12=∇gux12∇u⋅e1+e2e1+e2⋅e1+e2e1+e2+∇gux−12∇u⋅−e1+e2−e1+e2⋅−e1+e2−e1+e2=12∂x1gux12ux1+ux2+∂x2gux12ux1+ux2−∂x1gux−12−ux1+ux2+∂x2gux−12−ux1+ux2=12∂x1gux12+gux−12ux1+∂x2gux12+gux−12ux2+∂x2gux12−gux−12ux1+∂x1gux12−gux−12ux2=12∇⋅gux12+gux−12ux1ux2+gux12−gux−12ux2ux1=12∇⋅gux12+gux−12gux12−gux−12gux12−gux−12gux12+gux−12ux1ux2.
Besides, the matrixD∇u has two eigenvalues λ± as follows:(14)λ±=12gux1+gux2+gux12+gux−12±gux1−gux22+gux12−gux−122.
Sinceg>0 and(15)λ−=2gux1gux2+gux12gux−12+gux1+gux2gux12+gux−122λ+,
then we deduce that the matrixD∇u is symmetric positive definite, which completes the proof.
By applying the steepest descent method, the Euler–Lagrange equation (10) can be solved considering the steady state of the following evolution equation:(16)∂u∂t=∇⋅D∇u∇u−ρu−u0,inΩ×0,T,D∇u∇,n=0,on∂Ω×0,T,ux;0=u0x,inΩ.In [20], the existence and uniqueness of weak solutions for problem (16) were proven if the weighting parameter ρ=0.


## 3. Mathematical Investigation of the Proposed Model

### 3.1. Description of the Proposed Model

Our model’s main objective is to allow strong directional smoothing within the areas whereux1, ux2, ux12, or ux−12 is small and prevent blurring boundaries, contours, or corners that separate neighboring areas, where one or a combination of these differential operators has substantial value. Some sufficient assumptions have to be made to ensure that our nonlinear diffusion meets acceptable behavior along these lines. Therefore, we can assume, as exhibited by Perona–Malik [16], some minimal hypotheses on g:(17)H2:g:ℝ+⟶ℝisC1−function such thatgs>0,g′s≤0,lims⟶0gs>0,lims⟶∞gs=0.Besides, due to the decreasing positivity of the functiong, it is evident that our model’s behavior (16) encourages smoothing along edges in the e1, e2, e1+e2, or −e1+e2 directions and ceases across them.Moreover, we may provide another investigation analysis of our model by examining the eigenvalues-eigenvectors of the matrixD∇u. The matrix D∇u has two eigenvalues λ± with θ± being the corresponding eigenvectors:(18)θ±=η±η±,where(19)η±=gux1−gux2±gux1−gux22+gux12−gux−122gux12−gux−12,provided that ux12≠ux−12. We can then expand the first equation of (16) into(20)∂u∂t=∇⋅λ+θ+θ+T∇u+∇⋅λ−θ−θ−T∇u−ρu−u0.Accordingly, the diffusion caused by (20) is measured by the λ± values and oriented toward θ±. Specifically, it is clear from the expression of λ± that λ+≥λ−>0, which means that the diffusion toward θ+ is privileged over θ−. Furthermore, we can deduce the following results:(i)
In flat areas, we have(21)ux1=ux2,ux12=ux−12⟺ux1=ux2=ux12=ux−12=0,which means thatλ+=λ−=2g0>0. Then, the diffusion is isotropic and linear, and the smoothing effect is the same in all directions.(ii)
At straight edges:λ+>λ−=0.For instance, we can assume thatux1=ux2 and ux1≫0. Then, we obtain ux−12=0 and ux12≫0, which implies that(22)λ+=g0,θ+=12−12T,λ−=0,θ−=−12−12T.Consequently, the diffusion process is anisotropic and oriented along the−e1+e2 direction.(iii) Corners:λ+≥λ−>0.According to the characteristics of the functiong, diffusion is anisotropic and oriented along θ+ and θ− directions.In fact, the differenceλ+−λ−2=gux1−gux22+gux12−gux−122 gives insights into our diffusion model’s anisotropic property. In other words, it indicates the isotropic diffusion for the zero value and the anisotropic diffusion for larger values.


### 3.2. New Adaptive Diffusion Function Using Hermite Spline

For effective image denoising and better control of the diffusion via the process (16), we will approximate numerically the unknown diffusion function g that meets the requirements H1+H2 by using Hermite spline [21]. Therefore, we can use cubic Hermite spline to interpolate numeric data specified at 0, k, and h≥10<k<h to obtain a function that meets the requirements H1+H2. To this end, we propose a function g as follows:(23)gs=p0P1,0ks+v0P2,0ks+pkP1,k0s+vkP2,k0s,s∈0,k,pkP1,khs+vkP2,khs+phP1,hks+vhP2,hks,s∈k,h,phgh,1s+vhgh,2s,s∈h,∞,where p⋅ and v⋅ are the coefficients used to define the position and the velocity vector at a specific point, k and h are two threshold parameters, and Pi,cd is a family of the basis functions composed of polynomials of degree 3 used on the interval c,d such that(24)P1,cds=s−d22s+d−3cd−c3,P2,cds=s−d2s−cd−c2.And, we may consider(25)gh,1s=hlogh+22slogs+1−hloghs2,gh,2s=h2logh+2slogs+1−hlogh+1s2.Hence, we can reformulate the expression ofg as follows:(26)gs=∑i=03A0k,isi,s∈0,k,∑i=03Akh,isi,s∈k,h,Ah,2logs+1s+Ah,11s2,s∈h,∞,By considering the continuity ofφ at k and h, and using φ0=0, it follows then(27)φs=∑i=03A0k,ii+2si+2,s∈0,k,∑i=03Akh,ii+2si+2+A0k,i−Akh,ii+2ki+2,s∈k,h,Ah,2slogs+Ah,1logs+Ah,0,s∈h,∞,where(28)Ah,0=∑i=03Akh,ii+2hi+2+A0k,i−Akh,ii+2ki+2−Ah,2hlogh−Ah,1logh.The values of the coefficientsA0k,. and Akh,. are determined experimentally under the conditions of H1+H2 on the intervals 0,k and k,h. Besides, we may introduce some sufficient conditions on the coefficients Ah,. to ensure that the functions φ and g satisfy H1+H2 on h,∞:(29)Ah,2>0,Ah,1<hAh,2,Ah,1≥−hlogh2Ah,2.Now, we establish a useful growth condition for the functionφ.Proposition 2.
For alls∈ℝ+, we can prove that(30)φs≥As+B,where A>0 and B∈ℝ.Proof.
We know thatg′s≤0. Then,s∈0,k: φ′s=sgs>gks⟹φs≥gks2/2+C≥gks−gk/2+C with C∈ℝ.s∈k,h: φ′s=sgs>kgh⟹φs≥kghs+B with B∈ℝ.s∈h,∞: we have lims⟶+∞φs/s=+∞. Then, we may find an l=h for an A>0 such that φs>As, for all s≥l.


### 3.3. Existence and Uniqueness

In this section, we investigate the existence and uniqueness of weak solutions of the Euler–Lagrange equation associated with the energy functionalEu defined in (5). For ρ>0 and u0∈L2Ω, we consider(31)ρu−u0−divD∇u∇u=0,inΩ,D∇u∇u,n=0,on∂Ω.First, we introduce a new functional spaceLlogLhΩ:(32)LlogLhΩ=u:Ω⟶ℝ|∫Ω∩u≥hulogudx<∞.Next, we define a weak solution for problem (31).Definition 1.
A functionu∈L2Ω∩W1,1Ω with ∂xiu∈LlogLhΩ for i=1,2 is called a weak solution for problem (31) if, for any v∈C1Ω¯ we have(33)∫Ωρu−u0vdx+∫ΩD∇u∇u⋅∇vdx=0.
And whenv is a constant function, we obtain(34)∫Ωudx=∫Ωu0dx.
Before stating our main theorem, we need first to introduce a useful lemma.Lemma 1 (uniform integrability and weak convergence [22]).
Assume thatΩ is bounded and let umm=1∞ be a sequence of functions in L1Ω satisfying(35)supmumL1Ω<∞.
Suppose also(36)liml⟶∞supm∫Ω∩um≥lumdx=0.
Then, there exist a subsequenceumjj=1∞ and u˜∈L1Ω such that(37)umj⇀u˜weakly inL1Ω.
Now, we state our main theorem.Theorem 1.
There exists a unique weak solution for problem (31).Proof.
We consider the variational problem(38)minEu|u∈U,where(39)U=u∈L2Ω∩W1,1Ω|∂xiu∈LlogLhΩwithi=1,2,∫Ωudx=∫Ωu0dx,and the functional E as defined in (5). It is obvious that u0,Ω∈U knowing that u0,Ω=1/Ω∫Ωu0dx.
Since(40)0≤infu∈UEu≤Eu0,Ω=ρ2∫Ωu0,Ω−u02dx,
then we can construct a minimizing sequenceumm=1∞ in U such that Eum<Eu0,Ω+1 and(41)limm⟶∞Eum=infu∈UEu.
Besides,(42)∫Ωum2dx=∫Ωum−u0+u0−u0,Ω+u0,Ω2dx≤4∫Ωum−u02dx+4∫Ωu0−u0,Ω2dx+2∫Ωu0,Ω2dx≤8ρEum+Eu0,Ω+2u0,Ω2Ω≤8ρ2Eu0,Ω+1+2u0,Ω2Ω.
It follows then(43)supmumL2Ω<∞.
On the other hand, from Proposition2, we may find α>0 and β∈ℝ such that(44)Eu0,Ω+1>Eum≥∫Ωφ∂x1um+φ∂x2umdx≥α∇umL1Ω+β.
Then, we get(45)supm∇umL1Ω<∞.
Moreover, we have(46)lims⟶+∞φsslogs=Ah,2>0.
Then, fori=1,2 there exist l≥h≥1 and positive constants C and C1 (independent of l) such that(47)∫Ω∩∂xium≥lφ∂xiumdx≤C∫Ω∩∂xium≥h∂xiumlog∂xiumdx≤C1.
On the other hand, we know that(48)lims⟶+∞φss=+∞.
Then, givenϵ>0, let ηϵ=C1/ϵ and choose l≥h≥1 such that for all s≥l, we have φs>ηϵs. Hence, for i=1,2, we obtain(49)∫Ω∩∂xium≥l∂xiumdx≤1ηϵ∫Ω∩∂xium≥lφ∂xiumdx≤C1ηϵ=ϵ,and this is true for all m and arbitrary ϵ>0. It follows then that, for i=1,2, we have(50)liml⟶∞supm∫Ω∩∂xium≥l∂xiumdx=0.
From (43), (45)–(50) and by Lemma 1 and the weak compactness of L2Ω, we can find a subsequence umjj=1∞ of umm=1∞ and a function u˜∈L2Ω∩W1,1Ω such that(51)umj⇀u˜weakly inL2Ω,∇umj⇀∇u˜weakly inL1Ω.
Therefore, we have(52)∫Ωu˜dx=limj⟶∞∫Ωumjdx=∫Ωu0dx.
Additionally, knowing that the functionfs=slogs for s≥1 is increasing and convex, the function fs is also convex for all s≥1. Therefore, we obtain for i=1,2(53)f∂xiu˜≤f∂xiumj+f′∂xiu˜∂xiu˜−∂xiumj.
Integrating the above inequality overΩM∩h≤∂xiumj with ΩM=Ω∩h≤∂xiu˜≤M with M>0, we have(54)∫ΩMf∂xiu˜dx≤∫Ω∩h≤∂xiumjf∂xiumjdx+∫ΩM∩h≤∂xiumjf′∂xiu˜∂xiu˜−∂xiumjdx≤∫Ω∩h≤∂xiumjf∂xiumjdx+∫Ω∩h≤∂xiumjf′∂xiu˜χh≤∂xiu¯≤M∂xiu˜−∂xiumjdx.
Sincef′∂xiu˜χh≤∂xiu¯≤M∈L∞Ω and by letting j⟶∞, we get(55)∫ΩMf∂xiu˜dx≤liminfj⟶∞∫Ω∩h≤∂xiu˜f∂xiumjdx<∞.
Then, by lettingM⟶∞, we deduce(56)∫Ω∩h≤∂xiu˜f∂xiu˜dx<∞.
It follows then∂xiu˜∈LlogLhΩ. Besides, we havet(57)∫Ωu˜−u02dx≤liminfj⟶∞∫Ωumj−u02dx.
And by following the same reasoning as set out above, we know that the functionφs for s≥0 is increasing and convex. Then, we can easily deduce that(58)∫Ωφ∂xiu˜dx≤liminfj⟶∞∫Ωφ∂xiumjdx.
Therefore, by combining (57) with (58), we conclude(59)Eu˜≤liminfj⟶∞Eumj=infu∈UEu,which signifies that u˜∈U is a minimizer of the energy functional Eu, i.e.,(60)Eu˜=minu∈UEu.
Furthermore, for allv∈C1Ω¯ and for all t∈ℝ, we have u˜+tv−vΩ∈U with vΩ=1/Ω∫Ωvdx. Then, r0≤rt, where(61)rt=Eu˜+tv−vΩ.
Hence, we haver′0=0, which means(62)∫Ωρu˜−u0v−vΩdx+∫ΩD∇u˜∇u˜⋅∇vdx=0.
Because of (52), we have(63)∫Ωρu˜−u0vdx+∫ΩD∇u˜∇u˜⋅∇vdx=0.
We conclude then thatu˜ is a weak solution for problem (31).
Now, assuming that there is another minimizeru^ of E and using the fact that E is strictly convex (5) and (6), we have(64)minu∈UEu=12Eu˜+12Eu^>Eu˜+u^2≥minu∈UEu,a contradiction. Thus, there is only one minimizer, which completes the proof.


## 3.1. Description of the Proposed Model

Our model’s main objective is to allow strong directional smoothing within the areas whereux1, ux2, ux12, or ux−12 is small and prevent blurring boundaries, contours, or corners that separate neighboring areas, where one or a combination of these differential operators has substantial value. Some sufficient assumptions have to be made to ensure that our nonlinear diffusion meets acceptable behavior along these lines. Therefore, we can assume, as exhibited by Perona–Malik [16], some minimal hypotheses on g:(17)H2:g:ℝ+⟶ℝisC1−function such thatgs>0,g′s≤0,lims⟶0gs>0,lims⟶∞gs=0.Besides, due to the decreasing positivity of the functiong, it is evident that our model’s behavior (16) encourages smoothing along edges in the e1, e2, e1+e2, or −e1+e2 directions and ceases across them.Moreover, we may provide another investigation analysis of our model by examining the eigenvalues-eigenvectors of the matrixD∇u. The matrix D∇u has two eigenvalues λ± with θ± being the corresponding eigenvectors:(18)θ±=η±η±,where(19)η±=gux1−gux2±gux1−gux22+gux12−gux−122gux12−gux−12,provided that ux12≠ux−12. We can then expand the first equation of (16) into(20)∂u∂t=∇⋅λ+θ+θ+T∇u+∇⋅λ−θ−θ−T∇u−ρu−u0.Accordingly, the diffusion caused by (20) is measured by the λ± values and oriented toward θ±. Specifically, it is clear from the expression of λ± that λ+≥λ−>0, which means that the diffusion toward θ+ is privileged over θ−. Furthermore, we can deduce the following results:(i)
In flat areas, we have(21)ux1=ux2,ux12=ux−12⟺ux1=ux2=ux12=ux−12=0,which means thatλ+=λ−=2g0>0. Then, the diffusion is isotropic and linear, and the smoothing effect is the same in all directions.(ii)
At straight edges:λ+>λ−=0.For instance, we can assume thatux1=ux2 and ux1≫0. Then, we obtain ux−12=0 and ux12≫0, which implies that(22)λ+=g0,θ+=12−12T,λ−=0,θ−=−12−12T.Consequently, the diffusion process is anisotropic and oriented along the−e1+e2 direction.(iii) Corners:λ+≥λ−>0.According to the characteristics of the functiong, diffusion is anisotropic and oriented along θ+ and θ− directions.In fact, the differenceλ+−λ−2=gux1−gux22+gux12−gux−122 gives insights into our diffusion model’s anisotropic property. In other words, it indicates the isotropic diffusion for the zero value and the anisotropic diffusion for larger values.


## 3.2. New Adaptive Diffusion Function Using Hermite Spline

For effective image denoising and better control of the diffusion via the process (16), we will approximate numerically the unknown diffusion function g that meets the requirements H1+H2 by using Hermite spline [21]. Therefore, we can use cubic Hermite spline to interpolate numeric data specified at 0, k, and h≥10<k<h to obtain a function that meets the requirements H1+H2. To this end, we propose a function g as follows:(23)gs=p0P1,0ks+v0P2,0ks+pkP1,k0s+vkP2,k0s,s∈0,k,pkP1,khs+vkP2,khs+phP1,hks+vhP2,hks,s∈k,h,phgh,1s+vhgh,2s,s∈h,∞,where p⋅ and v⋅ are the coefficients used to define the position and the velocity vector at a specific point, k and h are two threshold parameters, and Pi,cd is a family of the basis functions composed of polynomials of degree 3 used on the interval c,d such that(24)P1,cds=s−d22s+d−3cd−c3,P2,cds=s−d2s−cd−c2.And, we may consider(25)gh,1s=hlogh+22slogs+1−hloghs2,gh,2s=h2logh+2slogs+1−hlogh+1s2.Hence, we can reformulate the expression ofg as follows:(26)gs=∑i=03A0k,isi,s∈0,k,∑i=03Akh,isi,s∈k,h,Ah,2logs+1s+Ah,11s2,s∈h,∞,By considering the continuity ofφ at k and h, and using φ0=0, it follows then(27)φs=∑i=03A0k,ii+2si+2,s∈0,k,∑i=03Akh,ii+2si+2+A0k,i−Akh,ii+2ki+2,s∈k,h,Ah,2slogs+Ah,1logs+Ah,0,s∈h,∞,where(28)Ah,0=∑i=03Akh,ii+2hi+2+A0k,i−Akh,ii+2ki+2−Ah,2hlogh−Ah,1logh.The values of the coefficientsA0k,. and Akh,. are determined experimentally under the conditions of H1+H2 on the intervals 0,k and k,h. Besides, we may introduce some sufficient conditions on the coefficients Ah,. to ensure that the functions φ and g satisfy H1+H2 on h,∞:(29)Ah,2>0,Ah,1<hAh,2,Ah,1≥−hlogh2Ah,2.Now, we establish a useful growth condition for the functionφ.Proposition 2.
For alls∈ℝ+, we can prove that(30)φs≥As+B,where A>0 and B∈ℝ.Proof.
We know thatg′s≤0. Then,s∈0,k: φ′s=sgs>gks⟹φs≥gks2/2+C≥gks−gk/2+C with C∈ℝ.s∈k,h: φ′s=sgs>kgh⟹φs≥kghs+B with B∈ℝ.s∈h,∞: we have lims⟶+∞φs/s=+∞. Then, we may find an l=h for an A>0 such that φs>As, for all s≥l.


## 3.3. Existence and Uniqueness

In this section, we investigate the existence and uniqueness of weak solutions of the Euler–Lagrange equation associated with the energy functionalEu defined in (5). For ρ>0 and u0∈L2Ω, we consider(31)ρu−u0−divD∇u∇u=0,inΩ,D∇u∇u,n=0,on∂Ω.First, we introduce a new functional spaceLlogLhΩ:(32)LlogLhΩ=u:Ω⟶ℝ|∫Ω∩u≥hulogudx<∞.Next, we define a weak solution for problem (31).Definition 1.
A functionu∈L2Ω∩W1,1Ω with ∂xiu∈LlogLhΩ for i=1,2 is called a weak solution for problem (31) if, for any v∈C1Ω¯ we have(33)∫Ωρu−u0vdx+∫ΩD∇u∇u⋅∇vdx=0.
And whenv is a constant function, we obtain(34)∫Ωudx=∫Ωu0dx.
Before stating our main theorem, we need first to introduce a useful lemma.Lemma 1 (uniform integrability and weak convergence [22]).
Assume thatΩ is bounded and let umm=1∞ be a sequence of functions in L1Ω satisfying(35)supmumL1Ω<∞.
Suppose also(36)liml⟶∞supm∫Ω∩um≥lumdx=0.
Then, there exist a subsequenceumjj=1∞ and u˜∈L1Ω such that(37)umj⇀u˜weakly inL1Ω.
Now, we state our main theorem.Theorem 1.
There exists a unique weak solution for problem (31).Proof.
We consider the variational problem(38)minEu|u∈U,where(39)U=u∈L2Ω∩W1,1Ω|∂xiu∈LlogLhΩwithi=1,2,∫Ωudx=∫Ωu0dx,and the functional E as defined in (5). It is obvious that u0,Ω∈U knowing that u0,Ω=1/Ω∫Ωu0dx.
Since(40)0≤infu∈UEu≤Eu0,Ω=ρ2∫Ωu0,Ω−u02dx,
then we can construct a minimizing sequenceumm=1∞ in U such that Eum<Eu0,Ω+1 and(41)limm⟶∞Eum=infu∈UEu.
Besides,(42)∫Ωum2dx=∫Ωum−u0+u0−u0,Ω+u0,Ω2dx≤4∫Ωum−u02dx+4∫Ωu0−u0,Ω2dx+2∫Ωu0,Ω2dx≤8ρEum+Eu0,Ω+2u0,Ω2Ω≤8ρ2Eu0,Ω+1+2u0,Ω2Ω.
It follows then(43)supmumL2Ω<∞.
On the other hand, from Proposition2, we may find α>0 and β∈ℝ such that(44)Eu0,Ω+1>Eum≥∫Ωφ∂x1um+φ∂x2umdx≥α∇umL1Ω+β.
Then, we get(45)supm∇umL1Ω<∞.
Moreover, we have(46)lims⟶+∞φsslogs=Ah,2>0.
Then, fori=1,2 there exist l≥h≥1 and positive constants C and C1 (independent of l) such that(47)∫Ω∩∂xium≥lφ∂xiumdx≤C∫Ω∩∂xium≥h∂xiumlog∂xiumdx≤C1.
On the other hand, we know that(48)lims⟶+∞φss=+∞.
Then, givenϵ>0, let ηϵ=C1/ϵ and choose l≥h≥1 such that for all s≥l, we have φs>ηϵs. Hence, for i=1,2, we obtain(49)∫Ω∩∂xium≥l∂xiumdx≤1ηϵ∫Ω∩∂xium≥lφ∂xiumdx≤C1ηϵ=ϵ,and this is true for all m and arbitrary ϵ>0. It follows then that, for i=1,2, we have(50)liml⟶∞supm∫Ω∩∂xium≥l∂xiumdx=0.
From (43), (45)–(50) and by Lemma 1 and the weak compactness of L2Ω, we can find a subsequence umjj=1∞ of umm=1∞ and a function u˜∈L2Ω∩W1,1Ω such that(51)umj⇀u˜weakly inL2Ω,∇umj⇀∇u˜weakly inL1Ω.
Therefore, we have(52)∫Ωu˜dx=limj⟶∞∫Ωumjdx=∫Ωu0dx.
Additionally, knowing that the functionfs=slogs for s≥1 is increasing and convex, the function fs is also convex for all s≥1. Therefore, we obtain for i=1,2(53)f∂xiu˜≤f∂xiumj+f′∂xiu˜∂xiu˜−∂xiumj.
Integrating the above inequality overΩM∩h≤∂xiumj with ΩM=Ω∩h≤∂xiu˜≤M with M>0, we have(54)∫ΩMf∂xiu˜dx≤∫Ω∩h≤∂xiumjf∂xiumjdx+∫ΩM∩h≤∂xiumjf′∂xiu˜∂xiu˜−∂xiumjdx≤∫Ω∩h≤∂xiumjf∂xiumjdx+∫Ω∩h≤∂xiumjf′∂xiu˜χh≤∂xiu¯≤M∂xiu˜−∂xiumjdx.
Sincef′∂xiu˜χh≤∂xiu¯≤M∈L∞Ω and by letting j⟶∞, we get(55)∫ΩMf∂xiu˜dx≤liminfj⟶∞∫Ω∩h≤∂xiu˜f∂xiumjdx<∞.
Then, by lettingM⟶∞, we deduce(56)∫Ω∩h≤∂xiu˜f∂xiu˜dx<∞.
It follows then∂xiu˜∈LlogLhΩ. Besides, we havet(57)∫Ωu˜−u02dx≤liminfj⟶∞∫Ωumj−u02dx.
And by following the same reasoning as set out above, we know that the functionφs for s≥0 is increasing and convex. Then, we can easily deduce that(58)∫Ωφ∂xiu˜dx≤liminfj⟶∞∫Ωφ∂xiumjdx.
Therefore, by combining (57) with (58), we conclude(59)Eu˜≤liminfj⟶∞Eumj=infu∈UEu,which signifies that u˜∈U is a minimizer of the energy functional Eu, i.e.,(60)Eu˜=minu∈UEu.
Furthermore, for allv∈C1Ω¯ and for all t∈ℝ, we have u˜+tv−vΩ∈U with vΩ=1/Ω∫Ωvdx. Then, r0≤rt, where(61)rt=Eu˜+tv−vΩ.
Hence, we haver′0=0, which means(62)∫Ωρu˜−u0v−vΩdx+∫ΩD∇u˜∇u˜⋅∇vdx=0.
Because of (52), we have(63)∫Ωρu˜−u0vdx+∫ΩD∇u˜∇u˜⋅∇vdx=0.
We conclude then thatu˜ is a weak solution for problem (31).
Now, assuming that there is another minimizeru^ of E and using the fact that E is strictly convex (5) and (6), we have(64)minu∈UEu=12Eu˜+12Eu^>Eu˜+u^2≥minu∈UEu,a contradiction. Thus, there is only one minimizer, which completes the proof.


## 4. Numerical Implementation and Experimental Results

### 4.1. Consistency and Stability of Finite Difference Approximation

For a consistent and stable discretization scheme, one can use the following accurate finite difference scheme approximation: at timetn=nδt, n≥0, and mesh points xi=iδ, yj=jδ0≤i≤N+1and0≤j≤M+1, we denote by ui,jn the finite difference approximation of uxi,yj;tn and exhibit the time-space derivatives discretization as follows:(65)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 assumingδ=1 and denoting(66)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 the solution for problem (16) by the above scheme to obtain the following discrete diffusion filter:(67)ui,jn+1=ui,jn+δtgNΔNu+gEΔEu+gSΔSu+gWΔWu+gNEΔNEu+gSEΔSEu+gSWΔSWu+gNWΔNWu2i,jn−ρδtui,jn−ui,j0,where ui,j0 is the degraded image, 1≤i≤N, 1≤j≤M, and n≥0. Furthermore, we use the discrete Neumann boundary condition:(68)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.Then, using the filter (67) on every initial image u0 yields a unique sequence unn∈ℕ [19]. Besides, due to the continuity of the function g, for every finite n, un depends continuously on u0. Thus, under a specific condition, equation (67) satisfies the following maximum-minimum principle, which describes a stability property for the discrete scheme.Theorem 2.
(discrete extremum principle). For an iteration step satisfying(69)0<δt<1ρ+6maxs∈ℝ+gs,scheme (67) satisfies(70)mini,jui,j0≤ui,jn≤maxi,jui,j0,for all 1≤i≤N, 1≤j≤M, and n∈ℕ.Proof.
For all1≤i≤N, 1≤j≤M, and n∈ℕ, we have(71)ui,jn+1=δtgNi,jnui,j+1n+gEi,jnui+1,jn+gSi,jnui,j−1n+gWi,jnui−1,jn+δt2gNEi,jnui+1,j+1n+gSEi,jnui+1,j−1n+gSWi,jnui−1,j−1n+gNWi,jnui−1,j+1n+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρui,jn+ρδtui,j0.
Since0<δt<1/ρ+6maxs∈ℝ+gs, ρ∈ℝ+∗, and g>0, the following inequality holds:(72)1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρ>0.
It follows then(73)ui,jn+1<δtgNi,jn+gEi,jn+gSi,jn+gWi,jnmaxi,jui,jn+δt2gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jnmaxi,jui,jn+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρmaxi,jui,jn+ρδtmaxi,jui,j0=1−ρδtmaxi,jui,jn+ρδtmaxi,jui,j0.
Similarly, we have(74)ui,jn+1>δtgNi,jn+gEi,jn+gSi,jn+gWi,jnmini,jui,jn+δt2gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jnmini,jui,jn+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρmini,jui,jn+ρδtmini,jui,j0=1−ρδtmini,jui,jn+ρδtmini,jui,j0.
Consequently, the result (70) is true for n=1 because(75)mini,jui,j0=1−ρδtmini,jui,j0+ρδtmini,jui,j0<ui,j1<1−ρδtmaxi,jui,j0+ρδtmaxi,jui,j0=maxi,jui,j0.
Therefore, by induction, it is easy to prove the result (70), which completes the proof.


### 4.2. Experimental Procedures and Results

#### 4.2.1. Experimental Procedures

This section is devoted to comparing our model (16), as an image denoising algorithm, with the ones proposed by Wang and Zhou (Wang–Zhou) [10] and Maiseli [15]. All the experiments are conducted under Windows 10 and GNU Octave version 6.2.0, running on a laptop with an Intel® CoreTM i7-10510U (8 MB cache, 4 Core, up to 4.9 GHz), 16 GB memory (LPDDR3, Dual-channel, 2133 MHz), and 512 GB storage (PCIe, SSD, 3×4). The experiments were done on three actual MRI scans (Figure 1) affected by different σ2-values of zero-mean white Gaussian noise and restored using our filter (67) with the boundary-initial conditions (68), provided that the diffusion function (23) satisfies the assumptions H1+H2.Figure 1
Brain MRI scans (512 × 512). (a) Patient 30: sagittal T1 of a 30-years-old female patient [23]. (b) Patient 50: coronal T2 of a 50-year-old female patient [24]. (c) Patient 55: axial T2 of a 55-year-old patient [25].
(a)(b)(c)We also used the discrete diffusion filter as described in [15] for Wang-Zhou and Maiseli models with the following diffusion functions:Wang–Zhou [10] diffusion function:(76)gs=1s+1+logs+1s,s≥0.Maiseli [15] diffusion function is a combination between the Perona–Malik diffusion function [16] and the Charbonnier diffusion function [8]:(77)gs=11+s/k12,s<k1,11+s/k22,s≥k1,where k1 and k2 are positive constants.Besides, in order to evaluate the quality of the restored images from different image denoising methods, we used two image quality metrics:(i)
Peak signal-to-noise ratioPSNR [26]:(78)PSNR=10log102552MNu−I22,whereI is the original or the uncorrupted image and u is the distorted or the restored image. PSNR is one of the oldest image quality metrics evaluating an image’s signal strength relative to noise, and it is always positive. We evaluate the PSNR metric by using the Octave built-in function “psnr.” However, due to its limitations and its failure in some circumstances as an adequate measure of visual quality [27], we used another metric that significantly correlates with the human visual perception.(ii)
Structural similarity indexSSIM [28]:(79)SSIM=2μuμI+c12σuI+c2μu2+μI2+c1σu2+σI2+c2,wherec1 and c2 are tuning parameters. μu, σu2, and σuI stand for the mean, variance, and covariance, respectively. It is a method for measuring the similarity between a degraded image and a perfect one, and it is bounded between zero and one. For a good similarity between the original and the restored images, we need higher values of SSIM-index. In all experiments, we estimate the SSIM-index using “ssim.m” with automatic downsampling, downloaded from the website: https://ece.uwaterloo.ca/z70wang/research/ssim/.In all experiments and to prioritizeSSIM-index over PSNR, we used a combined metric(80)IQM=PSNR+100×SSIM,to quantify the quality of the restored images. The iterations are stopped when the IQM value reaches the maximum.The image denoising process using the proposed model (67) can be implemented as shown in Algorithm 1. First, we consider I as one of the three original brain MRI scans (Figure 1) and generate in it a Gaussian white noise with zero-mean and σ2-value, using the Octave built-in function “imnoise.” Next, the initial value u0 is set to be the noisy image in the while-loop. After that, a while-loop is executed, and in each iteration, our filter is used to build a new image with a new IQM value. Finally, when the stopping criterion is confirmed, we obtain the restored image with the corresponding iteration.Algorithm 1: Proposed image denoising algorithm.
input:I.output:u, iter, IQM.(1)
Initialize:I;(2)
u0⟵imnoiseI,″gaussian″,0,σ2;(3)
IQM0⟵PSNRI,u0+100×SSIMI,u0;(4)
u⟵u0;(5)
iter⟵0;(6)
IQM⟵IQM0;(7)
n⟵0;(8)
While true do(9)
ΔNun⟵zeros1,M;un1:N−1,:−un2:N,:;(10)
ΔSun⟵−ΔNun2:N,:;zeros1,M;(11)
ΔEun⟵un:,2:M−un:,1:M−1zerosN,1;(12)
ΔWun⟵zerosN,1−ΔEun:,1:M−1;(13)
ΔNEun⟵zeros1,M;un1:N−1,2:M−un2:N,1:M−1zerosN−1,1;(14)
ΔSWun⟵zerosN−1,1−ΔNEun2:N,1:M−1;zeros1,M;(15)
ΔSEun⟵un2:N,2:M−un1:N−1,1:M−1zerosN−1,1;zeros1,M;(16)
ΔNWun⟵zeros1,M;zerosN−1,1−ΔSEun1:N−1,1:M−1;(17)
gαn⟵gαnΔαun,α∈Λ1=N,E,S,W;(18)
gβn⟵gβnΔβun/2,β∈Λ2=NE,SE,SW,NW;(19)
un+1⟵un+δt∑α∈Λ1gαnΔαun+∑β∈Λ2gβnΔβun/2−ρun−u0;(20)
IQMn+1⟵PSNRI,un+1+100×SSIMI,un+1;(21)
ifIQMn+1>IQMthen(22)
IQM⟵IQMn+1;(23)
iter⟵n+1;(24)
u⟵un+1;(25)
else(26)
Break;(27)
end(28)
n⟵n+1;(29)
end


#### 4.2.2. Experimental Results

Several tests have been performed on the images in Figure1 to obtain the best possible parameters:The time step:(81)δt=0.10917.The weighting parameter:(82)ρ=0.00577.The threshold parameters for Maiseli diffusion function (77):(83)k1=2.82261,k2=8.39061.The different parameters for our diffusion function (23):(84)p0=1.46493,v0=−0.00003,k=16.76863,pk=0.60493,vk=−0.00735,h=195.5899,ph=0.15094,vh=−0.00002.Tables1 and 2 show quantitative results on real images, corrupted with various white Gaussian noises, filtered by different filters. From Table 1, it is clear that the PSNR value of the restored images via Maiseli method has, in most cases, higher values with slight differences over the other methods. However, as shown in Table 2, when it comes to the SSIM-index, the proposed method reveals impressive results against the Wang-Zhou and Maiseli methods. Additionally, our method converges more rapidly to the solution than the others with the less iteration number.Table 1
PSNR values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedPSNRPSNRIterPSNRIterPSNRIterPatient 300.00523.440432.84363633.49281933.462870.01020.693530.91825231.50632731.440990.01519.092329.63026430.14883329.903411Patient 500.00523.655131.03132931.34031531.291350.01020.795629.04964229.21572229.066270.01519.140427.82085227.89222827.78358Patient 550.00524.032631.34623031.48971730.629660.01021.233929.19264329.19002428.410580.01519.572627.80865327.74392926.822910Table 2
SSIM-index values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedSSIMSSIMIterSSIMIterSSIMIterPatient 300.0050.70220.9089360.9230190.942070.0100.58890.8649520.8899270.916390.0150.51800.8282640.8619330.896011Patient 500.0050.78320.9139290.9157150.918250.0100.67670.8650420.8702220.873670.0150.60990.8285520.8367280.83978Patient 550.0050.71740.9065300.9155170.931360.0100.60410.8559430.8706240.891380.0150.53580.8138530.8327290.858710From a visual comparison, Figure2 shows that the denoising process using our diffusion function removes noises more efficiently and preserves essential image features. Besides, it can be seen from different images that all the restored images via Wang–Zhou and Maiseli methods are more affected by blocky artifacts than the proposed one, despite the higher values in PSNR. This can be attributed to the fact that the Wang–Zhou and Maiseli methods have lower SSIM-indices.Figure 2
Visual comparison on different actual MRI scans corrupted by zero-mean white Gaussian noise withσ2=0.015 and restored from various methods. (a) Original: patient 30. (b) Original: patient 50. (c) Original: patient 55. (d) Noisy: patient 30. (e) Noisy: patient 50. (f) Noisy: patient 55. (g) Wang–Zhou: patient 30. (h) Wang–Zhou: patient 50. (i) Wang–Zhou: patient 55. (j) Maiseli: patient 30. (k) Maiseli: patient 50. (l) Maiseli: patient 55. (m) Proposed: patient 30. (n) Proposed: patient 50. (o) Proposed: patient 55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 4.1. Consistency and Stability of Finite Difference Approximation

For a consistent and stable discretization scheme, one can use the following accurate finite difference scheme approximation: at timetn=nδt, n≥0, and mesh points xi=iδ, yj=jδ0≤i≤N+1and0≤j≤M+1, we denote by ui,jn the finite difference approximation of uxi,yj;tn and exhibit the time-space derivatives discretization as follows:(65)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 assumingδ=1 and denoting(66)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 the solution for problem (16) by the above scheme to obtain the following discrete diffusion filter:(67)ui,jn+1=ui,jn+δtgNΔNu+gEΔEu+gSΔSu+gWΔWu+gNEΔNEu+gSEΔSEu+gSWΔSWu+gNWΔNWu2i,jn−ρδtui,jn−ui,j0,where ui,j0 is the degraded image, 1≤i≤N, 1≤j≤M, and n≥0. Furthermore, we use the discrete Neumann boundary condition:(68)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.Then, using the filter (67) on every initial image u0 yields a unique sequence unn∈ℕ [19]. Besides, due to the continuity of the function g, for every finite n, un depends continuously on u0. Thus, under a specific condition, equation (67) satisfies the following maximum-minimum principle, which describes a stability property for the discrete scheme.Theorem 2.
(discrete extremum principle). For an iteration step satisfying(69)0<δt<1ρ+6maxs∈ℝ+gs,scheme (67) satisfies(70)mini,jui,j0≤ui,jn≤maxi,jui,j0,for all 1≤i≤N, 1≤j≤M, and n∈ℕ.Proof.
For all1≤i≤N, 1≤j≤M, and n∈ℕ, we have(71)ui,jn+1=δtgNi,jnui,j+1n+gEi,jnui+1,jn+gSi,jnui,j−1n+gWi,jnui−1,jn+δt2gNEi,jnui+1,j+1n+gSEi,jnui+1,j−1n+gSWi,jnui−1,j−1n+gNWi,jnui−1,j+1n+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρui,jn+ρδtui,j0.
Since0<δt<1/ρ+6maxs∈ℝ+gs, ρ∈ℝ+∗, and g>0, the following inequality holds:(72)1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρ>0.
It follows then(73)ui,jn+1<δtgNi,jn+gEi,jn+gSi,jn+gWi,jnmaxi,jui,jn+δt2gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jnmaxi,jui,jn+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρmaxi,jui,jn+ρδtmaxi,jui,j0=1−ρδtmaxi,jui,jn+ρδtmaxi,jui,j0.
Similarly, we have(74)ui,jn+1>δtgNi,jn+gEi,jn+gSi,jn+gWi,jnmini,jui,jn+δt2gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jnmini,jui,jn+1−δtgNi,jn+gEi,jn+gSi,jn+gWi,jn+gNEi,jn+gSEi,jn+gSWi,jn+gNWi,jn2+ρmini,jui,jn+ρδtmini,jui,j0=1−ρδtmini,jui,jn+ρδtmini,jui,j0.
Consequently, the result (70) is true for n=1 because(75)mini,jui,j0=1−ρδtmini,jui,j0+ρδtmini,jui,j0<ui,j1<1−ρδtmaxi,jui,j0+ρδtmaxi,jui,j0=maxi,jui,j0.
Therefore, by induction, it is easy to prove the result (70), which completes the proof.


## 4.2. Experimental Procedures and Results

### 4.2.1. Experimental Procedures

This section is devoted to comparing our model (16), as an image denoising algorithm, with the ones proposed by Wang and Zhou (Wang–Zhou) [10] and Maiseli [15]. All the experiments are conducted under Windows 10 and GNU Octave version 6.2.0, running on a laptop with an Intel® CoreTM i7-10510U (8 MB cache, 4 Core, up to 4.9 GHz), 16 GB memory (LPDDR3, Dual-channel, 2133 MHz), and 512 GB storage (PCIe, SSD, 3×4). The experiments were done on three actual MRI scans (Figure 1) affected by different σ2-values of zero-mean white Gaussian noise and restored using our filter (67) with the boundary-initial conditions (68), provided that the diffusion function (23) satisfies the assumptions H1+H2.Figure 1
Brain MRI scans (512 × 512). (a) Patient 30: sagittal T1 of a 30-years-old female patient [23]. (b) Patient 50: coronal T2 of a 50-year-old female patient [24]. (c) Patient 55: axial T2 of a 55-year-old patient [25].
(a)(b)(c)We also used the discrete diffusion filter as described in [15] for Wang-Zhou and Maiseli models with the following diffusion functions:Wang–Zhou [10] diffusion function:(76)gs=1s+1+logs+1s,s≥0.Maiseli [15] diffusion function is a combination between the Perona–Malik diffusion function [16] and the Charbonnier diffusion function [8]:(77)gs=11+s/k12,s<k1,11+s/k22,s≥k1,where k1 and k2 are positive constants.Besides, in order to evaluate the quality of the restored images from different image denoising methods, we used two image quality metrics:(i)
Peak signal-to-noise ratioPSNR [26]:(78)PSNR=10log102552MNu−I22,whereI is the original or the uncorrupted image and u is the distorted or the restored image. PSNR is one of the oldest image quality metrics evaluating an image’s signal strength relative to noise, and it is always positive. We evaluate the PSNR metric by using the Octave built-in function “psnr.” However, due to its limitations and its failure in some circumstances as an adequate measure of visual quality [27], we used another metric that significantly correlates with the human visual perception.(ii)
Structural similarity indexSSIM [28]:(79)SSIM=2μuμI+c12σuI+c2μu2+μI2+c1σu2+σI2+c2,wherec1 and c2 are tuning parameters. μu, σu2, and σuI stand for the mean, variance, and covariance, respectively. It is a method for measuring the similarity between a degraded image and a perfect one, and it is bounded between zero and one. For a good similarity between the original and the restored images, we need higher values of SSIM-index. In all experiments, we estimate the SSIM-index using “ssim.m” with automatic downsampling, downloaded from the website: https://ece.uwaterloo.ca/z70wang/research/ssim/.In all experiments and to prioritizeSSIM-index over PSNR, we used a combined metric(80)IQM=PSNR+100×SSIM,to quantify the quality of the restored images. The iterations are stopped when the IQM value reaches the maximum.The image denoising process using the proposed model (67) can be implemented as shown in Algorithm 1. First, we consider I as one of the three original brain MRI scans (Figure 1) and generate in it a Gaussian white noise with zero-mean and σ2-value, using the Octave built-in function “imnoise.” Next, the initial value u0 is set to be the noisy image in the while-loop. After that, a while-loop is executed, and in each iteration, our filter is used to build a new image with a new IQM value. Finally, when the stopping criterion is confirmed, we obtain the restored image with the corresponding iteration.Algorithm 1: Proposed image denoising algorithm.
input:I.output:u, iter, IQM.(1)
Initialize:I;(2)
u0⟵imnoiseI,″gaussian″,0,σ2;(3)
IQM0⟵PSNRI,u0+100×SSIMI,u0;(4)
u⟵u0;(5)
iter⟵0;(6)
IQM⟵IQM0;(7)
n⟵0;(8)
While true do(9)
ΔNun⟵zeros1,M;un1:N−1,:−un2:N,:;(10)
ΔSun⟵−ΔNun2:N,:;zeros1,M;(11)
ΔEun⟵un:,2:M−un:,1:M−1zerosN,1;(12)
ΔWun⟵zerosN,1−ΔEun:,1:M−1;(13)
ΔNEun⟵zeros1,M;un1:N−1,2:M−un2:N,1:M−1zerosN−1,1;(14)
ΔSWun⟵zerosN−1,1−ΔNEun2:N,1:M−1;zeros1,M;(15)
ΔSEun⟵un2:N,2:M−un1:N−1,1:M−1zerosN−1,1;zeros1,M;(16)
ΔNWun⟵zeros1,M;zerosN−1,1−ΔSEun1:N−1,1:M−1;(17)
gαn⟵gαnΔαun,α∈Λ1=N,E,S,W;(18)
gβn⟵gβnΔβun/2,β∈Λ2=NE,SE,SW,NW;(19)
un+1⟵un+δt∑α∈Λ1gαnΔαun+∑β∈Λ2gβnΔβun/2−ρun−u0;(20)
IQMn+1⟵PSNRI,un+1+100×SSIMI,un+1;(21)
ifIQMn+1>IQMthen(22)
IQM⟵IQMn+1;(23)
iter⟵n+1;(24)
u⟵un+1;(25)
else(26)
Break;(27)
end(28)
n⟵n+1;(29)
end


### 4.2.2. Experimental Results

Several tests have been performed on the images in Figure1 to obtain the best possible parameters:The time step:(81)δt=0.10917.The weighting parameter:(82)ρ=0.00577.The threshold parameters for Maiseli diffusion function (77):(83)k1=2.82261,k2=8.39061.The different parameters for our diffusion function (23):(84)p0=1.46493,v0=−0.00003,k=16.76863,pk=0.60493,vk=−0.00735,h=195.5899,ph=0.15094,vh=−0.00002.Tables1 and 2 show quantitative results on real images, corrupted with various white Gaussian noises, filtered by different filters. From Table 1, it is clear that the PSNR value of the restored images via Maiseli method has, in most cases, higher values with slight differences over the other methods. However, as shown in Table 2, when it comes to the SSIM-index, the proposed method reveals impressive results against the Wang-Zhou and Maiseli methods. Additionally, our method converges more rapidly to the solution than the others with the less iteration number.Table 1
PSNR values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedPSNRPSNRIterPSNRIterPSNRIterPatient 300.00523.440432.84363633.49281933.462870.01020.693530.91825231.50632731.440990.01519.092329.63026430.14883329.903411Patient 500.00523.655131.03132931.34031531.291350.01020.795629.04964229.21572229.066270.01519.140427.82085227.89222827.78358Patient 550.00524.032631.34623031.48971730.629660.01021.233929.19264329.19002428.410580.01519.572627.80865327.74392926.822910Table 2
SSIM-index values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedSSIMSSIMIterSSIMIterSSIMIterPatient 300.0050.70220.9089360.9230190.942070.0100.58890.8649520.8899270.916390.0150.51800.8282640.8619330.896011Patient 500.0050.78320.9139290.9157150.918250.0100.67670.8650420.8702220.873670.0150.60990.8285520.8367280.83978Patient 550.0050.71740.9065300.9155170.931360.0100.60410.8559430.8706240.891380.0150.53580.8138530.8327290.858710From a visual comparison, Figure2 shows that the denoising process using our diffusion function removes noises more efficiently and preserves essential image features. Besides, it can be seen from different images that all the restored images via Wang–Zhou and Maiseli methods are more affected by blocky artifacts than the proposed one, despite the higher values in PSNR. This can be attributed to the fact that the Wang–Zhou and Maiseli methods have lower SSIM-indices.Figure 2
Visual comparison on different actual MRI scans corrupted by zero-mean white Gaussian noise withσ2=0.015 and restored from various methods. (a) Original: patient 30. (b) Original: patient 50. (c) Original: patient 55. (d) Noisy: patient 30. (e) Noisy: patient 50. (f) Noisy: patient 55. (g) Wang–Zhou: patient 30. (h) Wang–Zhou: patient 50. (i) Wang–Zhou: patient 55. (j) Maiseli: patient 30. (k) Maiseli: patient 50. (l) Maiseli: patient 55. (m) Proposed: patient 30. (n) Proposed: patient 50. (o) Proposed: patient 55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 4.2.1. Experimental Procedures

This section is devoted to comparing our model (16), as an image denoising algorithm, with the ones proposed by Wang and Zhou (Wang–Zhou) [10] and Maiseli [15]. All the experiments are conducted under Windows 10 and GNU Octave version 6.2.0, running on a laptop with an Intel® CoreTM i7-10510U (8 MB cache, 4 Core, up to 4.9 GHz), 16 GB memory (LPDDR3, Dual-channel, 2133 MHz), and 512 GB storage (PCIe, SSD, 3×4). The experiments were done on three actual MRI scans (Figure 1) affected by different σ2-values of zero-mean white Gaussian noise and restored using our filter (67) with the boundary-initial conditions (68), provided that the diffusion function (23) satisfies the assumptions H1+H2.Figure 1
Brain MRI scans (512 × 512). (a) Patient 30: sagittal T1 of a 30-years-old female patient [23]. (b) Patient 50: coronal T2 of a 50-year-old female patient [24]. (c) Patient 55: axial T2 of a 55-year-old patient [25].
(a)(b)(c)We also used the discrete diffusion filter as described in [15] for Wang-Zhou and Maiseli models with the following diffusion functions:Wang–Zhou [10] diffusion function:(76)gs=1s+1+logs+1s,s≥0.Maiseli [15] diffusion function is a combination between the Perona–Malik diffusion function [16] and the Charbonnier diffusion function [8]:(77)gs=11+s/k12,s<k1,11+s/k22,s≥k1,where k1 and k2 are positive constants.Besides, in order to evaluate the quality of the restored images from different image denoising methods, we used two image quality metrics:(i)
Peak signal-to-noise ratioPSNR [26]:(78)PSNR=10log102552MNu−I22,whereI is the original or the uncorrupted image and u is the distorted or the restored image. PSNR is one of the oldest image quality metrics evaluating an image’s signal strength relative to noise, and it is always positive. We evaluate the PSNR metric by using the Octave built-in function “psnr.” However, due to its limitations and its failure in some circumstances as an adequate measure of visual quality [27], we used another metric that significantly correlates with the human visual perception.(ii)
Structural similarity indexSSIM [28]:(79)SSIM=2μuμI+c12σuI+c2μu2+μI2+c1σu2+σI2+c2,wherec1 and c2 are tuning parameters. μu, σu2, and σuI stand for the mean, variance, and covariance, respectively. It is a method for measuring the similarity between a degraded image and a perfect one, and it is bounded between zero and one. For a good similarity between the original and the restored images, we need higher values of SSIM-index. In all experiments, we estimate the SSIM-index using “ssim.m” with automatic downsampling, downloaded from the website: https://ece.uwaterloo.ca/z70wang/research/ssim/.In all experiments and to prioritizeSSIM-index over PSNR, we used a combined metric(80)IQM=PSNR+100×SSIM,to quantify the quality of the restored images. The iterations are stopped when the IQM value reaches the maximum.The image denoising process using the proposed model (67) can be implemented as shown in Algorithm 1. First, we consider I as one of the three original brain MRI scans (Figure 1) and generate in it a Gaussian white noise with zero-mean and σ2-value, using the Octave built-in function “imnoise.” Next, the initial value u0 is set to be the noisy image in the while-loop. After that, a while-loop is executed, and in each iteration, our filter is used to build a new image with a new IQM value. Finally, when the stopping criterion is confirmed, we obtain the restored image with the corresponding iteration.Algorithm 1: Proposed image denoising algorithm.
input:I.output:u, iter, IQM.(1)
Initialize:I;(2)
u0⟵imnoiseI,″gaussian″,0,σ2;(3)
IQM0⟵PSNRI,u0+100×SSIMI,u0;(4)
u⟵u0;(5)
iter⟵0;(6)
IQM⟵IQM0;(7)
n⟵0;(8)
While true do(9)
ΔNun⟵zeros1,M;un1:N−1,:−un2:N,:;(10)
ΔSun⟵−ΔNun2:N,:;zeros1,M;(11)
ΔEun⟵un:,2:M−un:,1:M−1zerosN,1;(12)
ΔWun⟵zerosN,1−ΔEun:,1:M−1;(13)
ΔNEun⟵zeros1,M;un1:N−1,2:M−un2:N,1:M−1zerosN−1,1;(14)
ΔSWun⟵zerosN−1,1−ΔNEun2:N,1:M−1;zeros1,M;(15)
ΔSEun⟵un2:N,2:M−un1:N−1,1:M−1zerosN−1,1;zeros1,M;(16)
ΔNWun⟵zeros1,M;zerosN−1,1−ΔSEun1:N−1,1:M−1;(17)
gαn⟵gαnΔαun,α∈Λ1=N,E,S,W;(18)
gβn⟵gβnΔβun/2,β∈Λ2=NE,SE,SW,NW;(19)
un+1⟵un+δt∑α∈Λ1gαnΔαun+∑β∈Λ2gβnΔβun/2−ρun−u0;(20)
IQMn+1⟵PSNRI,un+1+100×SSIMI,un+1;(21)
ifIQMn+1>IQMthen(22)
IQM⟵IQMn+1;(23)
iter⟵n+1;(24)
u⟵un+1;(25)
else(26)
Break;(27)
end(28)
n⟵n+1;(29)
end


## 4.2.2. Experimental Results

Several tests have been performed on the images in Figure1 to obtain the best possible parameters:The time step:(81)δt=0.10917.The weighting parameter:(82)ρ=0.00577.The threshold parameters for Maiseli diffusion function (77):(83)k1=2.82261,k2=8.39061.The different parameters for our diffusion function (23):(84)p0=1.46493,v0=−0.00003,k=16.76863,pk=0.60493,vk=−0.00735,h=195.5899,ph=0.15094,vh=−0.00002.Tables1 and 2 show quantitative results on real images, corrupted with various white Gaussian noises, filtered by different filters. From Table 1, it is clear that the PSNR value of the restored images via Maiseli method has, in most cases, higher values with slight differences over the other methods. However, as shown in Table 2, when it comes to the SSIM-index, the proposed method reveals impressive results against the Wang-Zhou and Maiseli methods. Additionally, our method converges more rapidly to the solution than the others with the less iteration number.Table 1
PSNR values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedPSNRPSNRIterPSNRIterPSNRIterPatient 300.00523.440432.84363633.49281933.462870.01020.693530.91825231.50632731.440990.01519.092329.63026430.14883329.903411Patient 500.00523.655131.03132931.34031531.291350.01020.795629.04964229.21572229.066270.01519.140427.82085227.89222827.78358Patient 550.00524.032631.34623031.48971730.629660.01021.233929.19264329.19002428.410580.01519.572627.80865327.74392926.822910Table 2
SSIM-index values of the images in Figure 1 are affected by different σ2 values of zero-mean white Gaussian noise and restored from different methods with their corresponding iteration number.
σ2NoisyWang and Zhou [10]Maiseli [15]ProposedSSIMSSIMIterSSIMIterSSIMIterPatient 300.0050.70220.9089360.9230190.942070.0100.58890.8649520.8899270.916390.0150.51800.8282640.8619330.896011Patient 500.0050.78320.9139290.9157150.918250.0100.67670.8650420.8702220.873670.0150.60990.8285520.8367280.83978Patient 550.0050.71740.9065300.9155170.931360.0100.60410.8559430.8706240.891380.0150.53580.8138530.8327290.858710From a visual comparison, Figure2 shows that the denoising process using our diffusion function removes noises more efficiently and preserves essential image features. Besides, it can be seen from different images that all the restored images via Wang–Zhou and Maiseli methods are more affected by blocky artifacts than the proposed one, despite the higher values in PSNR. This can be attributed to the fact that the Wang–Zhou and Maiseli methods have lower SSIM-indices.Figure 2
Visual comparison on different actual MRI scans corrupted by zero-mean white Gaussian noise withσ2=0.015 and restored from various methods. (a) Original: patient 30. (b) Original: patient 50. (c) Original: patient 55. (d) Noisy: patient 30. (e) Noisy: patient 50. (f) Noisy: patient 55. (g) Wang–Zhou: patient 30. (h) Wang–Zhou: patient 50. (i) Wang–Zhou: patient 55. (j) Maiseli: patient 30. (k) Maiseli: patient 50. (l) Maiseli: patient 55. (m) Proposed: patient 30. (n) Proposed: patient 50. (o) Proposed: patient 55.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 5. Conclusion

This paper presented a new anisotropic diffusion model for image restoration, eliminating the corruption caused by white Gaussian noise. The existence and uniqueness of weak solutions of functions in Orlicz–Sobolev space that minimize the energy functionalE have been proven, provided that the functions φ and g satisfy some specific conditions H1+H2. Besides, we have established a consistent and stable numerical model for the denoising process and used the cubic Hermite spline to approximate the best possible diffusion function that revealed its efficiency regarding the optimal image denoising via (67). We have also proved that our method provides better results compared to the Wang–Zhou and Maiseli methods.The next stage of the research attempts to prove the existence and uniqueness of weak solutions for the evolution problem (16) and use different numerical methods, such as finite element and mixed finite element methods.



---

*Source: 6691795-2021-04-21.xml*