# Perona-Malik Model with Diffusion Coefficient Depending on Fractional Gradient via Caputo-Fabrizio Derivative

Authors: Gustavo Asumu Mboro Nchama; Angela Leon Mecias; Mariano Rodriguez Ricard
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

The Perona-Malik (PM) model is used successfully in image processing to eliminate noise while preserving edges; however, this model has a major drawback: it tends to make the image look blocky. This work proposes to modify the PM model by introducing the Caputo-Fabrizio fractional gradient inside the diffusivity function. Experiments with natural images show that our model can suppress efficiently the blocky effect. Also, our model has good performance in visual quality, high peak signal-to-noise ratio (PSNR), and lower value of mean absolute error (MAE) and mean square error (MSE).


---

## Body

## 1. Introduction and Some Basic Definitions

Image processing based on partial differential equations (PDEs) is mainly used for smoothing and restoration purposes. Typical PDE techniques for image smoothing regard the original image as initial states of a parabolic process and extract filtered versions from its temporal evolution. In this paper,ux,y,t represents the image intensity values in the position x,y∈Ω⊂ℝ2 for a time t, ∂Ω is the smooth boundary, and u0x,y is the original image. It is a classical result that for any bounded u∈Cℝ2, the linear diffusion process
(1)∂ux,y,t∂t=Δux,y,t,Ω×0,T,(2)ux,y,0=u0x,y,∂Ω,possesses the solution
(3)ux,y,t=u0x,yt=0,k2t∗u0x,yt>0.For more details, see [1]. However, k can be a smooth kernel or a blur kernel. In 1960, Gabor remarked that the difference between the original and the blurred image is roughly proportional to its Laplacian [2]. To formalize this remark, we have to notice that k is spatially concentrated and that we may introduce a scale parameter for k, namely [2],
(4)khx=1hkxh1/2,then
(5)k∗u0x−u0xh→Δu0x.Here,x=x,y∈ℝ2. When h gets smaller, the blurring process looks more and more like the heat equation [2]. The linear diffusion equation (1) does not only smooth noise, but it also blurs important features such as edges and, thus, makes them harder to identify. To remove the noise while preserving the edges at best, Perona and Malik proposed the following nonlinear diffusion method (called by them anisotropic diffusion) [3]:
(6)∂ux,y,t∂t−divg∇ux,y,t∇ux,y,t=0,Ω×0,T,(7)∂ux,y,t∂n=0,∂Ω×0,T,(8)ux,y,0=u0x,y,Ω,where ux,y,0 is the original image intensity function, ux,y,t is the smoothed image intensity function in time t, ∇ denotes the gradient, div· is the divergence operator, ∂ux,y,t/∂n denotes the derivative in normal direction to the boundary, and g∣∇ux,y,t∣ is the diffusion coefficient function, also called edge-stopping function. The edge-stopping function gs is a nonnegative decreasing function satisfying two conditions. One is that gs=1 as s→0, so that the rate of diffusion is high within uniform or inner regions, and the other one is gs→0 as s→∞, so that the diffusion is totally zero across boundaries. The important property of edge functions is that they should have an insignificant value for those gradients that correspond to edges [3]. In the literature, different edge-stopping functions have been proposed. For example, in [3], Perona-Malik used
(9)g1s=e−s/K2,(10)g2s=11+s/K2,where s=∣∇ux,y,t∣ and K is the gradient magnitude threshold parameter that decides the amount of diffusion that take place (see for instance a way to estimate it in [4]). Other expressions for the diffusion coefficient can be found in [1, 5, 6].The flow functionϕ defined as
(11)ϕs=sgs,represents the sum of the brightness flow that is generated. The maximum flow is generated at locations where ∣∇u∣=K [7].The Perona-Malik equation (6) and the large amount of its modifications [5, 8–11] have demonstrated to be able to achieve a good trade-off between noise removal and edge preservation. Unfortunately, edge-stopping functions lead to backward-forward problems that are ill-posed [1]; for this reason, in [8, 12], authors introduced a modification in the diffusion coefficient to obtain a regularized version as follows:
(12)∂ux,y,t∂t−divg∇Gσ∗ux,y,t∇ux,y,t=0,Ω×0,T,(13)∂ux,y,t∂n=0,∂Ω×0,T,(14)ux,y,0=u0x,y,Ω,where Ω is a bounded domain of ℝ2 with an appropriately smooth boundary, n denotes the unit outer normal to Ω, and T>0. The term ∇Gσ∗ux,y,t is the regularized version of ∇ux,y,t, and it is the gradient of the solution at time t=σ2/2 of the linear heat equation (1). Despite the ill-posedness, the practical results seem to be quite good in most cases. As reported by Weickert in [1], the reason is that the numerical scheme used by Perona and Malik does not correspond to their equation but rather to a time-regularized one which is well-posed this time. Numerically, the mainly observable instability is the so-called staircase (blocky) effect, where a sigmoid edge evolves into piecewise linear segments that are separated by jumps. This effect is visually unpleasant and is likely to cause a computer vision system to falsely recognize as edges the boundaries of different blocks that belong to the same smooth area in the original image [13].To avoid staircase effects while achieving a good trade-off between noise removal and edge preservation, many authors have proposed to use fourth-order partial differential equations [13]. For example, You and Kaveh proposed in [13–17] a fourth-order PDE for noise removal and an algorithm that can remove speckle effect efficiently. In [15], the authors studied the fourth-order telegraph-diffusion equation for image restoration in which the numerical method not only preserves the edge but also avoids the staircase effects, and the existence, uniqueness, and stability were stabilized. The experiments show that all these fourth-order models can improve the peak signal-to-noise ratio, preserve texture, and eliminate the staircase (blocky) effect efficiently (see [18] and the references therein). The problem with the use of fourth-order equations is that they tend to leave the image with isolated black and white speckles (so-called speckle effect). This effect is characterized as pixels whose intensity values are either much larger or much smaller than those of the neighboring pixels [13, 19]. To eliminate the undesirable speckle effect, in [20], the energy functional associated with the Perona-Malik equation is redefined as an increasing function of the absolute value of the image intensity fractional derivative function (defined in the Fourier domain). The corresponding Euler-Lagrange equation is then a fractional-order anisotropic diffusion equation. The proposed pseudo-PDEs will lead to an interpolation between Perona-Malik equations and fourth-order anisotropic diffusion equations in [13]. In [18], a fractional-order Perona-Malik diffusion (FOPMD) equation is proposed. In this model, the integer-order derivative concerning spatial variables of the Perona-Malik diffusion is replaced with the fractional-order Grunwald-Letnikov derivatives (Polubny). The FOPMD model is given by
(15)∂u∂tx,y,t=divαg∣∇αux,y,t∣∇αux,y,t,Ω×0,T,(16)∂ux,y,t∂n=0,∂Ω×0,T,(17)ux,y,0=u0x,y,Ω,where the fractional order is α (0≤α≤2) and ux,y,t is the smooth gray scale image at time t. The fractional-order gradient vector with α order is defined as
(18)∇αux,y,t=∇xαux,y,t,∇yαux,y,t,where ∇pαux,y,t represents the partial-order derivative of ux,y,t with respect to the variable p whose order is α. Although model (15) has reported good performance of preserving edges and suppressing staircase and speckle effect, the resulting images still have some artifacts. Therefore, in [21], a new diffusion model named regularized fully spatial fractional-order Perona-Malik diffusion (RFS-FOPMD) is proposed, given by
(19)∂ux,y,t∂t=divαgGσ·∇αux,y,t∇αux,y,t,Ω×0,T,(20)∂ux,y,t∂n=0,∂Ω×0,T,(21)ux,y,0=u0x,y,Ω.As observed, equation (19) is obtained by substituting in equation (12), the integer spatial derivative for fractional-order derivative according to Grunwald-Letnikov. The fractional-order α is taken as 0.5≤α≤1.5. With this method, results obtained in [18] are improved, partly by improving the performance of edges locating by regularization [21]. Despite this, there are still errors in locating the position of the edges. In [22], it is said that the basic reason for locating edge positions falsely is that the fractional-order gradient module cannot be used as an edge indicator; therefore, the authors in [22] adopted a so-called external fractional-order gradient vector Perona-Malik diffusion by only replacing integer-order derivatives of the external gradient vector to their fractional-order counterparts while keeping first-order derivatives for diffusion coefficient. The model is the same as [9] except for the derivative used in the diffusion coefficient. A novel fractional PDE model is given by Guidotti and Longo in [23]; they address the well-posedness of the following fourth-order model for noise removal, using fractional derivatives defined by Fourier transform.
(22)∂ux,y,t∂t+Δg∇1−εux,y,tΔux,y,t=0,Ω×0,T,(23)∂ux,y,t∂n=0,∂Ω×0,T,(24)ux,y,0=u0x,y,Ω,with u periodic and ε∈0,1. Model (22) is a modification of the You and Kaveh model [13]. As numerical approximation of the equation, they used a scheme based on the Krylov subspace spectral method.There are many definitions of fractional derivatives (three popular definitions were given by Grunwald-Letnikov (G-L), Riemann-Liouville (R-L), and Caputo). These have been used in numerous fields of science such as study of the anomalous diffusion phenomenon [24–26], circuit theory [27–29], and image processing [30, 31], among other applications [11, 32–48]. Given the discussion above, we consider that using anisotropic diffusion models to eliminate noise in an image, preserving both strong and weak edges and without phenomena such as staircase, speckle, or any type of artifact, is a subject where much remains to be investigated. In image processing, integer-order differentiation operators are often used in edge detection, especially the first order for the gradient (e.g., Roberts, Prewitt, and Sobel) and second order for the Laplacian (e.g., Laplacian of Gaussian). However, the first-order derivative methods generally cause thicker edges, resulting in the loss of image details. The second-order derivative methods have a stronger response to fine detail, but they are more sensitive to noise [49]. To solve this problem, the fractional-order derivative has been introduced in the edge detection methods, with the capability to preserve more low-frequency contour features in the smooth areas, maintain high-frequency marginal features, and also enhance medium-frequency texture details [49, 50]. Many fractional-order operators are used for edge detection, such as the fractional-order Sobel operator [31], fractional-order CRONE operator [51], and quaternion fractional differential operator [52].Inspired by these ideas, and taking into account that the use of fractional derivatives inside the diffusivity function has proven to be robust in the presence of noise as said in [53], we present a modification to the Perona-Malik model. Instead of the diffusion coefficient controlling the diffusion from the threshold of the integer-order gradient, it does so according to the fractional-order gradient threshold. Our goal by incorporating the fractional gradient into the diffusion coefficient is to avoid undesirable artifacts such as blocky effect while removing noise and preserving edges. To the best of our knowledge, from the consulted literature, this model using the recently defined Caputo-Fabrizio derivative has not been used. Caputo and Fabrizio introduced a new fractional derivative in [54], intending to describe structures with different scales.We propose to insert the Caputo-Fabrizio fractional gradient∇αux,y,t into the diffusivity function as follows:
(25)∂u∂tx,y,t=div∇ux,y,t·g∇αux,y,t,Ω×0,T,(26)∂ux,y,t∂n=0,Ω×0,T,(27)ux,y,0=u0x,y,Ω.The flow functionϕ corresponding to equation (25) is a function depending on two variables and takes the following form:
(28)ϕs,v=sgv,where g is given by (9) and s=∇ux,y,t and v=∇αux,y,t. This means that the graph of function (28) is a surface. For simplicity reasons, the integer-order gradient ∇ux,y,t will be fixed by a real positive value while the fractional-order gradient ∇αux,y,t will be taken as the independent variable. The simulation of function (28) with different values of K produces pictures related to Figure 1. From these simulations, we observe that for v→0, the flow function takes higher values. Otherwise, when v takes great values, we see a different behavior.Figure 1
Graphics of flow function. Integer-ordergradient=5. (a) Flux function for k=2. (b) Flux function for k=4. (c) Flux function for k=6. (d) Flux function for k=8.
(a)(b)(c)(d)


