Fully fractional anisotropic diffusion for image denoising

Mathematical and Computer Modelling 54 (2011) 729–741

Contents lists available at ScienceDirect

Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm

Fully fractional anisotropic diffusion for image denoising Marko Janev a,∗ , Stevan Pilipović b , Teodor Atanacković a , Radovan Obradović a , Nebojša Ralević a a

Faculty of Engineering, University of Novi Sad, Trg Dositeja Obradovića 6, 21000 Novi Sad, Serbia

b

Department of Mathematics and Informatics, University of Novi Sad, Trg Dositeja Obradovića 4, 21000 Novi Sad, Serbia

article

info

Article history: Received 6 January 2010 Received in revised form 25 January 2011 Accepted 15 March 2011 Keywords: Anisotropic diffusion Fractional derivatives Fractional ordinary differential equations Fractional linear multistep methods Fractional order differences Image denoising

abstract This paper introduces a novel Fully Fractional Anisotropic Diffusion Equation for noise removal which contains spatial as well as time fractional derivatives. It is a generalization of a method proposed by Cuesta which interpolates between the heat and the wave equation by the use of time fractional derivatives, and the method proposed by Bai and Feng, which interpolates between the second and the fourth order anisotropic diffusion equation by the use of spatial fractional derivatives. This equation has the benefits of both of these methods. For the construction of a numerical scheme, the proposed partial differential equation (PDE) has been treated as a spatially discretized Fractional Ordinary Differential Equation (FODE) model, and then the Fractional Linear Multistep Method (FLMM) combined with the discrete Fourier transform (DFT) is used. We prove that the analytical solution to the proposed FODE has certain regularity properties which are sufficient to apply a convergent and stable fractional numerical procedure. Experimental results confirm that our model manages to preserve edges, especially highly oscillatory regions, more efficiently than the baseline parabolic diffusion models. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Since the work of Perona and Malik [1], PDE methods have been used for image processing, especially for denoising and stabilizing edges (see [1,2]). Perona and Malik were the first to replace isotropic diffusion expressed through a linear heat equation with an anisotropic diffusion. Anisotropic diffusion is associated with an energy dissipating process that seeks the minima of an energy functional. For example, the well-known total variation (TV) minimization model [3,4] can be obtained in the case when the energy functional is equal to the TV norm of the image. Although these methods have been demonstrated to be able to achieve a good trade-off between noise removal and edge preservation, the resulting image in the presence of the noise often has a ‘‘blocky’’ look. It is caused by the use of a second order PDE modeling methods. In order to reduce the ‘‘blocky effect’’, while preserving sharp jump discontinuities, many other non-linear filters have been suggested in the literature (see [5–9]). In [5], You and Kaveh proposed a class of fourth order PDEs which are obtained by the minimization of a functional given as an increasing function of the edge detector △u. Since the second order derivatives are zero if the image intensity function is planar, the class of fourth order PDEs will evolve and settle down to a planar image if the image support is infinite. This is important since piecewise planar images look more natural than the step images which are stationary points of the particular non-convex energy functional [5], whose minimization (after the gradient descent is applied) leads to the second order diffusion. The problem with the use of fourth order equations is that

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (M. Janev).

0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.03.017

730

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

it tends to leave the image with isolated black and white speckles (so called ‘‘speckle effect’’) which may be characterized as pixels whose intensity values are either much larger or much smaller than those of the neighboring pixels, as explained in [5]. Recently, fractional order PDEs have been studied and applied to the problem of image denoising. Bai and Feng in [10] proposed the use of a nonlinear anisotropic fractional diffusion equation based on the Euler–Lagrange equation of a cost functional which is an increasing function of the absolute value of fractional derivative of the image intensity function. They managed to interpolate between second and fourth order nonlinear anisotropic diffusion equation in order to obtain a more natural images. Cuesta [11] proposed the use of a linear integro-differential equation in order to control the diffusion effect in image enhancement. He was able to interpolate between the heat equation, where diffusion effect is the most present, and the wave equation so that the edges could be better preserved. In this paper, two mentioned approaches are used in order to gain a novel fully fractional nonlinear PDE which can interpolate between the parabolic and the hyperbolic PDE and, at the same time, between the second and the fourth order PDE. Thus, better image enhancement is obtained with an interpolation using fractional spatial derivatives, as well as using the edge stopping function. Moreover, in order to preserve edges in a highly oscillatory regions (not preserved by the heatlike diffusion), fractional time derivatives are used. Moreover we use an appropriate numerical scheme for which we have proved that it is convergent and stable. The numerical scheme is given in the framework of Fractional Linear Multistep Method (FLMM) applied to the fractional ODE model obtained by spatial discretization of the PDE. The paper is organized as follows. Basic definitions of fractional order integrals and derivatives are presented in Section 2. Anisotropic diffusion by means of spatial fractional derivatives is presented in Section 3 while in Section 4, an additional fractional order time derivative is introduced. In Section 5 we turn our attention to the discretization of the proposed continuous model. A numerical method that resolves the discretized equation is presented, and the consistency and the stability of the numerical scheme are proven. Experimental results are presented in Section 6. They show that in the proposed approach, the choice of the fractional parameters can lead to a better image enhancement in comparison to the baseline anisotropic diffusions. 2. Fractional order derivatives and integrals It is well known that for f ∈ L1 ([a, b]), an n-fold integral Jan has the representation: Jan f (t ) =

t

∫

1

(n − 1)!

(t − τ )(n−1) f (τ )dτ

(1)

a

