# Image Structure-Preserving Denoising Based on Difference Curvature Driven Fractional Nonlinear Diffusion

Authors: Xuehui Yin; Shangbo Zhou
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The traditional integer-order partial differential equations and gradient regularization based image denoising techniques often suffer from staircase effect, speckle artifacts, and the loss of image contrast and texture details. To address these issues, in this paper, a difference curvature driven fractional anisotropic diffusion for image noise removal is presented, which uses two new techniques, fractional calculus and difference curvature, to describe the intensity variations in images. The fractional-order derivatives information of an image can deal well with the textures of the image and achieve a good tradeoff between eliminating speckle artifacts and restraining staircase effect. The difference curvature constructed by the second order derivatives along the direction of gradient of an image and perpendicular to the gradient can effectively distinguish between ramps and edges. Fourier transform technique is also proposed to compute the fractional-order derivative. Experimental results demonstrate that the proposed denoising model can avoid speckle artifacts and staircase effect and preserve important features such as curvy edges, straight edges, ramps, corners, and textures. They are obviously superior to those of traditional integral based methods. The experimental results also reveal that our proposed model yields a good visual effect and better values of MSSIM and PSNR.


---

## Body

## 1. Introduction

A challenging task in designing image denoising model is to retain significant features, such as edges and texture details, while eliminating noise. Among a variety of integer-order partial differential equations based image denoising techniques [1–3], energy minimization based variational technique in two well-known forms of the nonlinear diffusion with anisotropic diffusion proposed by Perona and Malik [4] (Perona-Malik or PM) and total variation (TV) regularization proposed by Rudin et al. [5] (Rudin-Osher-Fatemi or ROF) are the ones that have shown a good edge preservation capability. Most of these nonlinear diffusion models are developed from the following model: (1)∂u∂t-∇·φx,y∇u=0,by changing the form of φx,y, which decides the degree of denoising and preservation of singularities (e.g., noise, edges, and texture). The isotropic diffusion expressed by a linear heat equation has been replaced by the anisotropic diffusion since the classical model was proposed by Perona and Malik in 1990.It is well known that the PM model yields very favorable results for eliminating image noise while retaining object edges and reducing oscillations. However, it also possesses some unsatisfactory properties in some cases; that is, it easily loses the contrast and texture details and produces staircase effects and so on [6, 7]. To remedy the loss of contrast and texture details, some scholars proposed an iterative regularization method for image restoration [8], and L2 norm in the fidelity term in (1) was replaced by L1 norm in [9–12]. Staircase effect is one of the unsatisfactory properties of the PM model. On the other hand, in [13–17], some models for image noise removal have been introduced to prevent staircase effect, which are based on higher order derivatives to energy norm or integrate high-order derivatives to original model. In addition, a total generalized variation (TGV) [18, 19] has been presented, which contains order greater than or equal to 2, unlike TV, which only considers first-order derivatives. It has some effects on image denoising, which not only leads to piecewise polynomial intensities but also obtains almost perfect edges while efficiently avoiding staircase effect.Recently, geometric attributes of images have been incorporated into the image denoising processing. The most commonly used local geometrical features are the curvature of the image level curve [20, 21], the Gauss curvature of the image surface [22–25], and the mean curvature of the image surface [26]. Although the aforementioned models have some effects on retaining image contrast, textures, and removing staircase effects, the results are still unsatisfactory. In short, the previously proposed models are integer-order differential based model in essence, they may produce blur edges, and their effect for textures preservation is not well. In particular, the variational method based high-order derivatives of the denoised image often lead to speckle artifacts, that is, isolated black and white speckles.Thus, in order to preserve textures and achieve a good tradeoff between eliminating speckle artifacts and restraining staircase effect, as a compromise between second-order PDEs and high-order PDEs, fractional-order PDEs based methods have been investigated and applied in the fields of computer vision. Unlike our familiar derivative operator, that is the integer-order derivative operator, the fractional-order derivative operator expands the integer-order derivative operator; it possesses the “nonlocal” property because the fractional-order derivative at a point depends upon the characteristics of the entire function and not just the values in the vicinity of the point [27–30], which is beneficial to improve the performance of textures preservation. Bai and Feng (BF) [31] present an anisotropic diffusion equation based on fractional-order space derivatives for image denoising. Their proposed model is interpolated between the PM equation and the fourth-order anisotropic diffusion equation, which is modeled by(2)∂u∂t=-Dxα∗cDxαu2Dxαux,y-Dyα∗cDyαu2Dyαux,y,where cs=1/1+s and Dxα∗ and Dyα∗ are the adjoints of Dxα and Dyα, respectively. The numerical results implemented fractional developmental equation in frequency domain and showed that both of the staircase effect and the speckle artifacts can be eliminated effectively. Guidotti and Lambers [32] extended gradient operator of energy norm from first-order to fractional-order. It has some effects on image denoising. Zhang et al. [33, 34] generalized the classical integer-order TV regularization model for image denoising, which is based on the Grünwald-Letnikov fractional-order derivative. Numerical examples show that the new models can eliminate both staircase effect and speckle artifacts when the parameters are chosen appropriately. Moreover, they are able to retain textures better than the TV model. For the case of calculation, Pu et al. [35] scheduled a class of fractional-order differential masks and illustrated that fractional-order differentiation can handle well fine structures such as texture. With regard to the tomographic reconstruction of the blurry and noisy binary images, a variational approach based on adapted fractional-order Hilbert spaces has been proposed by Bergounioux and Trélat [36]. A class of fully fractional anisotropic diffusion models was presented by Janev et al. [37] for noise removal, where spatial as well as time fractional-order derivatives are employed. Cuesta et al. [38] deduced a Volterra matrix-valued equation for image noise removal from a generalization of the linear fractional integral equation. However, they still have many great potential drawbacks. They simply used the first-order gradient operator of the energy norm to displace fractional-order one and the algorithms do not consider the fractional power of the energy norm, so the loss of image contrasts is not kept well after denoising.In this paper, we aim to use two new techniques, fractional calculus and difference curvature to image denoising. In comparison to the Gauss curvature and the mean curvature, the difference curvature can more effectively distinguish between ramps and edges. The fractional-order derivatives are able to better deal with the textures of the image and achieve a good tradeoff between eliminating speckle artifacts and restraining staircase effect on the entire image. We combine the difference curvature and fractional calculus frames for one diffusion. Near object edges, the difference curvature is large and the total variation term approximates to theL1 norm of fractional-order gradient flow for the image in order to preserve the edges; away from the edges (in ramp and flat regions), the difference curvature is small and the total variation term approximates to the L2 norm of fractional-order gradient flow for the image in order to prevent staircase effect. It is confirmed that the difference curvature is able to distinguish effectively between ramps and edges. Therefore the proposed model can effectively distinguish between ramps and edges, while improving the performance of texture preservation and preventing staircase effect and speckle artifacts.The remainder of this paper is organized as follows. In Section2, the fractional-order derivatives and the new edge indicator (difference curvature) are briefly reviewed. In Section 3, based on fractional-order derivatives and the difference curvature the proposed model and algorithm are presented. Section 4 develops a numerical algorithm to implement our proposed model. Section 5 is devoted to experimental examples. Comparing the simulation results obtained by our proposed denoising model with the baseline denoising methods reports a noticeable improvement in the quality of the denoised images evaluated qualitatively and quantitatively. Section 6 is a conclusion.