### 1.1. Some Basic Definitions in Fractional Calculus

Leta,b,p∈ℝ, such that a<b and 1≤p≤∞. We denote by C0a,b the space of all continuous functions on a,b with compact support. We denote C0ka,b as
(29)C0ka,b=Cka,b∩C0a,b,k∈ℤ+∪∞.The Sobolev spaceW1,pa,b is defined by W1,pa,b=u∈Lpa,b;∃g∈Lpa,bsuchthat∫abuφ′=−∫abgφ,∀φ∈C0∞a,b.We set(30)H1a,b=W1,2a,b.Three popular definitions of fractional calculus were given by Grunwald-Letnikov (G-L), Riemann-Liouville (R-L), and Caputo. Of these, G-L and R-L are the most popular definitions used in digital image processing. Recently, Caputo and Fabrizio introduced a new fractional derivative in [54]. Considering the spatial variable, it is defined as
(31)D0xαux=11−α∫0xe−α/1−αx−ξu′ξdξ,with u∈H1a,b,b>a and α∈0,1. For more details, see [4, 54–59].The remainder of this paper is organized as follows: Section2 presents discretization of numerical schemes. Experimental results are provided in Section 3. Finally, conclusions are summarized in Section 4.


## 1.1. Some Basic Definitions in Fractional Calculus

