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

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 signalto-noise ratio (PSNR), and lower value of mean absolute error (MAE) and mean square error (MSE).


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, uðx, y, tÞ represents the image intensity values in the position ðx, yÞ ∈ Ω ⊂ ℝ 2 for a time t, ∂Ω is the smooth boundary, and u 0 ðx, yÞ is the original image. It is a classical result that for any bounded u ∈ Cðℝ 2 Þ, the linear diffusion process For more details, see [1]. However, k can be a smooth kernel or a blur kernel. In 1960, Gabor remarked that the dif-ference 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], then 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]: where uðx, y, 0Þ is the original image intensity function, u ðx, y, tÞ is the smoothed image intensity function in time t, ∇ denotes the gradient, div ð·Þ is the divergence operator, ∂ uðx, y, tÞ/∂n denotes the derivative in normal direction to the boundary, and gð|∇uðx, y, tÞ | Þ is the diffusion coefficient function, also called edge-stopping function. The edgestopping function gðsÞ is a nonnegative decreasing function satisfying two conditions. One is that gðsÞ = 1 as s → 0, so that the rate of diffusion is high within uniform or inner regions, and the other one is gðsÞ → 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 where s = |∇uðx, 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 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][9][10][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: 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 σ * uðx, y, tÞÞ is the regularized version of ð∇uðx, 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][14][15][16][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 where the fractional order is α (0 ≤ α ≤ 2) and uðx, y, tÞ is the smooth gray scale image at time t. The fractional-order gradient vector with α order is defined as 2 Abstract and Applied Analysis where ∇ α p uðx, y, tÞ represents the partial-order derivative of uðx, 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 As observed, equation (19) is obtained by substituting in equation (12), the integer spatial derivative for fractionalorder 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 wellposedness of the following fourth-order model for noise removal, using fractional derivatives defined by Fourier transform.
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][25][26], circuit theory [27][28][29], and image processing [30,31], among other applications [11,[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][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 ∇ α uðx, y, tÞ into the diffusivity function as follows: The flow function ϕ corresponding to equation (25) is a function depending on two variables and takes the following form: where g is given by (9) and s = ∇uðx, y, tÞ and v = ∇ α uðx, y, tÞ . This means that the graph of function (28) is a surface. For 3 Abstract and Applied Analysis simplicity reasons, the integer-order gradient ∇uðx, y, tÞ will be fixed by a real positive value while the fractional-order gradient ∇ α uðx, 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.
1.1. Some Basic Definitions in Fractional Calculus. Let a, b, p ∈ ℝ, such that a < b and 1 ≤ p ≤ ∞. We denote by C 0 ða, bÞ the space of all continuous functions on ½a, b with compact support. We denote C k 0 ða, bÞ as The Sobolev space W 1,p ða, bÞ is defined by W 1,p ða, bÞ = fu ∈ L p ða, bÞ;∃g ∈ L p ða, bÞ such that We set 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 with u ∈ H 1 ða, bÞ, b > a and α ∈ ð0, 1Þ. For more details, see [4,[54][55][56][57][58][59]. The remainder of this paper is organized as follows: Section 2 presents discretization of numerical schemes.

Abstract and Applied Analysis
Experimental results are provided in Section 3. Finally, conclusions are summarized in Section 4.

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.

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 interval ½0, 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 As in digital gray image uðx, 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 where An analogue expression is obtained for D α 0y u.

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 is given in [13], as where with 0 < Δt ≤ 1/4 for the numerical scheme to be stable.

Numerical
Method of the Proposed Model. As we said above, we propose to introduce the Caputo-Fabrizio fractional gradient ∇ α uðx, y, tÞ inside the diffusivity function of the Perona-Malik model where g ∇u j j ð Þ= is used in this work. Considering that ∇u j j= then the discretization of |∇u | , using the forward finite difference scheme, is given by The discrete fractional-order gradient ∇ α u is an 8dimensional vector, as in [18]: where T represents the transpose of the vector and each component is defined as 5 Abstract and Applied Analysis where C α,k is given by (34). The discretization of the usual gradient, following the idea in [3], in eight directions, is given by where ∇ N u, ∇ S u, ∇ E u, and ∇ W u are defined as in (37), (38), (39), and (40), respectively, and 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 where

Experimental Results
In this section, experimental results are obtained applying the proposed fractional-order model on the 512 × 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] where N and M denote the width and height of the image, respectively, x ik is the pixel in the filtered image, and o ik is the pixel in the original image. ∥·∥ 1 denotes L 1 norm and ∥· ∥ 2 denotes L 2 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.

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 Figures 2-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 edgestopping function.

Data Availability
To support this study, the USC-SIPI Image Database was used.

Conflicts of Interest
The authors declare that they have no conflict of interest.