## 2. Preliminaries

### 2.1. Fractional-Order Derivatives

As a generalization form of the integer-order derivative, the fractional-order derivative has a long history, but it has no uniform definition until now. Three class definitions, Grümwald-Letnikov definition, Riemann-Liouville definition, and Caputo definition [29], are generally carried out by some investigations in recent literatures. Let us introduce the representation of an n-fold integral Jn for f∈L10,+∞ as follows:(3)Jnft=1n-1!∫0tt-τn-1fτdτ.Replacing n∈N in (3), we obtain a definition of a fractional-order integral(4)Jαft=1Γα∫0tt-τα-1fτdτ,0<t<+∞,where Γα is Euler’s continuous gamma function. In particular, if α∈N, then Γα=α-1!. Let n=α, the operator Dα defined as(5)Dαft=DnJn-αftfor 0<t<+∞ is called the Riemann-Liouville (R-L) differential operator with order α. However, in this paper, we apply the frequency domain definition to extend the step from integer to fractional, because it is easier for calculating. Therefore, the fractional-order partial derivatives can be formulated as (6)Dxαux,y=12π∫R2eixω1+yω2iω1αu^ω1,ω2dω1dω2,Dyαux,y=12π∫R2eixω1+yω2iω2αu^ω1,ω2dω1dω2.It is natural that the 2D fractional-order gradient operator is able to be formulated as Dαgx,y=Dxαg,Dyαg, with the associated magnitude Dαg=Dxαg2+Dyαg2.The definition of fractional-order Hilbert spaceHsRn goes by using the Fourier transform as follows [31, 36]. For arbitrary number s>0, define(7)HsRn={f∈L2Rn∣1+ξ2s/2f^ξ∈L2Rn}endowed with the norm(8)fHsRn=∫Rn1+ξ2sf^ξ2dξ1/2,where f^ξ is the Fourier transform of fξ. We denote by HsΩ∗ the dual of HsΩ.


### 2.2. Difference Curvature