Replacing n ∈ N in (1) with β ∈ R+ , and by using the Euler’s gamma function instead of the factorial, we obtain the left β fractional integral operator of order α Ja defined on L1 ([a, b]) by Jaβ f (t ) =

1

Γ (β)

t

∫

(t − τ )(α−1) f (τ )dτ ,

t ∈ [a, b].

(2)

a

β

For β = 0, we set Ja0 = I. It can be shown [12,13] that if f ∈ L1 ([a, b]) and β > 0, then the integral Ja f (t ) exists for almost every t ∈ [a, b]. β If we set n = ⌈β⌉, the operator Da defined by Dβa f (t ) = Dn Jan−β f (t )

(3)

for x ∈ [a, b] is called the left Riemann–Liouville (RL) differential operator of order β . For β = 0 we β

set D0a

= I. If f is n-times

absolutely differentiable on [a, b] i.e. f ∈ AC n ([a, b]), then [12,13] Da f exists almost everywhere on [a, b]. More appropriate for use in the differential equations (initial conditions have clear physical interpretation, which we β explain in Section 4) is the left Caputo fractional differential operator C Da defined by C

Dβa f (t ) = Jan−β Dn f (t ),

t ∈ [a, b]

(4)

which we consider as the time derivative in Section 4 in our PDE model, where we fix the interval β ∈ [1, 2). 3. Anisotropic diffusion by means of spatial fractional derivatives The main drawback of the isotropic diffusion is that the smoothing can damage image features such as edges and lines. To avoid such damages, the amount and the direction of smoothing have to be controlled. To overcome this drawback, Perona and Malik [1] derived a nonlinear extension of the heat equation, leading to an anisotropic diffusion. They considered equation ut = ∇(c (|∇ u|2 )∇ u), where c : R → [0, 1] is a decreasing and continuous function (Edge stopping function or diffusion coefficient) vanishing on the edges (high gradients), and close to 1 on regular regions (low gradients).

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

731

This equation can be interpreted as an energy-dissipating process that seeks the minima of the energy functional E (u) =

∫ Ω

f (|∇ u|)dxdy.

(5)

In the sequel we assume that the boundary ∂ Ω of the region Ω is smooth enough. Moreover, we assume that f is regular enough, f (|∇ u|) ≥ 0 and that f is an increasing function of |∇ u| satisfying f ′ (|∇ u|) = 0 for |∇ u| = 0. Functions c and f are related by c (s) =

f′

√  √

s

s

,

s > 0.

(6)

Bai and Feng [10] generalized (5) introducing the fractional spatial derivative Dαx,y u = (Dαx , Dαy ), for α ∈ [1, 2], instead of the gradient vector. Here Dαx and Dαy , for α > 0, are defined as pseudo-differential operators with necessary assumptions on u, which we will always assume to be satisfied. Dαx u(x, y) =

1

∫

ei(xξ1 +yξ2 ) (iξ1 )α uˆ (ξ1 , ξ2 )dξ1 dξ2 2π R2 ∫ 1 Dαy u(x, y) = ei(xξ1 +yξ2 ) (iξ2 )α uˆ (ξ1 , ξ2 )dξ1 dξ2 , 2π R2

(7)

(x, y) ∈ R2 ,

where we assumed that u is extendable out of Ω and that its Fourier transform uˆ (ξ1 , ξ2 ) exists in R2 . The minimization of the fractional energy functional

∫ Ω

f (|Dαx,y u|)dxdy,

α ∈ [1, 2]

(8)

⃗ = 0 on ∂ Ω , (n⃗ is the unit outward normal by the use of the partial integration with the imposed boundary condition Dαx,y u · n to ∂ Ω ), leading to the Euler–Lagrange equation (Dαx )∗ (c (|Dαx,y u|2 )Dαx u) + (Dαy )∗ (c (|Dαx,y u|2 )Dαy u) = 0. α ∗

(9) α

Note that (Dx ) in (9) is the adjoint operator of the linear operator Dx . As in [10], Eq. (9) can be solved via the following one (gradient descent procedure) ut = −(Dαx )∗ (c (|Dαx,y u|2 )Dαx u) − (Dαy )∗ (c (|Dαx,y u|2 )Dαy u),

(10)

where the original image u(0) = u0 is the initial condition. We note that the Eq. (10) is actually ill-posed [14]. This means that the existence and the uniqueness of the solution to (10) are not guaranteed. However, by introducing mollifier in the argument of diffusivity, i.e. by replacing |Dαx,y u| with |Dαx,y uσ |, where uσ = Gσ ∗ u and Gσ is n-dimensional Gaussian kernel with standard deviation σ > 0, one obtains the well-posed problem, i.e. existence and the uniqueness of the solution (see [14] or [15] for details). 4. Fully fractional anisotropic diffusion In this paper we propose a novel nonlinear PDE, the Fully Fractional Anisotropic Diffusion (FFAD), which combines temporal and spatial fractional derivatives in order to interpolate between the nonlinear heat equation and the hyperbolic nonlinear equation so that the amount of diffusion can be controlled. Thus, the central equation of our work is a generalization of (10), of the form: C

β

Dt u = −(Dαx )∗ (c (|Dαx,y u|2 )Dαx u) − (Dαy )∗ (c (|Dαx,y u|2 )Dαy u),

(11)

where β ∈ [1, 2), which is obtained from (10) by letting the time derivative to be fractional of order β . Actually, the left Caputo time fractional derivative of order β ∈ [1, 2) is used instead of ut . We use Caputo instead of the RL time fractional derivative. Let us explain, why?. We recall (see [13]) that in the case of Caputo time fractional derivative of order β ∈ [1, 2), one has to impose the initial conditions up to the second order. Thus, we put u(0) = u0 ,