Leta,b,p∈ℝ, such that a<b and 1≤p≤∞. We denote by C0a,b the space of all continuous functions on a,b with compact support. We denote C0ka,b as
(29)C0ka,b=Cka,b∩C0a,b,k∈ℤ+∪∞.The Sobolev spaceW1,pa,b is defined by W1,pa,b=u∈Lpa,b;∃g∈Lpa,bsuchthat∫abuφ′=−∫abgφ,∀φ∈C0∞a,b.We set(30)H1a,b=W1,2a,b.Three popular definitions of fractional calculus were given by Grunwald-Letnikov (G-L), Riemann-Liouville (R-L), and Caputo. Of these, G-L and R-L are the most popular definitions used in digital image processing. Recently, Caputo and Fabrizio introduced a new fractional derivative in [54]. Considering the spatial variable, it is defined as
(31)D0xαux=11−α∫0xe−α/1−αx−ξu′ξdξ,with u∈H1a,b,b>a and α∈0,1. For more details, see [4, 54–59].The remainder of this paper is organized as follows: Section2 presents discretization of numerical schemes. Experimental results are provided in Section 3. Finally, conclusions are summarized in Section 4.


## 2. Numerical Schemes

The main purpose of this section is to give an explicit finite difference scheme of the proposed model. For it, we will start discretizing the Caputo-Fabrizio fractional derivative. Next, we will recall the numerical approximation for the PM equation.