It is generally appreciated that the gradient of an image cannot effectively distinguish between the image edges and the image ramps. It is well known that the second-order derivatives are able to effectively distinguish edges and ramps in the case of one-dimensional signal [14]. Similarly, in the two-dimensional case of images, in order to effectively distinguish edges and ramps, we will introduce a new feature descriptor calleddifference curvature (DC). We define(9)DC=uηη-uξξ,where uηη and uξξ represent the second-order derivatives in the direction of the gradient ∇u and in the direction perpendicular to ∇u, respectively, and the operator · denotes the absolute value. Here (10)uηη=ux2uxx+2uxuyuxy+uy2uyyux2+uy2,uξξ=uy2uxx-2uxuyuxy+ux2uyyux2+uy2.The value of DC is large at the object edges, while it is small in flat and ramp regions [39], so the edges can be distinguished from flat and ramp regions when the DC is computed for an image. The property analysis of the new feature descriptor is as follows.(a)
At image edges,uηη is large and uξξ is small, so DC is large.(b)
In flat and ramp regions,uηη and uξξ are both small, so DC is small.(c)
At isolated noise,uηη and uξξ are both large but they are almost equal, so DC is small.


## 2.1. Fractional-Order Derivatives

As a generalization form of the integer-order derivative, the fractional-order derivative has a long history, but it has no uniform definition until now. Three class definitions, Grümwald-Letnikov definition, Riemann-Liouville definition, and Caputo definition [29], are generally carried out by some investigations in recent literatures. Let us introduce the representation of an n-fold integral Jn for f∈L10,+∞ as follows:(3)Jnft=1n-1!∫0tt-τn-1fτdτ.Replacing n∈N in (3), we obtain a definition of a fractional-order integral(4)Jαft=1Γα∫0tt-τα-1fτdτ,0<t<+∞,where Γα is Euler’s continuous gamma function. In particular, if α∈N, then Γα=α-1!. Let n=α, the operator Dα defined as(5)Dαft=DnJn-αftfor 0<t<+∞ is called the Riemann-Liouville (R-L) differential operator with order α. However, in this paper, we apply the frequency domain definition to extend the step from integer to fractional, because it is easier for calculating. Therefore, the fractional-order partial derivatives can be formulated as (6)Dxαux,y=12π∫R2eixω1+yω2iω1αu^ω1,ω2dω1dω2,Dyαux,y=12π∫R2eixω1+yω2iω2αu^ω1,ω2dω1dω2.It is natural that the 2D fractional-order gradient operator is able to be formulated as Dαgx,y=Dxαg,Dyαg, with the associated magnitude Dαg=Dxαg2+Dyαg2.The definition of fractional-order Hilbert spaceHsRn goes by using the Fourier transform as follows [31, 36]. For arbitrary number s>0, define(7)HsRn={f∈L2Rn∣1+ξ2s/2f^ξ∈L2Rn}endowed with the norm(8)fHsRn=∫Rn1+ξ2sf^ξ2dξ1/2,where f^ξ is the Fourier transform of fξ. We denote by HsΩ∗ the dual of HsΩ.


## 2.2. Difference Curvature

It is generally appreciated that the gradient of an image cannot effectively distinguish between the image edges and the image ramps. It is well known that the second-order derivatives are able to effectively distinguish edges and ramps in the case of one-dimensional signal [14]. Similarly, in the two-dimensional case of images, in order to effectively distinguish edges and ramps, we will introduce a new feature descriptor calleddifference curvature (DC). We define(9)DC=uηη-uξξ,where uηη and uξξ represent the second-order derivatives in the direction of the gradient ∇u and in the direction perpendicular to ∇u, respectively, and the operator · denotes the absolute value. Here (10)uηη=ux2uxx+2uxuyuxy+uy2uyyux2+uy2,uξξ=uy2uxx-2uxuyuxy+ux2uyyux2+uy2.The value of DC is large at the object edges, while it is small in flat and ramp regions [39], so the edges can be distinguished from flat and ramp regions when the DC is computed for an image. The property analysis of the new feature descriptor is as follows.(a)
At image edges,uηη is large and uξξ is small, so DC is large.(b)
In flat and ramp regions,uηη and uξξ are both small, so DC is small.(c)
At isolated noise,uηη and uξξ are both large but they are almost equal, so DC is small.


## 3. Difference Curvature-Based Fractional Anisotropic Diffusion Model

### 3.1. Description of Our Proposed Model