u′ (0) = 0,

(12)

where these data have a physical interpretations: u(0) = u0 is actually an observed noisy image, and u (0) = 0 since the deformation of the image is zero, at t = 0. In the case of the RL time fractional derivative of order [1, 2), the initial conditions, that we have to impose on (11) are Dβ−1 u(0) = b1 and Dβ−2 u(0) = b2 , for some b1 , b2 ∈ R, without a clear meaning what b1 and b2 should be. Since we will use the approximations of the spatial fractional derivatives (imposed in Section 5), we replace the ‘‘no-flux’’ ¯ ⊂ R2 taken instead type boundary condition with the periodicity assumption for u on some closed rectangular domain Π on a general bounded domain Ω . ′

732

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

5. Approximations and the numerical procedure To our knowledge, contrary to the Eq. (10), where the well-posedness could be obtained by using mollifiers (as it is explained in Section 3), there are no results on the well-posedness of (11). Thus, similarly to Discrete Nonlinear Filtering in [15], instead to use and analyze directly Eq. (11), we use the approximate Spatially Semi-Discretized Equation, that is the equation obtained from (11) by applying the spatial discretization. Thus, we obtain the FODE system, for which we prove that it is well-posed, i.e. that there exists a unique solution. Since in (11) the time derivative is fractional, the model that we obtain requires completely different numerical approach than the one presented in [15]. In addition, we prove that our analytical solution to the FODE possesses certain regularity properties, so that we can state the convergent and stable fractional numerical scheme. 5.1. Spatial discretization We proceed by discretizing Eq. (11) in the spatial domain, similarly as it was done in [10]. We sample the original continuous image ucont of size D × D, where D ∈ R+ , by the M × M uniform grid, and obtain u(p, l) = ucont (p∆x, l∆y) for p, l ∈ {0, 1, . . . , M − 1}, where ∆x = ∆y = D/M. It is a finite, discrete, 2D signal, which represents the actual image pixels in M × M resolution. Also, the 2D DFT is used in order to obtain approximations of the spatial fractional derivatives. ¯ ⊂ R2 . Because of that we use the periodicity assumption for the image on the closed rectangular domain Π α Fractional central differences (see [10]) are used for the approximation of Dx u as follows: α u = F −1  D x dft

[

1−e

−2π iω1 M

α

e

iπαω1 M

]

Fdft [u(p, l)] ,

(13)

where ω1 ∈ {0, 1, . . . , M − 1} are DFT frequencies that correspond to p; the spatial coordinate is x and Fdft is DFT operator. αx has the form D αx = F −1 ◦ K1 ◦ Fdft , where Operator D dft K1 = diag



1−e

−2π iω1 M

α

e

iπαω1 M

 (14) ∗

αx of D αx is expressed [10] by is a diagonal matrix and ◦ stands for the matrix multiplication. The adjoint operator D ∗

αx = (F −1 ◦ K1 ◦ Fdft )∗ = F −1 ◦ K1∗ ◦ Fdft . D dft dft

(15)

∗