### 2.1. Discretization of the Caputo-Fabrizio Derivative

To solve numerically the new fractional Perona-Malik model, we propose to discretize the Caputo-Fabrizio derivative based on the forward finite difference scheme in the interval0,x (analogously 0,y). Let us take a partition of N nodes of the interval 0,x, with step Δx=x/N and α∈0,1; then, we obtain
(32)D0xαu=11−α∫0xe−α/1−αx−ξu′ξdξ=11−α∑k=0N−1∫kx/Nk+1x/Ne−α/1−αx−ξu′ξkdξ≈11−α∑k=0N−1ukx+x/N−ukx/NΔx·∫kx/Nk+1x/Ne−α/1−αx−ξdξ=1α∑k=0N−1uk+1Δx−ukΔxΔx·e−α/1−αN−k+1Δx−e−α/1−αN−kΔx.As in digital gray imageux,y, the shortest distance of 2-D image on x and y coordinates is one pixel; then, we put Δx=1, and from (32), we obtain
(33)D0xαu≈1αe−α/1−αNeα/1−α−1∑k=0N−1eα/1−αk·uk+1−uk=∑k=0N−1Cα,k·uk+1−uk,where
(34)Cα,k=1αe−α/1−αNeα/1−α−1eα/1−αk.An analogue expression is obtained forD0yαu.