Let us recall a classic example of noise removal scheme driven by curvature, which is an isotropic diffusion scheme, where the smoothing process is formulated by a partial differential equation. Letu denote an image intensity function, t is the time (scale parameter), and c· denotes the speed function which is dependent on the curvature κ of the level curve; then the diffusion equation can be formulated as (11)∂u∂t=cκ∇u.If cκ=κ, (11) is the standard curvature evolution equation [21].With introducing a new coordinatez and assigning the grey level of the image to the coordinate, we obtain a level function ψx,y,z=ux,y-z. Then the image ux,y can be viewed as its zero level set x,y,z:ψx,y,z=0 corresponds to the surface z=ux,y. In [23], the mean curvature of image surface evolution for noise removal is proposed, which is modeled by(12)∂u∂t=∇·11+∇u2∇u=∇·11+ux2+uy2∇ux,y.A modified mean curvature motion for image smooth is proposed in [24]. Chan et al. [25] developed a diffusion equation which uses the mean curvature as the conductance term, that is, the mean curvature driven diffusion (MCDD) in which the amount of diffusion is dependent on the mean curvature, which can be written as(13)∂u∂t=∇·M∇ux,y,where M is the mean curvature.Lee and Seo [26] proposed a Gauss curvature driven diffusion (G-CDD) model in which the amount of diffusion is determined by the Gauss curvature; that is,(14)∂u∂t=∇·ϕuxxuyy-uxy21+ux2+uy22∇ux,y,where ϕ:R→[0,+∞] is a decreasing and continuous function, ϕ∇u≥0 with ϕ0=0, and uxxuyy-uxy2/1+ux2+uy22 is the Gauss curvature of the image surface z=ux,y.It is well known that the mean curvature and the Gauss curvature are the average and product ofuηη and uξξ, respectively, which can be defined as(15)MC=uηη+uξξ2,GC=uηη·uξξ,where MC is mean curvature and GC is Gauss curvature.Let us discuss the performance comparisons of gradient and three curvatures, that is, mean curvature, Gauss curvature, and difference curvature. Compared with the flaw of anisotropic diffusion which has a tendency to estimate the false edges during the evolution, the difference curvature driven diffusion for images is more appropriate than the gradient driven diffusion to detect the intensity oscillations because the difference curvature can distinguish the edges and ramp regions; however the gradient difference between edges and ramps is not significant.Although the mean curvature and the difference curvature can effectively distinguish edges and ramps, the mean curvature cannot distinguish edges and isolated noise in a noisy image. For the Gaussian curvature, especially when the image is with isolated noise, the values at corners are large, so the denoising based on the Gaussian curvature driven diffusion cannot retain corners. Meanwhile, the gradient and three curvatures both have a flaw of forming piecewise constant images and indistinctness of texture details because they use the integer-order derivative.Inspired by fractional-order anisotropic diffusion and the (mean/Gauss) curvature driven diffusion by introducing fractional calculus instead of the (mean/Gauss) curvature, we propose a novel diffusion model for image denoising which uses the difference curvature in the conductance term to describe the intensity variations in images, which is named difference curvature driven fractional-order anisotropic diffusion (DCFAD) and modeled by(16)∂u∂t=-Dxα∗ϕDCDxαux,y-Dyα∗ϕDCDyαux,y,where Dxα∗ and Dyα∗ are the adjoints of Dxα and Dyα, respectively, ϕDC=exp⁡-DC/K is conductance (diffusivity) function based on Laplacian kernel [40], and K is a threshold constant. As shown in Figure 1, compared with two conductance functions of PM diffusion model, the diffusion speed of conductance function based on Laplacian kernel is between them at object edges (where DC is large), and it is smaller than them in flat and ramp regions (where DC is large).Figure 1
Comparison of three conductance functions.


## 3.1. Description of Our Proposed Model

Let us recall a classic example of noise removal scheme driven by curvature, which is an isotropic diffusion scheme, where the smoothing process is formulated by a partial differential equation. Letu denote an image intensity function, t is the time (scale parameter), and c· denotes the speed function which is dependent on the curvature κ of the level curve; then the diffusion equation can be formulated as (11)∂u∂t=cκ∇u.If cκ=κ, (11) is the standard curvature evolution equation [21].With introducing a new coordinatez and assigning the grey level of the image to the coordinate, we obtain a level function ψx,y,z=ux,y-z. Then the image ux,y can be viewed as its zero level set x,y,z:ψx,y,z=0 corresponds to the surface z=ux,y. In [23], the mean curvature of image surface evolution for noise removal is proposed, which is modeled by(12)∂u∂t=∇·11+∇u2∇u=∇·11+ux2+uy2∇ux,y.A modified mean curvature motion for image smooth is proposed in [24]. Chan et al. [25] developed a diffusion equation which uses the mean curvature as the conductance term, that is, the mean curvature driven diffusion (MCDD) in which the amount of diffusion is dependent on the mean curvature, which can be written as(13)∂u∂t=∇·M∇ux,y,where M is the mean curvature.Lee and Seo [26] proposed a Gauss curvature driven diffusion (G-CDD) model in which the amount of diffusion is determined by the Gauss curvature; that is,(14)∂u∂t=∇·ϕuxxuyy-uxy21+ux2+uy22∇ux,y,where ϕ:R→[0,+∞] is a decreasing and continuous function, ϕ∇u≥0 with ϕ0=0, and uxxuyy-uxy2/1+ux2+uy22 is the Gauss curvature of the image surface z=ux,y.It is well known that the mean curvature and the Gauss curvature are the average and product ofuηη and uξξ, respectively, which can be defined as(15)MC=uηη+uξξ2,GC=uηη·uξξ,where MC is mean curvature and GC is Gauss curvature.Let us discuss the performance comparisons of gradient and three curvatures, that is, mean curvature, Gauss curvature, and difference curvature. Compared with the flaw of anisotropic diffusion which has a tendency to estimate the false edges during the evolution, the difference curvature driven diffusion for images is more appropriate than the gradient driven diffusion to detect the intensity oscillations because the difference curvature can distinguish the edges and ramp regions; however the gradient difference between edges and ramps is not significant.Although the mean curvature and the difference curvature can effectively distinguish edges and ramps, the mean curvature cannot distinguish edges and isolated noise in a noisy image. For the Gaussian curvature, especially when the image is with isolated noise, the values at corners are large, so the denoising based on the Gaussian curvature driven diffusion cannot retain corners. Meanwhile, the gradient and three curvatures both have a flaw of forming piecewise constant images and indistinctness of texture details because they use the integer-order derivative.Inspired by fractional-order anisotropic diffusion and the (mean/Gauss) curvature driven diffusion by introducing fractional calculus instead of the (mean/Gauss) curvature, we propose a novel diffusion model for image denoising which uses the difference curvature in the conductance term to describe the intensity variations in images, which is named difference curvature driven fractional-order anisotropic diffusion (DCFAD) and modeled by(16)∂u∂t=-Dxα∗ϕDCDxαux,y-Dyα∗ϕDCDyαux,y,where Dxα∗ and Dyα∗ are the adjoints of Dxα and Dyα, respectively, ϕDC=exp⁡-DC/K is conductance (diffusivity) function based on Laplacian kernel [40], and K is a threshold constant. As shown in Figure 1, compared with two conductance functions of PM diffusion model, the diffusion speed of conductance function based on Laplacian kernel is between them at object edges (where DC is large), and it is smaller than them in flat and ramp regions (where DC is large).Figure 1
Comparison of three conductance functions.