αy and D αy (see [10]), by replacing DFT frequencies ω1 that appears in (13) and (14), with ω2 The same is obtained for D that corresponds to l; the spatial coordinate is y. Relations (13)–(15) are used for the approximation of the Euler–Lagrange α u is real  (EL) part (right-hand side) of Eq. (11). As is proven in [10], when the image size M is an odd integer, then D x   −2π iω1 /M α valued. When M is an even integer, the domain of function p(ω1 ) = 1 − e is not symmetric with respect to α u (see also [10]). However, when image size is enough large  ω1 = ⌊M /2⌋, so usually, there is a complex component in D x (for example M = 512, as we use 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. As to images that are expected to appear in actual applications, they are usually much larger (e.g. 10 megapixels) and thus in their case the complex component is significantly smaller. Nevertheless, if M is even, and we want to eliminate this small complex component, we can extend the observed image u0 as follows [10]: uext 0 (p, l) = u0 (p, l),

0 ≤ p, l < M

uext 0 (p, M ) = u0 (p, M − 1), uext 0 (M , l) = u0 (M − 1, l), uext 0

(16)

0≤p
(M , M ) = u0 (M − 1, M − 1)

and use which has size (M + 1) × (M + 1), where M + 1 is odd, instead of u0 . We refer to [10] for the detail analysis concerning the spatial discretization. In our considerations we will also assume that c (s2 ) = f ′ (s)/s is analytic. This is true for the majority of the existing edge uext 0 ,

√

stopping functions in use (for example, Perona–Malik: c (s2 ) = (2/K )e−s /K , Minimal Surfaces: c (s2 ) = 2/ 1 + s2 GemanMcClure: c (s2 ) = 1/(1 + s2 ), etc.). Note that we use the Minimal Surfaces edge stopping function in the all experiments of Section 6. 2

5.2. Formulation of the spatially discretize fode system We proceed with the formulation of our spatially discretized FODE system. Let u(t ) be a vector obtained by concatenating all the columns of the matrix [u(t , p, l)] for some fixed t ∈ [0, T ], where T is the end of the time interval of interest and

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

733

αx , D αy defined by (13) in the form of matrices that act on the vector u(t ), p, l ∈ {0, . . . , M − 1}. In addition, we rearrange D αx )∗ , (D αy )∗ . Denote instead on the matrix [u(t , p, l)]. We do the same with (D α u(t )|2 ) = diag c (D α u(t ))2   c (|D x ,y x ,y j

 



j=1,...,N

,

(17)

α u(t ))2 = ((D  αx u(t ))j )2 + ((D αy u(t ))j )2 and N = M 2 . where (D x ,y j

αx , (D αx )∗ , D αy , (D αy )∗ to Eq. (11), the Euler–Lagrange part of the After applying approximate (finite) linear operators D equation is approximated as: 

G1 (u(t ))



 .. αx )∗ [hx (t )] + (D αy )∗ [hy (t )].  = (D . GM 2 (u(t ))

G(u(t )) = 

(18)



Terms hx (t ) and hy (t ) are given by α u(t )|2 )D  αx u(t ), hx (t ) = c (|D x ,y

α u(t )|2 )D  αy u(t ). hy (t ) = c (|D x ,y

(19)

We now state the approximate spatially discretized FODE system for (11): C

β

Dt u(t ) = G(u(t )),

u(t ) ∈ RN ,

u(0) = u0 ∈ RN ,

u′ (0) = 0,

(20)

where G is defined by (18) and (19). The following result is used for the formulation of our numerical procedure: Proposition 1. Function u : [0, T ] ⊂ R → RN is a solution to the spatially discretized FODE system (20) iff it is a solution to the following Volterra integral equation of the second kind: u( t ) = u0 +

1

Γ (β)

t

∫

(t − s)β−1 G(u(s))ds,

u( t ) ∈ R N ,

β ∈ [1, 2).

(21)

0

Moreover, there exists a unique continuous solution to (21). Proof. Function G : RN → RN given by (18) and (19) is Lipshitz continuous on every compact set K ⊂ RN since it is αx , (D αx )∗ , D αy , (D αy )∗ and by using the composition with the analytic edge obtained by applying bounded linear operators D stopping function c (s2 ). This implies that there exists a rectangular neighborhood U of u0 so that G is Lipshitz continuous ¯ As all the conditions for the application of Theorem 4.2.3 in [12], Chapter 4.2. are fulfilled, it follows that u ∈ C ([0, T ]) on U. is a solution of the nonlinear fractional ODE system iff it is a solution to the appropriate Volterra equation. The second part of the theorem follows directly from Lemma 4.2.5 in [12], where it is stated that under the same conditions as in Theorem 4.2.3 there exists a unique solution u ∈ C ([0, T ]) to (21).  Now we apply the result of Lubich [16,12] to the Volterra equation (22), in order to gain information on the behavior of the solution to (22) near the origin: Proposition 2. There exists a function Υ (z1 , z2 ) analytic in some neighborhood of (0, 0) ∈ C2 , so that u(t ) = Υ (t , t β ) for t ∈ [0, T ]. Proof. Eq. (21) is a Volterra equation of the form β−1

u( t ) = f ( t ) + J 0

[K (·, u(·))] (t ),

u(t ) ∈ RN ,

β ∈ [1, 2),

(22)

where f (t ) = u0 , and the kernel K (t , u) = G(u) is given by (18) and (19). The same arguments related to G given in αx )∗ , D αy , (D αy )∗ are bounded linear operators and c (s2 ) αx , (D Proposition 1 imply that G is an analytic function on RN , since D is analytic on RN . It follows that there exists a neighborhood U of u0 so that G is analytic in U. Moreover, adapting notation to Chapter 7 in [12] with f (t ) = F (t , t β ) = u0 , if we put F (x1 , x2 ) = u0 , it is obviously analytic in U. It follows that all the conditions for the application of the Theorem 4.2.7, Chapter 7 [12] are fulfilled. Thus, there exists a function Υ (z1 , z2 ) analytic in some neighborhood of (0, 0) ∈ C2 so that the solution to (21) can be represented as u(t ) = Υ (t , t β ).  Proposition 2 implies certain regularity of the solution, i.e. that the solution u to (21) has an asymptotic expansion in mix powers of t and t β of an arbitrary order, as t → 0, which we use in the next section. 5.3. Formulation of the numerical procedure We propose a convergent and stable implicit numerical scheme in order to obtain a numerical solution to (21), which is by Proposition 1, equivalent to the solution of the FODE system (20).

734

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

We follow [12,16]. Let ωj be given as the coefficients in the power expansion of ωβ (ξ ) = (ω(ξ ))β , where p − 1

(1 − ξ )k , ξ ∈ C, (23) k k =1 is the generating function of the Backward Difference (BD) [12,16] Linear Multi-step Method (LMM), and p ∈ {1, . . . , 6} (which we fix) is the order of convergence of BD. Let ωn,j be determined as the solution of the linear system of equations ω(ξ ) =

s −

ωn,j jγ =

j =0

n − Γ (1 + γ ) nγ −β − ωn−j jγ , Γ (1 + γ − β) j =1

γ ∈ A,

(24)

where A = {γ = k + jβ | k, j ∈ N0 , γ ≤ p − 1} and s = card(A) − 1. Recall that ωj are called convolution weights while ωn,j are called starting weights [12]. We obtain a numerical solution to (21) as a sequence n−1 −

un = hβ G(un ) −

ωn−j uj −

s −

j =0

−β

ωn,j uj + hβ

j =0

tn u0

Γ (1 − β)

,

(25)

where n = 1, . . . , N, tn = nh, N = T /h, and [0, T ] is interval of interest (see (21)). The proposed scheme (25) is an FLMM based on the result of Lubic presented in [16] (Fractional Backward Difference method), and additionally considered in [12] in Theorem 5.1.10, Chapter 5. The scheme (25) is obtained from the scheme presented in the Theorem 5.1.10, [12], by putting f (t , u) = G(u), where G is given by (18). Moreover, the last term in (25) is ∑⌈α⌉−1 obtained as the term k=0 bk tnk−α /Γ (k + 1 −α) in the quoted scheme, where we set ⌈α⌉ = 2, b0 = u0 and b1 = u′ (0) = 0. The mentioned theorem states, that if the analytical solution to (21) satisfies certain regularity assumptions, which in our case follow from Proposition 1), then the scheme given by (25) is convergent of order p. It means that the global error of the method on [0, T ] satisfies max{‖u(tn ) − un ‖2 |tn = nh, n = 1, . . . , N } = O(hp ),

