# Image Denoising Based on Adaptive Fractional Order with Improved PM Model

Authors: Jimin Yu; Rumeng Zhai; Shangbo Zhou; LiJian Tan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In order to improve the image quality, in this paper, we propose an improved PM model. In the proposed model, we introduce two novel diffusion coefficients and a residual error term and replace the integer differential operator with the fractional differential operator in the PM model. The diffusion coefficients can be used effectively for edge detection and noise removal. The residual error term can help to prevent image distortion. Fractional order differential operator has a good characteristic that it can enhance image texture information while removing image noise. Additionally, in the two new diffusion coefficients, a novel method is proposed for automatically setting parameterk, and it does not need to do any experiments to get the value of k. For the computing fractional order diffusion coefficient, we employ the discrete Fourier transform, and an iterative scheme is carried out in the frequency domain. In the proposed model, not only is the integer differential operator replaced with the fractional differential operator, but also the order of the fractional differentiation is determined adaptively with the local variance. Comparing with some existing models, the experimental results show that the proposed algorithm can not only better suppress noise, but also better preserve edge and texture information. Moreover, the running time is greatly reduced.


---

## Body

## 1. Introduction

Image denoising is a critical task in image processing, and the denoising algorithm is widely researched in recent years. The anisotropic diffusion filter has received much attention since it was first proposed by Perona and Malik in 1990, which is called Perona–Malik model (PM model)‎ [1]. In the classical PM model, the gradient value is used in the direction of east, west, south, and north to distinguish variations which are caused by noise or edge in a corrupted image ‎[2]. The shortcomings of PM model are that the model is easy to lose contrasting information and texture information and produce staircase effects ‎[3]. To solve these shortcomings, some improved algorithms have been proposed. For image denoising, Ma and Nie ‎[4] proposed an edge fusion scheme based on anisotropic diffusion models. Zhao et al. proposed an improved PM model based on local entropy ‎[5]. H. Zhang and Y. Zhang ‎[6] proposed an adaptive diffusion coefficient scheme based on anisotropic diffusion. Prasath and Delhibabu ‎[7] proposed a fuzzy diffusion coefficient which takes into account local pixel variability for better denoising and selective smoothing of edges. Sun et al. ‎[8] have investigated the edges and details blurring issues. To avoid diffusion perpendicular to edge direction, Wang et al. ‎[9] introduced a modified PM model using directional Laplacian. Prasath and Vorotnikov ‎[10] suggested a weighted anisotropic diffusion to reduce blurring and staircasing effects. Xu et al. ‎[11] suggested an adaptive thresholding in PM diffusion coefficient to better handle the diffusion as time elapsed. A new diffusion coefficient has been proposed by Tebini et al. [12, 13] for better control of the diffusion process in regions containing edge. Wang et al. ‎[14] proposed new second- and fourth-order anisotropic equations for denoising. All of these methods mentioned above are gradient dependent where the gradient controls the diffusion process and therefore degrades texture and fine details.At present, fractional calculus has been widely applied in many fields of science and engineering. As a result, differential equations with arbitrary orders have been large investigated for different applications in physics, fluid mechanics, physiology, engineering, potential theory, elasticity, and so on [15–18]. Yang et al. ‎[19] proposed an adaptive image denoising model of anisotropic diffusion based on fractional derivative. Bai and Feng ‎[20] proposed a fractional order anisotropic diffusion for image denoising. Yuan and Liu ‎[21] proposed an anisotropic diffusion model based on a new diffusion coefficient and fractional order differential for image denoising. Bai and Feng ‎[20] proposed a fractional order anisotropic diffusion model for image denoising, and the model can not only excellently remove noise but also save the edge information. To automatically select the optional order of fractional differentiation, Che et al. ‎[22] bring forward a denoising model in which the complexity of the local image texture was reflected by local variance and the order of fractional differentiation was decided adaptively. The above-mentioned image processing algorithms based on fractional order partial differential equation have made improvement on keeping detailed image information, texture information, good visual effects, and image denoising.The PM model mainly has three disadvantages as follows: staircase effects exist, the texture information is lacked, and the diffusion coefficient thresholdk is set by carrying out many experiments. The optimal order of traditional fractional differentiation is often preset manually. In this paper, in order to solve these problems, by combining adaptive fractional differential operator with improved PM model, a novel model is proposed. Additionally, in the proposed method, according to ‎[1, 25], we proposed two new diffuse coefficients, and automatically setting diffusion coefficient parameter k is developed. Based on above measurement, a denoising algorithm based on adaptive fractional order with improved PM model is put forward. In addition, in the improved model, the local image texture complexity and the fractional order can be reflected by the local variance.The paper is organized as follows. In Section2, the fractional order definition and PM model are introduced. In Section 3, analysis of image denoising is shown based on adaptive fractional order with PM model. In Section 4, the experimental results and analysis are illustrated. Finally, conclusions are drawn in Section 5.


## 2. Introduction of Fractional Order and PM Model

### 2.1. PM Denoising Model