### 2.2. Perona-Malik Scheme

Here, we recall the numerical approximation for the PM equation. LetΔx and Δy be the space steps such that Δx=Δy=1, and denoting the time step as Δt, the discretization of the PM equation
(35)∂ux,y,t∂t=divg∇u∇ux,y,t,is given in [13], as
(36)ui,j,t+1=ui,j,t+Δt·gN·∇Nui,j,t+gS·∇Sui,j,t+gE·∇Eui,j,t+gW·∇Wui,j,t,where
(37)∇Nui,j,t=ui−1,j,t−ui,j,t,(38)∇Sui,j,t=ui+1,j,t−ui,j,t,(39)∇Eui,j,t=ui,j+1,t−ui,j,t,(40)∇Wui,j,t=ui,j−1,t−ui,j,t,(41)gP=g∇Pui,j,P=N,S,E,W,with 0<Δt≤1/4 for the numerical scheme to be stable.


### 2.3. Numerical Method of the Proposed Model

As we said above, we propose to introduce the Caputo-Fabrizio fractional gradient∇αux,y,t inside the diffusivity function of the Perona-Malik model
(42)∂u∂tx,y,t=divg∇αux,y,t∇ux,y,t,where
(43)g∇u=11+∇u/K2,is used in this work. Considering that
(44)∇u=∂u∂x2+∂u∂y2,then the discretization of ∣∇u∣, using the forward finite difference scheme, is given by
(45)gi,j=11+1/K2ui+1j−uij2+uij+1−uij2.The discrete fractional-order gradient∇αu is an 8-dimensional vector, as in [18]:
(46)∇αui,j=∇Eαui,j,∇SEαui,j,∇Sαui,j,∇SWαui,j,∇Wαui,j,∇NWαui,j,∇Nαui,j,∇NEαui,jT,where T represents the transpose of the vector and each component is defined as
(47)∇Eαui,j≔∑k=0N−1Cα,k·ui,j+1,∇SEαui,j≔∑k=0N−1Cα,k·ui+k,j+k,(48)∇Sαui,j≔∑k=0N−1Cα,k·ui+k,j,∇SWαui,j≔∑k=0N−1Cα,k·ui+k,j−k,(49)∇Wαui,j≔∑k=0N−1Cα,k·ui,j−k,∇NWαui,j≔∑k=0N−1Cα,k·ui−k,j−k,(50)∇Nαui,j≔∑k=0N−1Cα,k·ui−k,j,∇NEαui,j≔∑k=0N−1Cα,k·ui−k,j+k,where Cα,k is given by (34). The discretization of the usual gradient, following the idea in [3], in eight directions, is given by
(51)∇ui,j=∇Eui,j,∇SEui,j,∇Sui,j,∇SWui,j,∇Wui,j,∇NWui,j,∇Nui,j,∇NEui,jT,where ∇Nu, ∇Su, ∇Eu, and ∇Wu are defined as in (37), (38), (39), and (40), respectively, and
(52)∇SEui,j,t=ui+1,j+1,t−ui,j,t,∇SWui,j,t=ui+1,j−1,t−ui,j,t,∇NWui,j,t=ui−1,j−1,t−ui,j,t,∇NEui,j,t=ui−1,j+1,t−ui,j,t.The subscripts N, S, E, W, NE, SW, NW, and SE denote the eight directions North, South, East, West, North-East, South-West, North-West, and South-East. Hence, the explicit discretization scheme for equation (42) is given by
(53)ui,j,t+1=ui,j,t+18·Δt·gEα·∇Eui,j,t+gSEα·∇SEui,j,t+gSα·∇Sui,j,t+gSW·∇SWui,j,t+gWα·∇Wui,j,t+gNWα·∇NWui,j,t+gNα·∇Nui,j,t+gNEα·∇NEui,j,t,where
(54)gEα=g∇Eαui,j,gSEα=g∇SEαui,j,gSα=g∇Sαui,j,gSWα=g∇SWαui,j,gWα=g∇Wαui,j,gNWα=g∇NWαui,j,gNα=g∇Nαui,j,gNEα=g∇NEαui,j.