as h → 0,

(26)

where u(tn ) is the exact solution to (21) at tn , while un is the numerical solution obtained by (25). It also states that the starting weights satisfy ωn,j = O(n−β−1 ), which we use in proving the stability of (25). We will give more details in order to explain how the regularity of the analytical solution to (21) effects the proposed scheme (25). Also, we will explain the role of the starting weights ωn,j . All those issues are explained in much more detail in [16], and in particular in [12], in Chapters 4 and 5. Theorem 4.3.12, [12] states the following: Let a∑ convergent LMM of order p ≥ 1 be given by its generating function ω(ξ ˜ ) = σ (ξ −1 )/ρ(ξ −1 ). Let also ω˜ β (ξ ) = (ω(ξ ˜ ))β = ∞ ˜ j ξ j for β > 0, j =0 ω and ωn,j satisfy (24), where A = ∪Lj=0 Aj , Aj = {γ = k + zj | k ∈ N0 , γ ≤ p − 1}. If the function y is regular enough, i.e., if it satisfies y(t ) =

L −

t zj vj (t ),

0 ≤ zj ≤ p − 1,

vj = C p [0, T ],

j = 0 , . . . , L,

(27)

j =0

then ωn,j = O(nβ−1 ), and the error of approximating the integral J β y(tn ) =

tn

∫

1

Γ (β)

(t − s)β−1 y(s)ds

(28)

0

with the convolution quadrature β Jh y(tn ) = hβ

n −

ω˜ n−j yj + hβ

j =0

s −

ω n ,j y j

(29)

j=0

behaves as O(hp ) when h → 0 (for β ≥ 1, which is our case). Note that Theorem 5.1.10, [12] is obtained as a reformulation of the Theorem 4.3.12, [12], by taking ω ˜ = ω given by (23) and replacing β with −β , while keeping the same necessary condition (27) on the regularity of the solution. Thus, condition (27) is also assumed in order to apply Theorem 5.1.10, [12] and to obtain the convergence of (25) as its direct consequence. We provide condition (27) as the consequence of Proposition 2. By Proposition 2, we obtain for arbitrary p ∈ N: u(t ) = Υ (t , t β ) =

L1 − L2 −

fk,j t k+jβ + O(hp )

k=0 j=0

=

L1 − k=0

=

L − k=0

fk t k

L2 −

gj t jβ + O(hp )

j =0

vj (t )t zj ,

zj = jβ

(30)

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

735

where vj , j = 1, . . . , L are polynomials of order p, v0 ∈ C p [0, T ], and zj ∈ A, so that the necessary regularity condition (27) in Theorems 4.3.12 and 5.1.10, [12] is satisfied. We recall (see proof of Theorem 4.3.12, [12]) that the degree of the regularity of the solution p that appears in (30) is actually the maximal order p ∈ {1, . . . , 6} of the underlying BD that could be used so that, in general, the regularity of the solution strongly affects the order of the convergence of the scheme (25). In our case, as p from Proposition 2 is arbitrary, we could obtain convergence of any order 1 ≤ p ≤ 6. We also recall, that the actual restriction 1 ≤ p ≤ 6 of the BD, is based on the fact that the BD is stable only up to order p = 6 (see [12] or [17]), and the stability of the underlying LMM method considered in Theorem 4.3.12, [12] is the necessary condition for the theorem to β hold. Moreover, the system (24) provides (J β − Jh )t γ = 0, γ = k + jβ ∈ A, so that, assuming the regularity assumption (27), the starting weights ωnj ‘‘compensate ’’ nonzero initial condition u(0) = u0 . This is crucial in making the Lubich Fractional Backward Difference method (and so the proposed scheme (25)) convergent in the case u0 ̸= 0. This is our case because u0 represents the observed noisy image. We refer to [16,12] for the more detailed analysis concerning quoted facts (especially we refer to the prof of the Theorem 4.3.12, Chapter 4, [12]). We recall that the Lyapunov stability of the numerical scheme means that the numerical solution remains bounded for arbitrary bounded perturbations of the initial conditions. In the proof of the next proposition we need the discretized version of Gronwall Lemma (see for example [18]): Let {yn }n∈N0 and {gn }n∈N0 be nonnegative sequences, and c > 0. If ∑  ∑ yn ≤ c + 0≤k 0, and U be the open bounded neighborhood of u0 in which (21) has the unique solution (see Proposition 1). Let un be obtained by (25), with uj = aj , j = 0, . . . , s, for s = card(A) − 1, and u˜ n be obtained by (25), with the perturbed a˜ j = aj + ej , j = 0, . . . , s, for ‖ej ‖ < δ . Then there exist K > 0, so that ‖en ‖ = ‖un − u˜ n ‖ < K δ as n → ∞, for ¯ any h > 0, such that hβ < 1/R, where R > 0 is a Lipschitz constant of G restricted on U. Proof. From (25) it follows that