In the classical Perona–Malik equation ‎[1], proposed by Perona and Malik, the diffusion process is defined by a partial differential equation. The PM equation is defined as follows:(1)∂ux,y,t∂t=div⁡c∇u2·∇u,ux,y,0=u0,where div is the divergence operator,∇ is the gradient operator, · is the amplitude, c(∇u) is the diffusion coefficient equation, u0 is the original noise image, andt is the introduced time, used to exhibit diffusion duration. For the diffusion coefficient function c(∇u), Perona et al. gave two classic options as (2) and (3), respectively.(2)c∇u=11+∇u/k2,(3)c∇u=exp⁡-∇uk2,where k is diffusion coefficient threshold parameter, ∇u is the gradient of the image, and in PM model, the diffusion coefficient c(∇u) is a nonnegative function on the magnitude of local gradient information. Its value is generally inversely proportional to the gradient magnitude of the image. In the isotropic flat area, the gradient modulus is generally small, and the value of the diffusion coefficient is close to 1. The diffusion behavior of the model is close to the thermal diffusion pattern, which achieves the purpose of smoothing flat area of the image. In the edge region of the image, the gradient modulus is very large, and the diffusion coefficient is close to zero. The diffusion rate of the model is almost zero.


### 2.2. Introduction of Fractional Order

The commonly used fractional order differentiation is Grünwald-Letnikov (G-L) definition and Riemann-Liouville (R-L) definition [26, 27]. However, the definition in frequency domain is easier to implement and a simple definition, so we use the Fourier transform to compute the fractional derivative in this paper. According to the basic theory of signal processing, the Fourier transform of energy signal can be obtained; based on the Fourier transform, we can draw the amplitude frequency characteristic. Figure 1 shows the amplitude frequency characteristic curves of fractional order differentiation with different fractional order.Figure 1
Amplitude frequency characteristic curves of fractional order differential.Figure1 shows that, in the low-frequency field with 0<ω<1, corresponding to the smooth regions of image, as the order increases, fractional differential preserves the low-frequency (LF) portion of signal to a certain degree nonlinearly. In the high frequency section with ω>1, as the differential order and frequency value increase, the increase rate of the characteristic curves is shown in a nonlinear rapid growth. Therefore, by using fractional order differential, it can not only remove noise, but also preserve edge and detail information and avoid edge being oversmoothed.Considering the diffusion operator of PM model, it only diffuses along the orthogonal direction of image gradient, and there is not diffusion in the parallel direction; therefore the detail texture information will be lost. To effectively retain the detail texture information of a denoising image, the integer differential operator should be replaced with the fractional order differential operator.


## 2.1. PM Denoising Model

In the classical Perona–Malik equation ‎[1], proposed by Perona and Malik, the diffusion process is defined by a partial differential equation. The PM equation is defined as follows:(1)∂ux,y,t∂t=div⁡c∇u2·∇u,ux,y,0=u0,where div is the divergence operator,∇ is the gradient operator, · is the amplitude, c(∇u) is the diffusion coefficient equation, u0 is the original noise image, andt is the introduced time, used to exhibit diffusion duration. For the diffusion coefficient function c(∇u), Perona et al. gave two classic options as (2) and (3), respectively.(2)c∇u=11+∇u/k2,(3)c∇u=exp⁡-∇uk2,where k is diffusion coefficient threshold parameter, ∇u is the gradient of the image, and in PM model, the diffusion coefficient c(∇u) is a nonnegative function on the magnitude of local gradient information. Its value is generally inversely proportional to the gradient magnitude of the image. In the isotropic flat area, the gradient modulus is generally small, and the value of the diffusion coefficient is close to 1. The diffusion behavior of the model is close to the thermal diffusion pattern, which achieves the purpose of smoothing flat area of the image. In the edge region of the image, the gradient modulus is very large, and the diffusion coefficient is close to zero. The diffusion rate of the model is almost zero.


## 2.2. Introduction of Fractional Order

The commonly used fractional order differentiation is Grünwald-Letnikov (G-L) definition and Riemann-Liouville (R-L) definition [26, 27]. However, the definition in frequency domain is easier to implement and a simple definition, so we use the Fourier transform to compute the fractional derivative in this paper. According to the basic theory of signal processing, the Fourier transform of energy signal can be obtained; based on the Fourier transform, we can draw the amplitude frequency characteristic. Figure 1 shows the amplitude frequency characteristic curves of fractional order differentiation with different fractional order.Figure 1
Amplitude frequency characteristic curves of fractional order differential.Figure1 shows that, in the low-frequency field with 0<ω<1, corresponding to the smooth regions of image, as the order increases, fractional differential preserves the low-frequency (LF) portion of signal to a certain degree nonlinearly. In the high frequency section with ω>1, as the differential order and frequency value increase, the increase rate of the characteristic curves is shown in a nonlinear rapid growth. Therefore, by using fractional order differential, it can not only remove noise, but also preserve edge and detail information and avoid edge being oversmoothed.Considering the diffusion operator of PM model, it only diffuses along the orthogonal direction of image gradient, and there is not diffusion in the parallel direction; therefore the detail texture information will be lost. To effectively retain the detail texture information of a denoising image, the integer differential operator should be replaced with the fractional order differential operator.


## 3. Analysis of Improved PM Model

### 3.1. The Proposed Method