## 2.1. Discretization of the Caputo-Fabrizio Derivative

To solve numerically the new fractional Perona-Malik model, we propose to discretize the Caputo-Fabrizio derivative based on the forward finite difference scheme in the interval0,x (analogously 0,y). Let us take a partition of N nodes of the interval 0,x, with step Δx=x/N and α∈0,1; then, we obtain
(32)D0xαu=11−α∫0xe−α/1−αx−ξu′ξdξ=11−α∑k=0N−1∫kx/Nk+1x/Ne−α/1−αx−ξu′ξkdξ≈11−α∑k=0N−1ukx+x/N−ukx/NΔx·∫kx/Nk+1x/Ne−α/1−αx−ξdξ=1α∑k=0N−1uk+1Δx−ukΔxΔx·e−α/1−αN−k+1Δx−e−α/1−αN−kΔx.As in digital gray imageux,y, the shortest distance of 2-D image on x and y coordinates is one pixel; then, we put Δx=1, and from (32), we obtain
(33)D0xαu≈1αe−α/1−αNeα/1−α−1∑k=0N−1eα/1−αk·uk+1−uk=∑k=0N−1Cα,k·uk+1−uk,where
(34)Cα,k=1αe−α/1−αNeα/1−α−1eα/1−αk.An analogue expression is obtained forD0yαu.


## 2.2. Perona-Malik Scheme

Here, we recall the numerical approximation for the PM equation. LetΔx and Δy be the space steps such that Δx=Δy=1, and denoting the time step as Δt, the discretization of the PM equation
(35)∂ux,y,t∂t=divg∇u∇ux,y,t,is given in [13], as
(36)ui,j,t+1=ui,j,t+Δt·gN·∇Nui,j,t+gS·∇Sui,j,t+gE·∇Eui,j,t+gW·∇Wui,j,t,where
(37)∇Nui,j,t=ui−1,j,t−ui,j,t,(38)∇Sui,j,t=ui+1,j,t−ui,j,t,(39)∇Eui,j,t=ui,j+1,t−ui,j,t,(40)∇Wui,j,t=ui,j−1,t−ui,j,t,(41)gP=g∇Pui,j,P=N,S,E,W,with 0<Δt≤1/4 for the numerical scheme to be stable.


## 2.3. Numerical Method of the Proposed Model