‖e n ‖ ≤

n −1 −

|ωn−j |‖ej ‖ +

j=0

s −

−β

|ωn,j |‖ej ‖ + hβ R‖en ‖ + hβ

j=0

tn ‖e0 ‖

Γ (1 − β)

.

(31)

¯ with the Lipschitz constat R > 0. Note that as tn = hn, the least term in (31) Here we use the fact that G is Lipschitz on U, equals hβ

h−β n−β ‖e0 ‖ Γ (1−β)

=

(1 − hβ R)zn ≤

n−β ‖e0 ‖ . Γ (1−β)

s −

Now for zn = ‖en ‖ ≥ 0 we have

|ωn,j |zj +

j =0

z0 n−β

Γ (1 − β)

+

n−1 −

|ωn−j |zj .

(32)

j =0

If hβ < 1/R, we obtain 0 < 1/(1 − hβ R) ≤ C1 , for some C1 > 0. Also, as |ωn,j | ≤ C2 n−β−1 for some C2 > 0, and zj ≤ δ, j = 0, . . . , s, the sum of the first and the second term in (32) is bounded above with C3 δ for some C3 > 0 for n ≥ 1. We can now write zn ≤ C δ +

n −1 −

gj zj ,

where gj ≤ C˜ (n − j)−β−1

(33)

j=0

for some C˜ > 0, where C = C1 C3 . It is the consequence of the fact that the binomial coefficients satisfy

z  k

=O



1 Re(z )+1



for z ∈ C, z → ∞, so that ωj = O(j−β−1 ), as those are the coefficients in expansion for (23). Using the discretized version of Gronwall Lemma and the fact that β ≥ 1, from (33) we obtain

 zn ≤ C δ exp

n −1 −

 gj



≤ C δ exp C

j =0

which completes the proof.

n −1 − 1 ˜ j =0

jβ+1

 < K δ,

n ∈ N,

(34)



As the scheme (25) is implicit, in every time step n = 1, . . . , N we need to solve a nonlinear equation in order to obtain un . For that reason we have to prove the structural stability of the proposed method. This means that the numerical solution stays bounded as n → ∞, for an arbitrary bounded perturbation that results from inaccuracies in solving system (25) for un . Proposition 4. Assume that un is an exact solution to the system (25) and that u˜ n is obtained with the defect δn in the discrete time step n ∈ N. Moreover, assume that ‖δn ‖ < δ for all n ∈ N, and that ‖ej ‖ = ‖˜uj − uj ‖ < δ for j = 1, . . . , s. Then, there exist K > 0, so that ‖en ‖ = ‖˜un − un ‖ < K δ , as n → ∞.

736

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

Fig. 1. Original (a, c) and noisy (b, d) test images with PSNR = 15.0.

Proof. From (25) we obtain

‖en ‖ ≤ ‖δn ‖ +

n−1 −

|ωn−j |‖ej ‖ +

j =0

s −

|ωn,j |‖ej ‖.

(35)

j =0

As ‖ej ‖, ‖δn ‖ < δ for j = 1, . . . , s, and ωn,j = O(n−β−1 ), we obtain zn ≤ C δ +

n−1 −

|ωn−j |zj ,

where zn = ‖en ‖

j =0

and by applying discretized version of Gronwall Lemma

 zn ≤ C δ exp

n −1 −

 |ωn−j |

j =0



for some C > 0. Using the fact that ωj = O(j−β−1 ), for some C1 > 0, and that β ≥ 1, we obtain zn ≤ C δ exp C1 K δ, n ∈ N, which completes the proof.

∑∞

1 j=0 jβ+1



<



We note that the function G restricted on U¯ (see Proposition 3), obtained by using (17)–(19) is Lipschitz, for every α > 0. Thus, the fractional order α ∈ [1, 2] of the spatial derivatives that we use in the actual denoising, does not effect either the Lyapunov, or the stability of the proposed method (25).  structural  The matrix aij = [jγ ] , γ ∈ A of system (24) is an exponential Vandermonde matrix with real exponents [12,19]. Consequently it is regular, assuring the unique solution for the starting weights ωn,j . In the general case this system is an

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

737

Fig. 2. Results on the ‘‘Barbara’’ test image, processed by the equations proposed by: Perona (a), Bai (b), Cao (c), and by the novel FFAD (d).

ill-conditioned problem (see [19]). Thus, for the experiments presented in the next section, we select the simplest p = 1 in order to avoid possible numerical pitfalls that could be caused by an inaccurate computation of the starting weights. β With p = 1, from ∑n (23) we get ω(ξ ) = (1 − ξ ) . Now by (24) we obtain s = 0 with just one starting weight ωn,0 = 1/Γ (n − β) − j=0 ωj . Now, from (25) we obtain β

un = h G(un ) −

n −1 −

 ωn−j uj −

j=0

n−β

Γ (n − β)

−

n − j =0

 ωj u0 + hβ

−β

tn u0

Γ (1 − β)

.

(36)

Convolution weights ωj are computed recursively as:

ω0 = 1,



ωj = 1 −

β +1 j



ωj−1 ,

j ∈ N.

(37)