To solve the shortcomings that the classical PM model is tending to cause staircase effect and lose texture information, we propose an image denoising algorithm based on adaptive fractional order of PM model:(4)∂u∂t=div⁡cDαu2Dαu+λu-u0ux,y,0=u0x,y,where α is the fractional order, Dαu is the fractional derivative operator, Dαu=(Dαxu,Dαyu) and Dαu=Dαx2u+Dαy2u, c(Dαu2) is the diffusion coefficient equation, λ(u-u0) is residual error, and λ is an adjustment parameter, which can be employed to control the fidelity between the noisy image and the denoising image to preserve main characteristics of the noise image.The diffusion coefficientc(·) is very important in the PM model; in ‎[23], the diffusion coefficient c(·) is proposed based on adaptive fractional order anisotropic diffusion; however, the equality of image denoising is affected by noise, and the running time is long, so the efficiency of the algorithm greatly reduces. In ‎[23], the diffusion coefficient is given as (5)cDαu=11+Dαxu2+Dαyu2.The typical disadvantages of the Perona–Malik model are that it can not effectively protect sharp edges and texture details during denoising, but this drawback can be overcame by using the appropriate diffusion coefficient. In order to better preserve detailed and texture information and greatly reduce the running time of the algorithm to improve the operational efficiency, according to [1, 25], we proposed two new diffuse coefficients as (6) and (7). Additionally, in the two new diffusion coefficients, a novel method is proposed for automatically setting parameterk; it does not need to do experiments to getk.(6)c1Dαu=11+Dαxu2+Dαyu2/K2,(7)c2Dαu=11+Dαxu2+Dαyu2/kβDαxu2+Dαyu2,where β((Dαxu2+Dαyu2))=2-2/(1+(Dαxu2+Dαyu2)/k2).The diffuse coefficient parameterk is used to adjust the gradient. By using the adaptive factork, it can expect that the improved PM equation has a better effect on noise removal, edge detection, and image enhancement. C1 is the classical diffusion coefficient. C2 is an adaptive operator based on the variable exponent; it mainly used the edge indicator to segment the image into inside regions and the region’s nearby boundaries, and the diffusion mode is adaptively determined by the edge indicator; C2 as diffuse coefficient not only removes noise in the inside regions, but also can preserve edge in the region’s nearby boundaries. From above analysis, by employing these two diffusion equations to improve the PM model, it will sufficiently preserve detail information and remove noise.The relationship curves between the different diffuse coefficient and the gradient magnitude are shown as in Figure1, in which we can visually see the enhancement and weakening of the diffusion effect controlled by the diffuse coefficient.From Figure2, it is obvious that, in the low magnitude part of 0<grad(u)<30, corresponding to the image of the smooth region, the diffusion rate of the model is almost 1; thus, the noise can be efficiently removed. In the edge region of the image, the gradient magnitude is very large, and the diffusion coefficient is close to zero; as the diffusion rate of the model is almost zero, the edge of the image can be better protected.Figure 2
The relationship curves between the different edge stopping function and the gradient magnitude (k=5).


### 3.2. Numerical Implementation of Proposed Denoising Model