## 4. Numerical Implementation

In this section, to solve (16), we introduce the frequency domain definition for the fractional-order derivative because it is easy to implement. We suppose the original discrete image has n×n pixels and is uniformly sampled from its continuous version starting at 0,0. It is a finite and discrete 2D signal, which stands for the actual image pixels in n×n resolution. We employ the 2D discrete Fourier transform (2D DFT) to approximate the spatial fractional-order difference. Note that we usually handle an image on a bounded domain, while the Fourier domain fractional-order derivative is defined on the whole space. Thus, the image should be prolonged onto the whole space in practical computations.Fractional central partial difference formulation (see [31]) can have the following form:(17)D~xαu=F-11-exp⁡-j2πω1nα×exp⁡-jπαω1nu^ω1,ω21-exp⁡-j2πω1nα,(18)D~yαu=F-11-exp⁡-j2πω2nα×exp⁡-jπαω2nu^ω1,ω21-exp⁡-j2πω2nα,where ω1,ω2∈0,1,…,n-1 are DFT frequencies that correspond to x and y, respectively. F-1 is the inverse 2D DFT operator.The adjoints of operatorsD~xαu and D~yαu are calculated by the following formulation, respectively:(19)D~xα∗u=F-1conj1-exp⁡-j2πω1nα×exp⁡-jπαω1nu^ω1,ω21-exp⁡-j2πω1nα,D~yα∗u=F-1conj1-exp⁡-j2πω2nα×exp⁡-jπαω2nu^ω1,ω21-exp⁡-j2πω2nα,where conj· is the complex conjugation.OperatorD~xαu has the form D~xαu=F-1∘K1∘F, where (20)K1=diag1-exp⁡-j2πω1nαexp⁡-jπαω1nstands for a diagonal matrix and ∘ represents the matrix multiplication. The adjoint operator D~xα∗u of D~xαu is calculated by the following formulation:(21)D~xα∗u=F-1∘K1∘F∗=F-1∘K1∗∘F,where(22)K1∗=diagconj1-exp⁡-j2πω1nαexp⁡-jπαω1n.The same is obtained forD~yαu and D~yα∗u. As is proven in [31], when the image size n is an odd integer, then D~xαu is real valued. When n is an even integer, the domain of function pω1=1-exp⁡-j2πω1/nα is not symmetric with respect to ω1=n/2, so usually there is a complex component in D~xαu. However, the image size is large enough (e.g., n=512, as used in the experiments); this complex component is very small due to the fact that the symmetry is minimally disturbed so that it can be ignored in computations.To solve (16), the finite difference schemes are applied as follows. Let Δt be the time step; we calculate the iterative scheme of the initial image ux,y and work in the Fourier domain, only returning to the spatial domain when calculating lxn=ϕDCDxαun and lyn=ϕDCDyαun, where ϕDC=exp⁡-DC/K.We suppose that the imageut is defined on 1h,2h,…,Nh×1h,2h,…,Nh, where h is the space step size. We use the notation ui,jn=unτ,ih,jh. We denote ui,jn=0 for i,j<1 or i,j>N.The discrete approximation is defined as follows:(23)Dxui,jn=ui+1,jn-ui,jnh,Dyui,jn=ui,j+1n-ui,jnh,Dxxui,jn=ui+1,jn+ui-1,jn-2ui,jnh,Dyyui,jn=ui,j+1n+ui,j-1n-2ui,jnh,Dxyui,jn=ui+1,j+1n+ui-1,j-1n-ui+1,j-1n+ui-1,j+1nh,uηηi,jn=Dxui,jn2Dxxui,jn+2Dxui,jnDyui,jnDxyui,jn+Dyui,jn2Dyyui,jn·Dxui,jn2+Dyui,jn2-1,uξξi,jn=Dyui,jn2Dxxui,jn-2Dxui,jnDyui,jnDxyui,jn+Dxui,jn2Dyyui,jn·Dxui,jn2+Dyui,jn2-1,DCi,jn=uηηi,jn-uξξi,jn.Our noise removal method for solving (16) is as follows.Step 1.
Let the input image beu and set n=1, un=u, k, Δt, and t=kΔt; calculate the 2D DFT u^n of un.Step 2.
Calculate theα-order central partial differences D~xαun and D~yαun using (18).Step 3.
Calculatelxn=ϕDCi,jnD~xαun and lyn=ϕDCi,jnD~yαun in the spatial domain.
Calculated^n=K1∗∘Flxn+K2∗∘Flyn.Step 4.
Calculateu^n+1=u^n-d^n×Δt and set n=n+1; if n=k, calculate the 2D IDFT of u^n and then stop; else go to Step 2.