We note that (36) is stable and convergent of order p = 1. Thus, the numerical error given by (26) behaves as O(h), when h → 0. It means that maxn {‖un − u(tn )‖2 } ≤ Kh, for some K > 0. Since it is difficult to find K , we use in experiments sufficiently small time steps h > 0 in order to obtain satisfactory visual results. In implementation, we use the BFGS quasiNewton iterative numerical method for finding the solution un to the nonlinear system (36). This method is applicable to (36) since ∇u G(u) is bounded for all u ∈ U¯ and α > 0. This follows from (13)–(15), (18) and (19) [20]. Moreover, we use the stopping criteria ‖ukn+1 − ukn ‖2 < ε , for some predefined ε > 0. As we have shown in Proposition 4, the method that we proposed is structurally stable. Thus, if we keep the perturbation, in computing solution to (36), small enough (i.e., if we keep ε small enough) we can achieve an arbitrary small computational error. We do this in the experimental results in the next section, so that the computational error could not be visually detected. Note that it would be too expensive in

738

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

Fig. 3. Enlargements for the ‘‘Barbara’’ test image, processed by the equations proposed by: Perona (a), Bai (b), Cao (c), and by the FFAD (d).

terms of required memory and computational time, to hold all samples from t = 0 to t = tn in the computation of un . In order to avoid this, we have used the short memory principle [13,12]. We use the fixed memory length l, so that only the convolution weights ω0 , . . . , ωl−1 are nonzero. In that way (see [13,12]) we introduce an additional error E that is bounded by E ≤ Ml−β /Γ (1 − β), where M = supu∈U¯ G(u) < ∞. 6. Experimental results In this section, experimental results are obtained by applying the proposed FFAD equation on the real images. Images were corrupted by the white Gaussian additive noise. The proposed method was compared with the several baseline denoising methods. One is the classical Perona–Malik (PM) [1] anisotropic diffusion. The second one is the previously mentioned method proposed by Bai and Feng (BF) [10] which gives very good results in preserving edges. The third one is the method proposed by Cao and Yin (CY) [21] which belongs to the class of parabolic-hyperbolic equations applied to image restoration. As well as FFAD, it also uses a partially hyperbolic nature of model in order to preserve edges. In order to understand the experiments, FFAD could be presented as the process where irregular features of the initial image are smoothed by the partially diffusive nature of the equation (for example, set β = 1.5). On the other hand, due to the partially wave-like nature of the equation, the image can be viewed as a regular two dimensional elastic manifold embedded in R3 that contracts and fluctuates, thus generating additional edge enhancement effects. By interpolating between the two behaviors (which is a trade-off), both effects, diffusive and elastic give their contribution to the enhancement process.

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

739

Fig. 4. Results on the ‘‘bicycle’’ test image, processed by the equations proposed by: Perona (a), Bai (b), Cao (c), and by the novel FFAD (d).

As the qualitative, i.e. objective measures of methods efficiency in denoising task, we chose Peak Signal to Noise Ratio (PSNR) measure which is widely spread in the image processing community, as well as the more profound Structural Similarity Index Measure (SSIM), which is actually designed in order to improve the standard PSNR measure (see [22]). For all figures, the optimal stopping time was chosen so that the best SSIM is obtained. On Fig. 1(a) and (c), original test images are presented, while Fig. 1(b) and (d) display noisy images. Fig. 2 presents the results for testing the proposed method in comparison to the previously mentioned baseline methods. The test images were chosen to possess some highly oscillatory regions, for example regions with black and white stripes densely concentrated, as can be seen in the test images. Those are the regions where the FFAD shows better edge enhancement results in comparison to the first [1] and second [10] baseline models which have a pure parabolic nature. It is due to the fact that in those regions contractions and fluctuations caused by its wave-like nature of FFAD come to the full extent. Also, it can be seen that the edge enhancement effect caused by the wave-like nature of the equation is more present in the case of FFAD, than in the case of the equation proposed by Cao and Yin. On Fig. 2(a)–(d) respectively, results from the methods given by Perona and Malik, Bai and Feng, Cao and Yin, and the proposed novel FFAD are presented on the ‘‘Barbara’’ testing image. The order of the spatial fractional derivative present in the FFAD model was set to be α = 1.8 and it is the best reported in [10]. The order of the time fractional derivative was set to be β = 1.5 in order to include diffusion as well as the wave-like behavior of the equation. The memory length of l = 20 with h = 0.1 is chosen. Also, the stopping criterion was determined by setting ε = 0.01. For the model proposed by Cao and Yin, parameter λ was set to 20 which is the best reported in [21], with the time step set to be τ = 0.1. As it can be seen, FFAD has saved the stripes on Barbara’s kerchief and other regions, while the first and the second baseline diffusions have failed to do it. They have over-smoothed the stripes. On the other hand, homogenous regions are equally smoothed (i.e. noise is removed) in the case of the FFAD, as in case of the first two baseline diffusions. The third baseline method showed better

740

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

Fig. 5. Enlargements for the ‘‘bicycle’’ test image, processed by the equations proposed by: Perona (a), Bai (b), Cao (c), and by the FFAD (d). Table 1 Results of denoising of the proposed FFAD method in comparison to the baseline methods are presented, on the ‘‘Barbara’’ test image, for three different values of Gaussian noise standard deviation σ . PSNR

FFAD PM BF CY

SSIM

σ = 15

σ = 20

σ = 25

σ = 15

σ = 20

σ = 25

25.01 24.12 23.98 24.80

22.51 19.83 20.71 22.07

18.31 17.02 17.65 18.10

0.82 0.79 0.80 0.81

0.77 0.73 0.75 0.76

0.70 0.67 0.67 0.68