In this subsection, in order to analyzeu in (4), the energy functional is proposed, and an anisotropic diffusion is then shown to be an energy-dissipating process which can seek the minimum of the energy functional.The energy functional can be defined as (8) in the continuous images space over the support Ω.(8)Eu=∫ΩfDαudΩ-λ2u-u02,where Ω is the image domain,f′s=c(s)s, c(s) is the diffuse coefficient, and f(·)≥0 is an increasing function. For solving minimum value problem of formula (8), the Euler-Lagrange equation can be employed. For any test function η∈C∞(Ω), we define (9)Φa=∫ΩfDαu+aηdxdy-λ2u-aη-u02.When a is equal to zero, for all η∈C∞(Ω), we can obtain(10)Φ′0=∫ΩcDαu2DαxuDαxη+cDαu2DαyuDαxη+λu-u0×ηdxdy=∫ΩDαx∗cDαu2Dαxu+Dαy∗cDαu2Dαyu+λu-u0×ηdxdy,where Dαx∗ is the adjoint matrix of Dαx, Dαy∗ is the adjoint matrix of Dαy, and λ is an adjustment parameter. Therefore, when the Φa satisfies the condition Φ′(0)=0, for arbitrariness of η, we can obtain the following equation:(11)Dαx∗cDαu2Dαxu+Dαy∗cDαu2Dαyu+λu-u0=0.By introducing time variable, (11) may be solved with the following gradient descent procedure: (12)∂u∂t=-Dαx∗cDαu2Dαxu-Dαy∗cDαu2Dαyu-λ×u-u0.The solution of (11) can be obtained by the limit of the solution of (12) as t→∞. We maybe stop the time evolution earlier to achieve an optimal tradeoff between noise removal and edgepreservation.We assume that the size of the original discrete imageu is m×m pixels. Then we use the 2D discrete Fourier transform (2D DFT) to compute the fractional order difference. And the discrete Fourier transform is easy to implement and has low computational cost. Thus, the fractional order differential equation (12) can be easy calculated using discrete Fourier transform (2D DFT).The definition of 2D DFT is as follows:(13)u^ω1,ω2=1m∑x,y=0m-1ux,yexp⁡-j2πω1x+ω2ym.The translation property of the 2D DFT is (14)ux-x0,y-y0⟷exp⁡-j2πω1x0+ω2y0mu^ω1,ω2.Thus, the definition of the first-order partial differential can be obtained as (15)ux,y-ux-Δx,y⟷1-exp⁡-j2πω1Δxmu^ω1,ω2.Therefore, when the smallest equal interval of u(x,y) is chosen as Δx=Δy=1, the first-order partial differential is Dxu=u(x,y)-u(x-1,y).Based on formula (15), the fractional order partial differential for x in the frequency domain can be defined as(16)Dαxu⟷1-exp⁡-j2πω1mαu^ω1,ω2,where α denotes the adaptive fractional order and fory, Dαyu can also be obtained in the same way as computing Dαxu.In this paper, we introduce the central difference method to calculate the fractional order differential. And the trigonometric interpolation can be automatically implemented by the translation property of the 2D DFT; thus, we can obtain(17)ux+α2,y⟷exp⁡jπαω1mu^ω1,ω2.In formula (17), α/2 is not only an integer; therefore, the continuous image u is clearly defined for all real number x, and we can obtain the central difference equation:(18)D~αxu=Dαxux+α2,y⟷1-exp⁡-j2πω1mα×exp⁡jπαω1mu^ω1,ω2.In addition, the inverse two-dimensional Fourier transform is introduced; the symbolF denotes the two-dimensional discrete Fourier transform operator; contrarily, F-1 represents the inverse operator, so we can obtain(19)D~αxu=F-11-exp⁡-j2πω1mα×exp⁡jπαω1mFu.Next we mainly calculate the adjoint matrixD~∗αxu of D~αxu; we define K1 as the pure diagonal matrix in frequency domain as(20)K1=diag⁡1-exp⁡-j2πω1mα×exp⁡jπαω1m.We get(21)D~∗αx=F-1∘K1∘F∗=K1∘F∗F-1∗=F∗∘K1∗∘F-1∗.Because the conjugate matrix of F is F-1and K1 is a pure diagonal matrix, K1∗ denotes the complex conjugate matrix of K1 and conj· is the complex conjugate operator, so (22)D~αx∗u=F-1∘K1∗∘F=F-1conj1-exp⁡-j2πω1mα×exp⁡jπαω1mFu.According to the algorithm,D~αy and D~αy∗ can also be obtained by using the same algorithm getting D~αx and D~αx∗, so the denoising algorithm can be obtained as follows:(23)g^=K1∗∘FcD~αu2D~αxu+K2∗∘FcD~αu2D~αyu+λ×u-u0,where D~αu=(D~αxu,D~αyu), K1∗=diag⁡conj1-exp⁡-j2πω1/mα×exp⁡jπαω1/m, and K2∗ = diag⁡conj1-exp⁡-j2πω2/mα×exp⁡jπαω2/m. The diffuse coefficient c(•) is the two new diffuse coefficients as. (6) and (7) proposed by us.In this paper, there are two important factors: the gradient parameterk and the fractional order α. For the first factork, the diffuse coefficient will increase as k increases, as increasing of diffuse coefficient will result in the loss of more image texture details. The parameter k that does not meet the requirements will significantly affect the denoising performance of the model. Thus, choosing the appropriate value is important for the simulation results. In most models,kpreset a constant, which cannot reflect the change of the image gradient accurately. Here, a new method for automatically setting k is proposed. The diffusion coefficient parameter k can be computed by (24)k=k0e-1/6nΔt,where k0=mean(Du0).The mean denotes the average operator,nis the iteration number, Δtis the time step, and u0 is noise image.For the second factor, the fractional orderα also plays an important role in the proposed model. In this paper, we propose an adaptive fractional order method. Considering the exponential function is an increasing function which is consistent with the change of the order; thus, the adaptive fractional order can be chosen as follows:(25)α=efσ+k2,where k2 is adjusting parameter and in the experiment we take k2=0.5 by carrying out a large number of experiments; according to the form of the local entropy function proposed by ‎[5], f(σ) can be calculated by following equation:(26)fσ=k0+k1×σx,y-min⁡σx,ymax⁡σx,y-min⁡σx,y,(27)μx,y=1m1×m2∑x,y∈Mux,y,σx,y=1m1×m2∑x,y∈Mu2x,y-μ2x,y,where μx,y denotes the local mean of each pixel in the image ux,y, σx,y denotes local variance of each pixel, max is the maximum operator, min is the minimum operator, k0=0, k1=0.693, and M is m1×m2 rectangular window of which x,y is taken as the center.Figure3 shows the denoising results with (6) as the diffusion coefficient. Figure 4 shows the comparison results of denoising images, where the diffusion coefficient is taken, respectively, as ‎[23] and the proposed equation (7) (σ=15). Table 1 lists the comparison results of PSNR, SSIM, and information entropy (ENTROY).Table 1
The comparison results of PSNR, SSIM, and information entropy (ENTROY) between [23] and (6), the Gaussian white noise with mean square deviation σ=15.
Image PSNR SSIM Running time(s) ENTROY Reference [23] Equation (6) Reference [23] Equation (6) Reference [23] Equation (6) Reference [23] Equation (6) Jinxiang 33.8955 35.0832 0.8712 0.9112 241.24 17.218 6.4016 6.4370 Satellite 29.8487 31.4236 0.8193 0.8573 243.60 13.475 6.8716 6.8849Figure 3
Denoising images with formula (6) as the diffusion coefficient.
(a)
Original (b)
Noise (σ=15) (c)
Denoised (σ=15) (d)
Noise (σ=30) (e)
Denoised (σ=30)Figure 4
Comparison results of denoising images between methods of [23] and the proposed formula (7).
(a)
Reference [23] (b)
Proposed equation (6) (c)
Reference [23] (d)
Proposed equation (6)From Figures3 and 4, it can be seen that, with formula (6) and formula (7) as the diffusion coefficient, the denoised images have a better visual effect. Therefore both of two diffusion coefficients can effectively remove noise. Additionally, from the data of Table 1, the values of PSNR, SSIM, and information entropy (ENTROY) of denoised image with formula (6) as diffuse coefficient are higher than that of the algorithm of ‎[23]. Furthermore, the running time is significantly shorter. Therefore, the performance of the proposed algorithm is better than that of the algorithm of ‎[23].


### 3.3. Numerical Implementation