## 5. Experimental Results

To analyze and show the good denoising capability of our proposed difference curvature driven fractional-order anisotropic diffusion (DCFAD) for image denoising, some experiments are carried out to compare the denoising result of the proposed method with that of several baseline denoising methods, including the PDE-based methods, the PM model [4], fourth-order PDEs [13], mean curvature-based variational model proposed by Zhu and Chan (ZC) [22], Gauss curvature driven diffusion (G-CDD) [26], the fractional anisotropic diffusion proposed by Bai and Feng (BF) [31], and the other class methods: nonlocal means filtering (NLMF) [41], the wavelet-based method (Bayes least-squares estimate of Gaussian scale mixture, BLS-GSM) [42], and the total generalized variation (TGV) [19]. The proper selection of the fractional-order α of our proposed model depends on the noisy image. We have tested the order from 0.1 to 3.0 with step size of 0.2 on every testing image and we find that α=1.8 can get the best results on the “Barbara” image, α=1.6 can get the best results on the “Bicycle” image, and α=1.8 can get the best results on the “Baboon” image. The time step is set to be Δt=4-α [35] and K=30. For our model, the peak signal to noise ratio (PSNR) and the mean structural similarity (MSSIM) [43] are adopted to illustrate denoising performance and structure preservation in terms of objective criterion and visual fidelity. We adopt a uniform stopping criterion for all the baseline denoising methods and the proposed method; that is, PSNRun+1<PSNRun or the iteration counter becomes larger than 5000 otherwise, where un denotes the denoised image in the nth iteration. The results of all methods are obtained by performing under MATLAB R2011a running on a desktop with Intel Core i3-2120 CUP 3.30 GHz and 8 GB of memory under 64-bit Windows 7.Figure2 shows the results of comparison of some denoising methods on a digital image “Barbara.” The values of PSNR and MSSIM on image “Barbara” are presented in Table 1. In Figure 2(a), the original image is presented and the noisy image with stand variance σ=20 is shown in Figure 2(b). Figures 2(c)–2(g) show the denoising results with the PM model, fourth-order PDEs, ZC, G-CDD, BF, NLMF, BLS-GSM, TGV, and our proposed DCFAD model, respectively. In order to perform better than the other models, we present the denoising results of a small patch of the original noisy image in Figure 3.Table 1
Comparison of the denoising performance for the testing images using different methods with three different values of Gaussian noise standard deviationσ. Results of a consistent PSNR (the first line) and MSSIM (the second line) are listed, respectively.
Image
σ
Denoising methods
PM
Fourth-order PDEs
ZC
G-CDD
BF
NLMF
BLS-GSM
TGV
Proposed
Barbara (512 × 512)
10
33.08
33.27
33.52
33.45
34.31
32.58
35.06
34.52
35.79
0.921
0.923
0.942
0.935
0.946
0.909
0.958
0.951
0.974
20
30.86
30.68
31.44
31.23
31.63
30.42
32.38
32.08
32.71
0.861
0.854
0.901
0.896
0.905
0.851
0.917
0.913
0.926
30
27.26
27.89
28.41
28.17
29.21
27.65
29.36
29.26
29.86
0.810
0.822
0.847
0.831
0.861
0.815
0.870
0.863
0.885
Bicycle (512 × 512)
10
29.13
29.28
30.32
30.20
30.65
30.26
31.10
30.84
31.62
0.845
0.851
0.869
0.858
0.870
0.866
0.876
0.871
0.889
20
24.74
24.68
26.82
25.88
27.27
26.05
27.31
27.30
27.52
0.723
0.721
0.781
0.731
0.801
0.777
0.808
0.806
0.826
30
22.67
22.81
24.89
23.11
25.17
23.14
25.23
25.21
26.04
0.648
0.651
0.720
0.655
0.735
0.656
0.749
0.742
0.761
Baboon (512 × 512)
10
29.12
29.27
29.89
29.68
30.11
28.99
30.18
30.13
30.26
0.873
0.876
0.888
0.882
0.882
0.831
0.890
0.884
0.895
20
24.86
25.00
25.46
25.29
25.86
24.48
26.03
25.96
26.12
0.743
0.744
0.768
0.759
0.766
0.607
0.773
0.769
0.780
30
22.80
22.85
23.32
23.16
23.75
21.13
23.88
23.80
24.14
0.645
0.647
0.663
0.658
0.661
0.456
0.670
0.664
0.692Figure 2
Comparison of the denoised results on “Barbara” image. (a) Original image. (b) Noisy image (additive Gaussian white noise, PSNR = 22.1 db). (c) PM model. (d) Fourth-order PDEs. (e) ZC model. (f) G-CDD. (g) BF. (h) NLMF. (i) BLS-GSM. (j) TGV. (k) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)Figure 3
Comparison of a local region of the cleaned images of “Barbara” image. (a) A region in the noisy image. (b) PM model. (c) Fourth-order PDEs. (d) ZC model. (e) G-CDD. (f) BF. (g) NLMF. (h) BLS-GSM. (i) TGV. (j) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)From Figure3 and Table 1, we can see the following performances. (1) A large scale cartoon part (face) shows staircase effect in the PM model because the flat regions have been smoothed very strongly. (2) The denoising result of the fourth-order PDEs shows speckle artifacts. (3) Small scale texture part (scarf) cannot be kept for the PM and fourth-order PDEs. (4) As seen from Figures 3(d) and 3(i), the noise of the ZC model at the high-frequency edges and textures regions (the stripes on Barbara’s scarf) is reserved; that is, it cannot distinguish edges and isolated noise in the noisy image, while it cannot keep small scale texture when the image includes a high level of noise. (5) In Figure 3(e), the result of the G-CDD model is able to distinguish object edges and ramps but retain the noise at the textures and corners, which will lead to the loss of small scale textures and corners when the image includes a high level of noise. (6) The results of Figures 3(f) and 3(i) exhibit that the BF model yields the false edges because the gradient indicator is sensitive to a high level of noise at the oscillations. (7) As shown in Figure 3(g), the denoising result obtained by the NLMF model is very good but smoothes the image so strongly that it leads to the blur of small scale textures. (8) From Figure 3(h), the wavelet-based method BLS-GSM yields artificial ringing artifacts when the model is reconstructed by wavelet. (9) Although we can see that the result of TGV model avoids blocky effects (staircase effects) near discontinuities in smooth regions and leads to a nature piecewise smooth, the TGV model cannot distinguish the edges and ramp regions because the gradient difference between edges and ramps is not significant. Meanwhile, the TGV contains order greater than or equal to two, which is not able to avoid speckle effects. (10) The values of MSSIM on image “Barbara” imply that our proposed model possesses good local structure preservation capability for the original image. (11) The denoising capability and the local structure-preserving capability of our proposed method can be improved when the fractional-order derivatives and the difference curvature are computed on the predenoised image. (12) Comparing with the denoising results of the other eight models, our proposed model has a good performance in terms of objective criterion and visual fidelity, which obtains better PSNR and is able to retain important features such as edges, ramps, corners, and textures.Figure4 shows the results of comparison of some denoising methods on “bicycle” image with many geometric shapes including disks and squares that possess significant amount of highly oscillatory regions (stripes of various thicknesses in the right upper corner, i.e., textures). The noisy image with stand variance σ=20 is shown in Figure 4(b). Table 1 presents the PSNR and MSSIM of Figure 4. In order to further study the performance of all models, we present the denoising results for a small patch of the original noisy image in Figure 5. It can be confirmed from the enlargements that our proposed model does not introduce other artifacts and has more distinct edges enhancement effect on texture regions than other methods. Thus the difference curvature calculated on the noisy image can improve the denoising performance. For the image denoising results, as it can be seen, PM method preserves edges and corners but has oversmoothed the stripes of various thicknesses in the right upper corner. The problem for fourth-order PDEs model is that speckle effects are not able to be avoided. The ZC model exhibits better visual results on the stripes of various thicknesses in the right upper corner in comparison to the first two methods, but some blurry regions appear in stripes of various thicknesses in the right upper corner in comparison to the proposed method. In Figure 5, the denoising result shows the G-CDD model can preserve the ramp but the denoising capability in the corners is much worse. Figure 5(f) shows the BF model loses a large number of high-frequency edges. As shown in Figure 5(g), the NLMF model smoothes strongly the stripes of various thicknesses in the right upper corner. The wavelet-based method BLS-GSM tends to show some artifacts in Figure 5(h). The TGV model exhibits some blurry regions on the stripes of various thicknesses in the left upper corner compared with the proposed method.Figure 4
Comparison of the denoised results on “Bicycle” image. (a) Original image. (b) Noisy image (additive Gaussian white noise, PSNR = 22.1 db). (c) PM model. (d) Fourth-order PDEs. (e) ZC model. (f) G-CDD. (g) BF. (h) NLMF. (i) BLS-GSM. (j) TGV. (k) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)Figure 5
Comparison of a local region of the cleaned images of “Bicycle” image. (a) A region in the noisy image. (b) PM model. (c) Fourth-order PDEs. (d) ZC model. (e) G-CDD. (f) BF. (g) NLMF. (h) BLS-GSM. (i) TGV. (j) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure6 shows the results on a “Baboon” image of size 512×512. The noisy image with stand variance σ=20 is shown in Figure 6(b). Figure 6(c) indicates that the PM model brings in a false edge at the nasolabial groove. The fourth-order PDEs method is able to effectively overcome the staircase effect, but speckle effects are introduced. The ZC model shows the corners preserving capability is worse. It can be observed that the G-CDD model cannot preserve edges at the nasolabial groove. Some blurry edges at the nasolabial groove appear in result of the BF model. The result obtained by the NLMF model smoothes strongly the flat regions of Baboon, which leads to the blur of the hairs on Baboon. The result from wavelet-based method BLS-GSM exhibits good denoising capability but creates some artifacts. As seen, our proposed model is able to better save the Baboon’s nasal surface in smooth areas and nonlinearly keep the bridge of the nose in the regions where the gray level changes greatly and also enhance Baboon’s hairs in those regions where the gray level does not change obviously. The TGV model exhibits speckle effects on Baboon’s nasal surface in smooth areas.Figure 6
Comparison of the denoised results on “Baboon” image. (a) Original image. (b) Noisy image (additive Gaussian white noise, PSNR = 22.1 db). (c) PM model. (d) Fourth-order PDEs. (e) ZC model. (f) G-CDD. (g) BF. (h) NLMF. (i) BLS-GSM. (j) TGV. (k) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)In order to perform better than the other models to preserve edges, ramps, corners, and textures, denoising results of a small patch of the original noisy image are demonstrated in Figure7. To show the practical applicability of our model, we test our proposed method on a 512×512 CT Scan medical image corrupted by the random white Gaussian additive noise with stand variance σ=20. In Figure 8, we show our model can retain the important features, such as curvy edges, straight edges, ramps, corners, and small scale textures.Figure 7
Comparison of a local region of the cleaned images of “Baboon” image. (a) A region in the noisy image. (b) PM model. (c) Fourth-order PDEs. (d) ZC model. (e) G-CDD. (f) BF. (g) NLMF. (h) BLS-GSM. (i) TGV. (j) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 8
Comparison of the denoised results on CT Scan image. (a) Original image. (b) Noisy image (additive Gaussian white noise, PSNR = 22.1 db). (c) PM model. (d) Fourth-order PDEs. (e) ZC model. (f) G-CDD. (g) BF. (h) NLMF. (i) BLS-GSM. (j) TGV. (k) Our proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)To illustrate theΔt does not depend on image size, we run our proposed algorithm on three benchmarked images, “Barbara,” “Lena,” and “Peppers,” with two different resolutions, 256×256 and 512×512. Table 2 reports the corresponding results on the MSSIM with σ=20. It can be confirmed from Table 2 that the MSSIM values obtained by the testing images with resolution of 256×256 are close to that of 512×512, respectively. It indicates that Δt in our adoptive algorithm does not depend on image size.Table 2
Comparison of MSSIM on three benchmarked images for different resolutions withσ=20.
Image size
Barbara
Lena
Peppers
256 × 256
0.9195
0.8567
0.8452
512 × 512
0.9261
0.8662
0.8538


## 6. Conclusions

In this paper, we present an introduction of two new techniques, fractional calculus and difference curvature, to the field of image processing and the implementation of fractional-order PDE, which can effectively distinguish between ramps and edges while preserving important structures such as edges, ramps, corners, and textures and preventing staircase effect and speckle artifacts. The experimental results show that our proposed method is capable of preserving the low frequency contour feature in the smooth area, nonlinearly maintaining the high-frequency edge and texture details in those areas where gray level changes greatly, and also without smoothing the important fine details in those areas where gray level has little changed while preserving corners and effectively distinguishing between ramps and edges.



---

*Source: 930984-2015-04-07.xml*