As we said above, we propose to introduce the Caputo-Fabrizio fractional gradient∇αux,y,t inside the diffusivity function of the Perona-Malik model
(42)∂u∂tx,y,t=divg∇αux,y,t∇ux,y,t,where
(43)g∇u=11+∇u/K2,is used in this work. Considering that
(44)∇u=∂u∂x2+∂u∂y2,then the discretization of ∣∇u∣, using the forward finite difference scheme, is given by
(45)gi,j=11+1/K2ui+1j−uij2+uij+1−uij2.The discrete fractional-order gradient∇αu is an 8-dimensional vector, as in [18]:
(46)∇αui,j=∇Eαui,j,∇SEαui,j,∇Sαui,j,∇SWαui,j,∇Wαui,j,∇NWαui,j,∇Nαui,j,∇NEαui,jT,where T represents the transpose of the vector and each component is defined as
(47)∇Eαui,j≔∑k=0N−1Cα,k·ui,j+1,∇SEαui,j≔∑k=0N−1Cα,k·ui+k,j+k,(48)∇Sαui,j≔∑k=0N−1Cα,k·ui+k,j,∇SWαui,j≔∑k=0N−1Cα,k·ui+k,j−k,(49)∇Wαui,j≔∑k=0N−1Cα,k·ui,j−k,∇NWαui,j≔∑k=0N−1Cα,k·ui−k,j−k,(50)∇Nαui,j≔∑k=0N−1Cα,k·ui−k,j,∇NEαui,j≔∑k=0N−1Cα,k·ui−k,j+k,where Cα,k is given by (34). The discretization of the usual gradient, following the idea in [3], in eight directions, is given by
(51)∇ui,j=∇Eui,j,∇SEui,j,∇Sui,j,∇SWui,j,∇Wui,j,∇NWui,j,∇Nui,j,∇NEui,jT,where ∇Nu, ∇Su, ∇Eu, and ∇Wu are defined as in (37), (38), (39), and (40), respectively, and
(52)∇SEui,j,t=ui+1,j+1,t−ui,j,t,∇SWui,j,t=ui+1,j−1,t−ui,j,t,∇NWui,j,t=ui−1,j−1,t−ui,j,t,∇NEui,j,t=ui−1,j+1,t−ui,j,t.The subscripts N, S, E, W, NE, SW, NW, and SE denote the eight directions North, South, East, West, North-East, South-West, North-West, and South-East. Hence, the explicit discretization scheme for equation (42) is given by
(53)ui,j,t+1=ui,j,t+18·Δt·gEα·∇Eui,j,t+gSEα·∇SEui,j,t+gSα·∇Sui,j,t+gSW·∇SWui,j,t+gWα·∇Wui,j,t+gNWα·∇NWui,j,t+gNα·∇Nui,j,t+gNEα·∇NEui,j,t,where
(54)gEα=g∇Eαui,j,gSEα=g∇SEαui,j,gSα=g∇Sαui,j,gSWα=g∇SWαui,j,gWα=g∇Wαui,j,gNWα=g∇NWαui,j,gNα=g∇Nαui,j,gNEα=g∇NEαui,j.


## 3. Experimental Results