The implementation of the proposed method consists of the following steps:(1) Setn=1, un(x,y)=u(x,y), Δt=0.1, k0=0, k1=0.693, and k2=0.8.(2) Compute the local mean and variance of the original imageu0(x,y), and then compute the order of the fractional differential operator using formula (27).(3) According to formula (24), calculate the diffusion coefficient parameterk.(4) Compute the 2D DFTu^n(x,y) of un(x,y), and according to the formula (19) calculate the partial differential D~αxun(x,y) and D~αyun(x,y), then use formula (6) and (7) as diffusion coefficient to compute the overall diffusion equation of formula (23).(5) Computeu^n+1(x,y) with the iterative format to get u^n+1(x,y)=u^n(x,y)-Δt×g^n(x,y).(6) Letn=n+1, u0(x,y)=u^n+1(x,y), return to step (2) until the peak signal-to-noise ratio (PSNR) reaches the maximum.


## 3.1. The Proposed Method

To solve the shortcomings that the classical PM model is tending to cause staircase effect and lose texture information, we propose an image denoising algorithm based on adaptive fractional order of PM model:(4)∂u∂t=div⁡cDαu2Dαu+λu-u0ux,y,0=u0x,y,where α is the fractional order, Dαu is the fractional derivative operator, Dαu=(Dαxu,Dαyu) and Dαu=Dαx2u+Dαy2u, c(Dαu2) is the diffusion coefficient equation, λ(u-u0) is residual error, and λ is an adjustment parameter, which can be employed to control the fidelity between the noisy image and the denoising image to preserve main characteristics of the noise image.The diffusion coefficientc(·) is very important in the PM model; in ‎[23], the diffusion coefficient c(·) is proposed based on adaptive fractional order anisotropic diffusion; however, the equality of image denoising is affected by noise, and the running time is long, so the efficiency of the algorithm greatly reduces. In ‎[23], the diffusion coefficient is given as (5)cDαu=11+Dαxu2+Dαyu2.The typical disadvantages of the Perona–Malik model are that it can not effectively protect sharp edges and texture details during denoising, but this drawback can be overcame by using the appropriate diffusion coefficient. In order to better preserve detailed and texture information and greatly reduce the running time of the algorithm to improve the operational efficiency, according to [1, 25], we proposed two new diffuse coefficients as (6) and (7). Additionally, in the two new diffusion coefficients, a novel method is proposed for automatically setting parameterk; it does not need to do experiments to getk.(6)c1Dαu=11+Dαxu2+Dαyu2/K2,(7)c2Dαu=11+Dαxu2+Dαyu2/kβDαxu2+Dαyu2,where β((Dαxu2+Dαyu2))=2-2/(1+(Dαxu2+Dαyu2)/k2).The diffuse coefficient parameterk is used to adjust the gradient. By using the adaptive factork, it can expect that the improved PM equation has a better effect on noise removal, edge detection, and image enhancement. C1 is the classical diffusion coefficient. C2 is an adaptive operator based on the variable exponent; it mainly used the edge indicator to segment the image into inside regions and the region’s nearby boundaries, and the diffusion mode is adaptively determined by the edge indicator; C2 as diffuse coefficient not only removes noise in the inside regions, but also can preserve edge in the region’s nearby boundaries. From above analysis, by employing these two diffusion equations to improve the PM model, it will sufficiently preserve detail information and remove noise.The relationship curves between the different diffuse coefficient and the gradient magnitude are shown as in Figure1, in which we can visually see the enhancement and weakening of the diffusion effect controlled by the diffuse coefficient.From Figure2, it is obvious that, in the low magnitude part of 0<grad(u)<30, corresponding to the image of the smooth region, the diffusion rate of the model is almost 1; thus, the noise can be efficiently removed. In the edge region of the image, the gradient magnitude is very large, and the diffusion coefficient is close to zero; as the diffusion rate of the model is almost zero, the edge of the image can be better protected.Figure 2
The relationship curves between the different edge stopping function and the gradient magnitude (k=5).


## 3.2. Numerical Implementation of Proposed Denoising Model