visual results on the stripes in comparison to the first two methods, but it can be seen that the FFAD makes them even more enhanced. On Fig. 3(a)–(d) close details of the ‘‘Barbara’’ test image are presented for all methods in a similar order. It can be confirmed from the enlargements that the FFAD (Fig. 3(d)) has more significant edge enhancement effect on the highly oscillatory regions than all three baseline methods. On Fig. 4, denoising results are obtained for ‘‘bicycle’’ test image that posses significant amount of highly oscillatory regions (stripes of various thickness in the right upper corner) and the results are presented in a similar order as on Fig. 2. On Fig. 5 the enlargements for the same test image are presented in a similar order. In Table 1, the results of denoising of the proposed FFAD method in comparison to the baseline methods are presented, on the ‘‘Barbara’’ test image, for three different values of Gaussian noise standard deviation σ . It can be seen that FFAD obtained better results for all values of σ , in the means of PSNR as well as SSIM, in comparison to all baseline methods. Optimal stopping time was chosen so that the best PSNR (for first tree columns) or SSIM (for other tree columns) is obtained. In Table 2, similar results are presented for the ‘‘bicycle’’ test image. The same can be concluded, as for the ‘‘Barbara’’ test image.

M. Janev et al. / Mathematical and Computer Modelling 54 (2011) 729–741

741

Table 2 Results of denoising of the proposed FFAD method in comparison to the baseline methods are presented, on the ‘‘bicycle’’ test image, for three different values of Gaussian noise standard deviation σ . PSNR

FFAD PM BF CY

SSIM

σ = 15

σ = 20

σ = 25

σ = 15

σ = 20

σ = 25

25.13 24.21 24.02 24.93

22.63 20.12 21.73 22.17

18.45 17.14 17.70 18.26

0.85 0.81 0.81 0.83

0.79 0.75 0.76 0.78

0.72 0.68 0.69 0.70

7. Conclusion The proposed method is a generalization of the methods proposed by Feng [10] and Cuesta [11] i.e., it interpolates derivatives appearing in the nonlinear PDE in the spatial as well as in the time domain. It thus contains the advantages of both methods. The FODE spatially discretized model is proposed including the FLMM numerical scheme that resolves it. The convergence and the stability of the scheme is also proven. The experiments presented in this work show that with a proper parameter settings, it manages to obtain good results regarding the PSNR and visual effect especially in the highly oscillatory regions. Acknowledgement This research was supported by the Serbian Ministry of Science, grant numbers 174005, 174024, 44003, 32035. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern. Anal. Mach. Intell. 12 (1990) 629–639. F. Catte, P.L. Lions, J.M. Morel, T. Coll, Image selective smoothing and edge detection by nonlinear diffusion, SIAM J. Numer. Anal. 29 (1992) 182–193. L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. C. Vogel, M. Oman, Iterative methods for total variation denoising, SIAM J. Sci. Statisc. Comput. 17 (1996) 227–238. Y.L. You, M. Kaveh, Fourth-order partial differential equation for noise removal, IEEE Trans. Image Process. 9 (2000) 1723–1730. M. Lysaker, A. Lundervold, X.C. Tai, Noise removal using fourth-order partial differential equation with aplication to medical magnetic resonance images in space and time, IEEE Trans. Image Process 12 (2003) 1579–1590. A. Chambole, P.L. Lions, Image recovery via total variation minimization and related problems, Numer. Math. 76 (1997) 167–188. P. Blomgren, P. Mulet, T.F. Chan, C.K. Wong, Total variation image restoration: numerical methods and extensions, in: Proc. Int. Conf. Image Process., vol. 3, 1997, pp. 384–387. T.F. Chan, A. Marquina, P. Mulet, High order total variation-based image restoration, SIAM J. Sci. Comput. 22 (2000) 503–516. J. Bai, X. Chu Feng, Fractional order anisotropic difusion for image denoising, IEEE Trans. Image Process. 16 (2007) 2492–2502. E. Cuesta, J. Finat, Image Processing by means of linear integro-differential equation, in: Proc. Int. Conf. Vizual. Imaging and Immage Process., 10.10.2003, Benalmadena, Spain. M. Weilber, Efficient numerical methods for fractional differential equations and their analytical background, Doctorial Dissertation, 09.06.2005. I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. G. Aubert, P. Kornprobst, Mathematical Problems in Image Processing: PDE’s and the Calculus of Variations, in: Applied Mathematical Sciences, vol. 147, Springer-Verlag, 2002. S. Didas, J. Weickert, B. Burgeth, Properties of higher order nonlinear diffusion filtering, J. Math Imaging Vis. 35 (2009) 208–226. C. Lubich, Fractional linear multistep methods for abel-volterra integral equations of the second kind, Math. Comp. 45 (172) (1985) 463–469. P. Henrici, Discrete variable methods in ordinary differential equations, in: Encyclopedia of Mathematics and its Applications, vol. 34, Camb. Univ. Pr., New York, 1968. Dieudonné, Foundations of Modern Analysis, Academic Press, 1960. K. Diethelm, J.M. Ford, N.J. Ford, M. Weilbeer, Pitfalls in fast numerical solvers for fractional differential equations, J. Comput. Appl. Math. 186 (2006) 482–503. Avriel Mordecai, Nonlinear Programming: Analysis and Methods, Dover Publishing, 2003. Y. Cao, J. Yin, G. Liu, M. Li, A class of nonlinear parabolic-hyperbolic equations applied to image restoration, Nonlinear Analysis. RWA 11 (1) (2010) 253–261. Z. Wang, A.C. Bovik, H.R. Sheikh, E.P. Simoncelli, Image quality assessment: from error to structural similarity, IEEE Trans. Image Process. 13 (4) (2004) 600–612.