In this section, experimental results are obtained applying the proposed fractional-order model on the512×512 images which were corrupted by the Gaussian noise with mean zero and variances v=0.01 and v=0.02. Gaussian noise is common in images acquired by cameras and telescopes, and it modifies all the pixels of the image. Furthermore, this noise type is one of the most studied in the literature on account of the frequency with which it occurs in practical situations. Our model is compared with the classical Perona-Malik (PM) anisotropic diffusion model. We set Δt=0.25 and used Δx=1 and Δy=1. To implement equation (53) numerically, in (47), we considered N=2, which means that in the approximation of fractional derivative, we take into account two points in each direction for experimental calculus. In the literature, there are some methodologies for the estimation of the contrast parameter K [61]. However, for reasons of simplification, we will use fixed contrast parameters which are K=1,2,3,4,5,6,7. The number of iterations used is 40. The performance of the two PDE filters has been assessed by using some well-known quality measures such as the PSNR (peak signal-to-noise ratio), MSE (mean square error), and MAE (mean absolute error) which are defined by [60]
(55)MSE=∑i=1N∑k=1Mxik−oik13×NM,MAE=∑i=1N∑k=1Mxik−oik223×NM,PSNR=20log10255MSE,where N and M denote the width and height of the image, respectively, xik is the pixel in the filtered image, and oik is the pixel in the original image. ∥·∥1 denotes L1 norm and ∥·∥2 denotes L2 norm (Euclidean distance). In the comparison of statistic parameters, it is important to note that the larger the PSNR value, the better the statistical result. On the contrary, the smaller the MSE and MAE, the better the result of the model.Figures2–7 show a comparison of denoised images obtained by our proposal and the Perona-Malik models. Tables 1–3 show the performance of our proposal compared to the Perona-Malik model using the quality measures MAE, MSE, and PSNR. All these tables have been obtained by using K=1,2,3,4,5,6,7. From Tables 1–3, it can be seen that the PSNR values of the proposed algorithm are higher than those of the Perona-Malik model; thus, it can efficiently delete the noise while preserving the original content of the image. Meanwhile, the values of MAE and MSE of the proposed algorithm are lower, which shows that the proposed method can better preserve texture details in the denoised image and has a better visual effects. Images of the third and fourth rows of Figures 2–7 have been obtained by using α=0.1 and α=0.2, respectively. After employing the proposed filter with the two previous mentioned α values combined with the above indicated threshold parameters K, we observed that the α parameter with the best statistic results is α=0.1 when K=5, as shown in Tables 1–3.Figure 2
Threshold parameterK=6 and iteration number=20. (a) Lena original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Figure 3
Threshold parameterK=7 and iteration number=20. (a) Lena original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Figure 4
Threshold parameterK=6 and iteration number=20. (a) Tiger original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Figure 5
Threshold parameterK=7 and iteration number=20. (a) Tiger original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Figure 6
Threshold parameterK=6 and iteration number=20. (a) Boat original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Figure 7
Threshold parameterK=7 and iteration number=20. (a) Boat original image in the left, noisy image with Gaussian noise with mean=0 and variance=0.01 in the middle, and noisy image with Gaussian noise with mean=0 and variance=0.02 in the right. (b) Processed by the PM scheme. (c) Processed by the proposed scheme for α=0.1. (d) Processed by the proposed scheme for α=0.2.
(a)(b)(c)(d)Table 1
Lena PSNR, MSE, and MAE.
Perona-Malik modelProposed model withα=0.1KPSNRMSEMAEKPSNRMSEMAE128.756086.592510.6559132.691534.98884.3134228.997081.91909.8173232.855733.69104.1568329.461773.60508.5643332.857033.68054.1559430.244661.46447.0600432.913933.24224.0919531.351647.63425.5919532.903733.32044.1114632.481936.71884.5154632.906633.29854.1086732.589635.81994.9641732.902433.33034.1126Table 2
Boat PSNR, MSE, and MAE.
Perona-Malik modelProposed model withα=0.1KPSNRMSEMAEKPSNRMSEMAE129.356975.40308.5066132.399537.42224.6925229.463573.57448.5799232.670435.15964.3854330.315660.46686.5120332.252038.71484.7724430.873253.18116.0197432.309338.20804.8410531.986841.15244.7247532.714134.80744.3401632.406137.36584.7866632.642835.38384.2737732.640535.40234.5418732.646035.35734.2702Table 3
Lion PSNR, MSE, and MAE.
Perona-Malik modelProposed model withα=0.1KPSNRMSEMAEKPSNRMSEMAE128.750486.703813.5416131.256148.69316.8764229.003281.801711.8001231.194849.38607.0452329.118579.658312.2037331.356647.57956.6767429.459873.642411.2282431.362447.51636.6653530.065864.04799.2298531.365147.48676.6602630.628656.26239.2298631.096850.51277.3495731.037351.20997.3767731.225049.04357.0367


## 4. Conclusion

In this paper, we proposed a new version of Perona-Malik diffusion (PMD) using the Caputo-Fabrizio fractional gradient inside the diffusivity function. Experiments showed that filtered images by the proposed method look better than results with the Perona-Malik filter. Our technique has demonstrated good performance in visual quality, higher peak signal-to-noise ratio (PSNR), and lower value of mean absolute error (MAE) and mean square error (MSE) than the Perona-Malik filter. Images of the third and fourth rows from Figures2–7 have been obtained by using α=0.1 and α=0.2, respectively. Experiments show that the best statistics are obtained for the fractional derivative of order α=0.1 when K=5, as can be seen in Tables 1–3. As future work, we want to introduce the convolution between the Caputo-Fabrizio fractional gradient and the Gaussian filter inside the edge-stopping function.



---

*Source: 7624829-2020-08-01.xml*