In this subsection, in order to analyzeu in (4), the energy functional is proposed, and an anisotropic diffusion is then shown to be an energy-dissipating process which can seek the minimum of the energy functional.The energy functional can be defined as (8) in the continuous images space over the support Ω.(8)Eu=∫ΩfDαudΩ-λ2u-u02,where Ω is the image domain,f′s=c(s)s, c(s) is the diffuse coefficient, and f(·)≥0 is an increasing function. For solving minimum value problem of formula (8), the Euler-Lagrange equation can be employed. For any test function η∈C∞(Ω), we define (9)Φa=∫ΩfDαu+aηdxdy-λ2u-aη-u02.When a is equal to zero, for all η∈C∞(Ω), we can obtain(10)Φ′0=∫ΩcDαu2DαxuDαxη+cDαu2DαyuDαxη+λu-u0×ηdxdy=∫ΩDαx∗cDαu2Dαxu+Dαy∗cDαu2Dαyu+λu-u0×ηdxdy,where Dαx∗ is the adjoint matrix of Dαx, Dαy∗ is the adjoint matrix of Dαy, and λ is an adjustment parameter. Therefore, when the Φa satisfies the condition Φ′(0)=0, for arbitrariness of η, we can obtain the following equation:(11)Dαx∗cDαu2Dαxu+Dαy∗cDαu2Dαyu+λu-u0=0.By introducing time variable, (11) may be solved with the following gradient descent procedure: (12)∂u∂t=-Dαx∗cDαu2Dαxu-Dαy∗cDαu2Dαyu-λ×u-u0.The solution of (11) can be obtained by the limit of the solution of (12) as t→∞. We maybe stop the time evolution earlier to achieve an optimal tradeoff between noise removal and edgepreservation.We assume that the size of the original discrete imageu is m×m pixels. Then we use the 2D discrete Fourier transform (2D DFT) to compute the fractional order difference. And the discrete Fourier transform is easy to implement and has low computational cost. Thus, the fractional order differential equation (12) can be easy calculated using discrete Fourier transform (2D DFT).The definition of 2D DFT is as follows:(13)u^ω1,ω2=1m∑x,y=0m-1ux,yexp⁡-j2πω1x+ω2ym.The translation property of the 2D DFT is (14)ux-x0,y-y0⟷exp⁡-j2πω1x0+ω2y0mu^ω1,ω2.Thus, the definition of the first-order partial differential can be obtained as (15)ux,y-ux-Δx,y⟷1-exp⁡-j2πω1Δxmu^ω1,ω2.Therefore, when the smallest equal interval of u(x,y) is chosen as Δx=Δy=1, the first-order partial differential is Dxu=u(x,y)-u(x-1,y).Based on formula (15), the fractional order partial differential for x in the frequency domain can be defined as(16)Dαxu⟷1-exp⁡-j2πω1mαu^ω1,ω2,where α denotes the adaptive fractional order and fory, Dαyu can also be obtained in the same way as computing Dαxu.In this paper, we introduce the central difference method to calculate the fractional order differential. And the trigonometric interpolation can be automatically implemented by the translation property of the 2D DFT; thus, we can obtain(17)ux+α2,y⟷exp⁡jπαω1mu^ω1,ω2.In formula (17), α/2 is not only an integer; therefore, the continuous image u is clearly defined for all real number x, and we can obtain the central difference equation:(18)D~αxu=Dαxux+α2,y⟷1-exp⁡-j2πω1mα×exp⁡jπαω1mu^ω1,ω2.In addition, the inverse two-dimensional Fourier transform is introduced; the symbolF denotes the two-dimensional discrete Fourier transform operator; contrarily, F-1 represents the inverse operator, so we can obtain(19)D~αxu=F-11-exp⁡-j2πω1mα×exp⁡jπαω1mFu.Next we mainly calculate the adjoint matrixD~∗αxu of D~αxu; we define K1 as the pure diagonal matrix in frequency domain as(20)K1=diag⁡1-exp⁡-j2πω1mα×exp⁡jπαω1m.We get(21)D~∗αx=F-1∘K1∘F∗=K1∘F∗F-1∗=F∗∘K1∗∘F-1∗.Because the conjugate matrix of F is F-1and K1 is a pure diagonal matrix, K1∗ denotes the complex conjugate matrix of K1 and conj· is the complex conjugate operator, so (22)D~αx∗u=F-1∘K1∗∘F=F-1conj1-exp⁡-j2πω1mα×exp⁡jπαω1mFu.According to the algorithm,D~αy and D~αy∗ can also be obtained by using the same algorithm getting D~αx and D~αx∗, so the denoising algorithm can be obtained as follows:(23)g^=K1∗∘FcD~αu2D~αxu+K2∗∘FcD~αu2D~αyu+λ×u-u0,where D~αu=(D~αxu,D~αyu), K1∗=diag⁡conj1-exp⁡-j2πω1/mα×exp⁡jπαω1/m, and K2∗ = diag⁡conj1-exp⁡-j2πω2/mα×exp⁡jπαω2/m. The diffuse coefficient c(•) is the two new diffuse coefficients as. (6) and (7) proposed by us.In this paper, there are two important factors: the gradient parameterk and the fractional order α. For the first factork, the diffuse coefficient will increase as k increases, as increasing of diffuse coefficient will result in the loss of more image texture details. The parameter k that does not meet the requirements will significantly affect the denoising performance of the model. Thus, choosing the appropriate value is important for the simulation results. In most models,kpreset a constant, which cannot reflect the change of the image gradient accurately. Here, a new method for automatically setting k is proposed. The diffusion coefficient parameter k can be computed by (24)k=k0e-1/6nΔt,where k0=mean(Du0).The mean denotes the average operator,nis the iteration number, Δtis the time step, and u0 is noise image.For the second factor, the fractional orderα also plays an important role in the proposed model. In this paper, we propose an adaptive fractional order method. Considering the exponential function is an increasing function which is consistent with the change of the order; thus, the adaptive fractional order can be chosen as follows:(25)α=efσ+k2,where k2 is adjusting parameter and in the experiment we take k2=0.5 by carrying out a large number of experiments; according to the form of the local entropy function proposed by ‎[5], f(σ) can be calculated by following equation:(26)fσ=k0+k1×σx,y-min⁡σx,ymax⁡σx,y-min⁡σx,y,(27)μx,y=1m1×m2∑x,y∈Mux,y,σx,y=1m1×m2∑x,y∈Mu2x,y-μ2x,y,where μx,y denotes the local mean of each pixel in the image ux,y, σx,y denotes local variance of each pixel, max is the maximum operator, min is the minimum operator, k0=0, k1=0.693, and M is m1×m2 rectangular window of which x,y is taken as the center.Figure3 shows the denoising results with (6) as the diffusion coefficient. Figure 4 shows the comparison results of denoising images, where the diffusion coefficient is taken, respectively, as ‎[23] and the proposed equation (7) (σ=15). Table 1 lists the comparison results of PSNR, SSIM, and information entropy (ENTROY).Table 1
The comparison results of PSNR, SSIM, and information entropy (ENTROY) between [23] and (6), the Gaussian white noise with mean square deviation σ=15.
Image PSNR SSIM Running time(s) ENTROY Reference [23] Equation (6) Reference [23] Equation (6) Reference [23] Equation (6) Reference [23] Equation (6) Jinxiang 33.8955 35.0832 0.8712 0.9112 241.24 17.218 6.4016 6.4370 Satellite 29.8487 31.4236 0.8193 0.8573 243.60 13.475 6.8716 6.8849Figure 3
Denoising images with formula (6) as the diffusion coefficient.
(a)
Original (b)
Noise (σ=15) (c)
Denoised (σ=15) (d)
Noise (σ=30) (e)
Denoised (σ=30)Figure 4
Comparison results of denoising images between methods of [23] and the proposed formula (7).
(a)
Reference [23] (b)
Proposed equation (6) (c)
Reference [23] (d)
Proposed equation (6)From Figures3 and 4, it can be seen that, with formula (6) and formula (7) as the diffusion coefficient, the denoised images have a better visual effect. Therefore both of two diffusion coefficients can effectively remove noise. Additionally, from the data of Table 1, the values of PSNR, SSIM, and information entropy (ENTROY) of denoised image with formula (6) as diffuse coefficient are higher than that of the algorithm of ‎[23]. Furthermore, the running time is significantly shorter. Therefore, the performance of the proposed algorithm is better than that of the algorithm of ‎[23].


## 3.3. Numerical Implementation

The implementation of the proposed method consists of the following steps:(1) Setn=1, un(x,y)=u(x,y), Δt=0.1, k0=0, k1=0.693, and k2=0.8.(2) Compute the local mean and variance of the original imageu0(x,y), and then compute the order of the fractional differential operator using formula (27).(3) According to formula (24), calculate the diffusion coefficient parameterk.(4) Compute the 2D DFTu^n(x,y) of un(x,y), and according to the formula (19) calculate the partial differential D~αxun(x,y) and D~αyun(x,y), then use formula (6) and (7) as diffusion coefficient to compute the overall diffusion equation of formula (23).(5) Computeu^n+1(x,y) with the iterative format to get u^n+1(x,y)=u^n(x,y)-Δt×g^n(x,y).(6) Letn=n+1, u0(x,y)=u^n+1(x,y), return to step (2) until the peak signal-to-noise ratio (PSNR) reaches the maximum.


## 4. Experimental Results and Analysis

In this section, the performance of the proposed method will be compared with that of some algorithms, the selected three standard gray scale images include “barbara,” ”Jinxiang,” “Lenna,” and “Satellite” with size (512×512). Figures 5, 6, 7, and 8, respectively, show the comparison result of denoised images between the proposed model and CSF7 Filter, pseudo-PDEs, Curvature Filter, FDC Model, Adaptive TV Mode, and ZHAO algorithm, and Tables 2, 3, and 4, respectively, show the comparison results of the denoised images of PSNR, SSIM, ENTROY between proposed model and CSF7 Filter, pseudo-PDEs, Curvature Filter, FDC Model, Adaptive TV Mode, and ZHAO algorithm, the Gaussian white noise with mean square deviation δ=15. All the experiments are implemented with MatlabR2012a for Windows 10.Table 2
Comparison results of the PSNR of the denoising image with different algorithms.
Image CSF7 Filter ‎[24] FDC Model Pseudo-PDEs Curvature Filter ZHAO Model Model in this paper Barbara 28.83 28.6852 28.9968 27.4634 27.8251 28.913 Jinxiang 34.17 31.5397 33.063 33.261 31.6419 35.083 Lena 32.25 30.9927 31.9991 32.172 30.6552 32.267 Satellite 31.18 29.6473 30.1351 29.9147 28.6496 31.423Table 3
Comparison results of the SSIM of the denoising image with different algorithms.
Image CSF7 Filter ‎[24] FDC Model Pseudo-PDEs Curvature Filter ZHAO Model Model in this paper Barbara 0.8207 0.6883 0.8522 0.8073 0.8086 0.8302 Jinxiang 0.9091 0.7762 0.8570 0.8651 0.8274 0.9112 Lena 0.9085 0.8322 0.8569 0.8582 0.8319 0.8549 Satellite 0.8445 0.7911 0.8082 0.7946 0.7559 0.8573Table 4
Comparison results of the ENTROY of the denoising image with different algorithms.
Image CSF7 Filter ‎[24] FDC Model Pseudo-PDEs Curvature Filter ZHAO Model Model in this paper Barbara 7.4461 7.4228 7.4728 7.4308 7.5315 7.4913 Jinxiang 6.3432 6.3737 6.3829 6.2256 6.3038 6.4370 Lena 7.4094 7.4368 7.4342 7.3757 7.4382 7.4951 Satellite 6.8675 6.6550 6.8632 6.7904 6.8498 6.8849Figure 5
Barbara image denoising with different algorithms.
(a)
Original image (b)
Noise image (c)
CSF7 Filter (d)
Adaptive TV Model (e)
Pseudo_PDE (f)
Curvature Filter (g)
ZhAO Model (h)
FDC Model (i)
Proposed equation (6) (j)
Proposed equation (7)Figure 6
Jinxiang image denoising with different algorithms.
(a)
Original image (b)
Noise image (c)
CSF7 Filtering (d)
Adaptive TV Model (e)
Pseudo_PDE (f)
Curvature Filter (g)
ZhAO Model (h)
FDC Model (i)
Proposed equation (6) (j)
Proposed equation (7)Figure 7
Lena image denoising with different algorithms.
(a)
Original image (b)
Noise image (c)
CSF7 Filter (d)
Adaptive TV Model (e)
Pseudo_PDE (f)
Curvature Filter (g)
ZhAO Model (h)
FDC Model (i)
Proposed equation (6) (j)
Proposed equation (7)Figure 8
Satellite image denoising with different algorithms.
(a)
Original image (b)
Noise image (c)
CSF7 Filter (d)
Adaptive TV model (e)
Pseudo_PDE (f)
Curvature Filter (g)
ZhAO Model (h)
FDC Model (i)
Proposed equation (6) (j)
Proposed equation (7)And all are on a PC with Intel (R) Core (TM) 2 Duo CPU and 4 GB of RAM. In all experiments, the noisy images can be represented as follows:(28)I=I0+σ∗randn,where I0 is an original image andIis a noisy image. The “randn” represents the normal random noise with mean zero and variance 1. σ is the mean standard deviation of noise. To eliminate the influence of isolate point noise, mean filtering is used to presmooth noisy images.In addition, we will use quantitative evaluation index to objectively evaluate the performance of the proposed denoising algorithm; for an image denoising algorithm, its denoising effect can be evaluated with the ratio of peak signal-to-noise ratio (PSNR). The higher the value of PSNR is, the better the denoising effect of the algorithm is. PSNR is an engineering term that represents the ratio between the maximum possible power of a signal and the destructive noise power. The peak signal-to-noise ratio (PSNR) is defined as follows:(29)MSE=1M×N∑i=1M∑j=1Nui,j-u0i,j2,PSNR=10lg⁡2552MSE,where the size of the original image is M×N, u0(i,j) represents the original image, and u(i,j) represents the denoising image.Structural similarity index measurement (SSIM) is always used to consider the similarity between two images. The higher the value is, the closer the image content is. Structural similarity index measurement is defined as follows:(30)SSIM=2μaμb+c12σab+c2μa2+μb2+c1σa2+σb2+c2,where a and b are different images, μa and μb are the mean of a and b, respectively. σa2 and σb2 denote the variance of a and b, respectively. σab represents the covariance of a and b. C1 and C2 denote constants used to maintain stability. The comparison results of PSNR, SSIM, and ENTORY and running time of five different denoising algorithms are shown as in Tables 2, 3, 4, and 5.Figures5–8 make a visual comparison between the CSF filter, pseudo-PDEs, Curvature Filter, FDC Model, Adaptive TV Model, and ZHAO Model and the proposed denoising algorithm. From Figures 5–8, it can be seen that CSF7 Filter, pseudo-PDEs, Adaptive TV Model, Curvature Filter, FDC Model, and ZHAO Model can effectively remove the image noises; however, such as FDC Model, the denoised Lena image is blurred so that it lost a lot of texture details information. But the proposed algorithm has better visual effect and also can weaken the staircase effect. Moreover, it can better preserve edges and retain more texture information.In order to objectively evaluate the proposed algorithm, it needs to compare the peak signal-to-noise ratio and the information entropy and the similarity to compare the image edge enhancement and image weak texture retention. From Tables 2–4, it can be seen that the PSNR, SSIM, and ENTORY values of proposed algorithm are higher than that of the CSF7 Filter, pseudo-PDEs, Curvature Filter, FDC Model, and ZHAO Model; thus, it can efficiently remove the noise while retaining the original content of the image. Meanwhile, the value of SSIM of the adaptive fractional order algorithm proposed in this paper is slightly higher, which can better preserve texture details in the denoised image and has a good visual effect. From Table5, the running time of the proposed algorithm is obviously shorter than that of the other algorithms. In ‎[23], where formula (5) acting as the diffusion coefficient leads to the running time to be longer, it is greatly reducing the efficiency of the algorithm. To overcome this problem, we use a new diffusion coefficient instead of formula (5), we use formula (6) and (7) as new diffusion coefficients, and the parameter k of new diffusion coefficients is automatically set; it does not need to do any experiments to get k value. From the data of Table 1, with formula (6) as the diffusion coefficient obviously has a better denoising effect than with formula (5), the PSNR,SSIM, and ENTORY values of the proposed algorithm also are large than those of the other algorithms, and the running time obviously shortened.Table 5
Comparison results of the running time of the denoising image with different algorithms (unit: second).
Image CSF7 Filter‎ [24] FDC Model Pseudo-PDEs Curvature Filter ZHAO Model Model in this paper Barbara 8.63 30.6762 182.2906 37.1060 30.6762 12.5917 Jinxiang 8.79 809.6865 174.4096 45.1872 809.6865 17.2187 Lena 8.76 1256.3321 191.5609 14.1970 256.3321 15.8083 Satellite 8.53 184.6692 196.6252 13.3099 184.6692 13.4750


## 5. Conclusion

In this paper, for image denoising, we propose an adaptive fractional order differential PM model. The PM model integer differential operator is replaced with the fractional differential operator for retaining the texture detail of the image, so it can weaken the staircase effect and preserve detail information of image. In addition, the fractional order is adaptively determined by the local variance, it can further improve the ability of denoising. Moreover, the adjustment parameterk of the diffusion coefficient is automatically determined. The proposed algorithm not only can shorten the running time, but also can enhance the denoising capability. From experimental results, it can be seen that the proposed image denoising algorithm can effectively preserve edges and avoid edges being oversmoothed. In addition, the proposed algorithm can better preserve detail and texture information, and it has a good visual effect; furthermore, the PSNR and SSIM value are also higher than that of the other algorithms.



---

*Source: 9620754-2018-05-20.xml*