## Advanced Theoretical and Applied Studies of Fractional Differential Equations 2013

View this Special IssueResearch Article | Open Access

Yi-Fei Pu, Ji-Liu Zhou, Patrick Siarry, Ni Zhang, Yi-Guang Liu, "Fractional Partial Differential Equation: Fractional Total Variation and Fractional Steepest Descent Approach-Based Multiscale Denoising Model for Texture Image", *Abstract and Applied Analysis*, vol. 2013, Article ID 483791, 19 pages, 2013. https://doi.org/10.1155/2013/483791

# Fractional Partial Differential Equation: Fractional Total Variation and Fractional Steepest Descent Approach-Based Multiscale Denoising Model for Texture Image

**Academic Editor:**Juan J. Trujillo

#### Abstract

The traditional integer-order partial differential equation-based image denoising approaches often blur the edge and complex texture detail; thus, their denoising effects for texture image are not very good. To solve the problem, a fractional partial differential equation-based denoising model for texture image is proposed, which applies a novel mathematical method—fractional calculus to image processing from the view of system evolution. We know from previous studies that fractional-order calculus has some unique properties comparing to integer-order differential calculus that it can nonlinearly enhance complex texture detail during the digital image processing. The goal of the proposed model is to overcome the problems mentioned above by using the properties of fractional differential calculus. It extended traditional integer-order equation to a fractional order and proposed the fractional Green’s formula and the fractional Euler-Lagrange formula for two-dimensional image processing, and then a fractional partial differential equation based denoising model was proposed. The experimental results prove that the abilities of the proposed denoising model to preserve the high-frequency edge and complex texture information are obviously superior to those of traditional integral based algorithms, especially for texture detail rich images.

#### 1. Introduction

Fractional calculus has been an important branch of mathematical analysis over the last 300 years [1–4]; however, it is still little known by many mathematicians and physical scientists in both the domestic and overseas engineering fields. Fractional calculus of the Hausdorff measure is not well established after more than 90 years studies [5, 6], whereas fractional calculus in the Euclidean measure seems more completed. So, Euclidean measure is commonly required in mathematics [5, 6]. In general, fractional calculus in the Euclidean measure extends the integer step to a fractional step. Random variable of physical process in the Euclidean measure can be deemed to be the displacement of particles by random movement; thus, fractional calculus can be used for the analysis and processing of the physical states and processes in principle [7–15]. Fractional calculus has one obvious feature, that is, that most fractional calculus is based on a power function and the rest is based on the addition or production of a certain function and a power function [1–6]. It is possible that this feature indicates some changing law of nature. Scientific research has proved that the fractional-order or dimensional mathematical approach provides the best description for many natural phenomena [16–19]. Fractional calculus in the Euclidean measure has been used in many fields, including diffusion process, viscoelasticity theory, and random fractal dynamics. Methods to apply fractional calculus to modern signal analysis and processing [18–30], especially to digital image processing [31–38], are an emerging branch to study, which has been seldom explored.

Integer-order partial differential equation-based image processing is an important branch in the field of image processing. By exploring the essence of image and image processing, people tend to reconstruct the traditional image processing approaches through strictly mathematical theories, and it will be a great challenge to practical-oriented traditional image processing. Image denoising is a significant research branch of integer-order partial differential equation-based image processing, with two kinds of denoising approach: the nonlinear diffusion-based method and the minimum energy norm-based variational method [39–42]. They have two corresponding models, which are the anisotropic diffusion proposed by Perona and Malik [43] (Perona-Malik or PM) and the total variation model proposed by Rudin et al. [44] (Rudin-Osher-Fatemi or ROF). The PM model simulates the denoising process by a thermal energy diffusion process and the denoising result is the balanced state of thermal diffusion, while the ROF model describes the same thermal energy by a total variation. In further study, some researchers have applied the PM model and the ROF model to color images [45, 46], discussed the selection of the parameters for the models [47–51], and found the optimal stopping point in iteration process [52, 53]. Rudin and his team proposed a variable time step method to solve the Euler-Lagrange equation [44]. Vogel and Oman proposed improving the stability of ROF model by a fixed point iteration approach [54]. Darbon and Sigelle decomposed the original problems into independent optimal Markov random fields by using level set methods and obtained globally optimal solution by reconstruction [55–57]. Wohlberg and Rodriguez proposed to solve the total variation by using an iterative weighted norm to improve the computing efficiency [58]. Meanwhile, Catté et al. proposed to perform a Gaussian smoothing process in the initial stage to improve the suitability of the PM model [59]. However, PM model and ROF model have some obvious defects in image denoising; that is, they can easily lose the contrast and texture details and can produce staircase effects [39, 60, 61]. Some improved models have been proposed to solve these problems. To maintain the contrast and texture details, some scholars have proposed to replace the norm with the norm [62–65], while Osher et al. proposed an iterative regularization method [66]. Gilboa et al. proposed a denoising method using a numerical adaptive fidelity term that can change with the space [67]. Esedoglu and Osher proposed to decompose images using the anisotropic ROF model and retaining certain edge directional information [68]. To remove the staircase effects, Blomgren et al. proposed to extend the total variation denoising model by changing it with gradients [69, 70]. Some scholars introduced high-order derivative to energy norm [71–76]. Lysaker et al. integrated high-order deductive to original ROF model [77, 78], while other scholars proposed a two-stage denoising algorithm, which smoothes the corresponding vector field first and then fits it by using the curve surface [79, 80]. The above methods have provided some improvements in maintaining contrast and texture details and removing the staircase effect, but they still have some drawbacks. First, the improved algorithms have greatly increased calculation complexity for real-time processing and excessive storage and computational requirements will lead them to be impractical. Second, the above algorithms are essentially integer-order differential based algorithm, and thus they may cause the edge field to be somewhat fuzzy and the texture-preservation effect to be less effective than expected.

We therefore propose to introduce a new mathematical method—fractional calculus to the denoising field for texture image and implementing a fractional partial differential equation to solve the above problems by the integer-order partial differential equation-based denoising algorithms [23, 33–38]. Guidotti and Lambers [81] and Bai and Feng [82] have pushed the classic anisotropic diffusion model to the fractional field, extended gradient operator of the energy norm from first-order to fractional-order, and numerically implemented the fractional developmental equation in the frequency domain, which has some effects on image denoising. However, the algorithm still has certain drawbacks. First, it simply took the gradient operator of the energy norm from first order to fractional order and still cannot essentially solve the problem of how to nonlinearly maintain the texture details via the anisotropic diffusion. Therefore, the texture information is not retained well after denoising. Second, the algorithm does not include the effects of the fractional power of the energy norm and the fractional extreme value on nonlinearly maintaining texture details. Third, the method does not deduce the corresponding fractional Euler-Lagrange formula according to fractional calculus features and directly replace it according to the complex conjugate transpose features of the Hilbert adjoint operator. It greatly increased the complex of the numerical implementation of the fractional developmental equation in frequency field. Finally, the transition function of fractional calculus in Fourier transform domain is . Its form looks simple, but the Fourier’s inverse transform of belongs to the first kind of Euler integral, which is difficult to calculate theoretically. The algorithm simply converted the first-order difference into the fractional-order difference in the frequency domain form and replaced the fractional differential operator, which has not solved the computing problem of the Euler integral of the first kind.

The properties of fractional differential are as follows [23, 24, 38]. First, the fractional differential of a constant is non-zero, whereas it must be zero under integer-order differential. Fractional calculus varies from a maximum at a singular leaping point to zero in the smooth areas where the signal is unchanged or not changed greatly; note that, by default, any integer-order differential in a smooth area is approximated to zero, which is the remarkable difference between the fractional differential and integer-order differential. Second, the fractional differential at the starting point of a gradient of a signal phase or slope is nonzero, which nonlinearly enhances the singularity of high-frequency components. With the increasing fractional order, the strengthening of the singularity of high-frequency components is also greater. For example, when , the strengthening is less than when . The integral differential is a special case of the fractional calculus. Finally, the fractional calculus along the slope is neither zero nor constant but is a nonlinear curve, while integer-order differential along slope is the constant. From this discussion, we can see that fractional calculus can nonlinearly enhance the complex texture details during the digital image processing. Fractional calculus can nonlinearly maintain the low-frequency contour features in smooth area to the furthest degree, nonlinearly enhance the high-frequency edge and texture details in those areas where gray level changes frequently, and nonlinearly enhance high-frequency texture details in those areas where gray level does not change obviously [23, 33–38].

A fractional partial differential equation-based denoising algorithm is proposed. The experimental results prove that it can not only preserve the low-frequency contour feature in the smooth area but also nonlinearly maintain the high-frequency edge and texture details both in the areas where gray level did not change obviously or change frequently. As for texture-rich images, the abilities for preserving the high-frequency edge and complex texture details of the proposed fractional based denoising model are obviously superior to the traditional integral based algorithms. The outline of the paper is as follows. First, it introduces three common-used definitions of fractional calculus, that is, Grümwald-Letnikov, Riemann-Liouville, and Caputo, which are the premise of the fractional developmental equation-based model. Second, we obtain fractional Green’s formula for two-dimensional image by extending the classical integer-order to a fractional-order and also fractional Euler-Lagrange formula. On the basis, a fractional partial differential equation is proposed. Finally, we show the denoising capabilities of the proposed model by comparing with Gaussian denoising, fourth-order TV denoising, bilateral filtering denoising, contourlet denoising, wavelet denoising, nonlocal means noise filtering (NLMF) denoising, and fractional-order anisotropic diffusion denoising.

#### 2. Related Work

The common-used definitions of fractional calculus in the Euclidean measure are Grümwald-Letnikov definition, Riemann-Liouville definition, and Caputo definition [1–6].

Grümwald-Letnikov defined that -order differential of signal can be expressed by where the duration of is and is any real number (fraction included). denotes Grümwald-Letnikov defined fractional-order differential operator, and is Gamma function. Equation (1) shows that Grümwald-Letnikov definition in the Euclidean measure extends the step from integer to fractional, and thus it extends the order from integer differential to fractional differential. Grümwald-Letnikov defined fractional calculus is easily calculated, which only relates to the discrete sampling value of that correlates to and irrelates to the derivative or the integral value.

Riemann-Liouville defined the -order integral when is shown as where represents the Riemann-Liouville defined fractional differential operator. As for -order differential when , satisfies . Riemann-Liouville defined -order differential can be given by Fourier transform of the is expressed as where denotes imaginary unit and represents digital frequency. If is causal signal, (4) can be simplified to read

#### 3. Theoretical Analysis for Fractional Partial Differential Equation: Fractional Total Variation and Fractional Steepest Descent Approach Based Multiscale Denoising Model for Texture Image

##### 3.1. The Fractional Green Formula for Two-Dimensional Image

The premise of implementing Euler-Lagrange formula is to obtain the proper Green’s formula [83]. We therefore extend the order of Green’s formula from the integer to a fractional first in order to implement fractional Euler-Lagrange formula of two-dimensional image.

Consider to be simply connected plane region, taking the piecewise smooth curve as a boundary; then the differintegrable functions and [1–6] are continuous in and , and the fractional continuous partial derivatives for and exist. If we consider to represent the first-order differential operator, then represents the -order fractional differential operator, denotes the first-order integral operator, and represents -order fractional integral operator when . Note that represents the -order integral operator of curve surface in the plane. is the -order integral operator in the section of curve along the direction of . is -order fractional integral operator in the closed curve along counter-clockwise direction. Consider that boundary is circled by two curves , , or , , , as shown in Figure 1.

As for differintegrable function [1–6], when , it follows that − . Thus, when − , it has Similarly, it has Fractional Green’s formula of two-dimensional image can be expressed by When and are reciprocal, that is, , it follows that . Then, (8) can be simplified to read We know from (9) the following. First, when , (9) can be simplied as , which is the expression of fractional Green’s formula in reference [84]. Second, when , it follows that . The classical integer-order Green’s formula is the special case of fractional Green formula.

##### 3.2. The Fractional Euler-Lagrange Formula for Two-Dimensional Image

To implement a fractional partial differential equation-based denoising model, we must obtain the fractional Euler-Lagrange formula first, and thus we furtherly deduce the fractional Euler-Lagrange formula for two-dimensional image based on the above fractional Green’s formula.

Consider the differintegrable numerical function in two-dimensional space to be and the differintegrable vector function to be [1–6]; the -order fractional differential operator is . Here, is a type of linear operator. When , then represents an equality operator, which is neither differential nor integral, where and , respectively, represent the unit vectors in the - and -directions. In general, the two-dimensional image region is a rectangular simply connected space, and thus the piecewise smooth boundary is also a rectangular closed curve, as shown in Figure 2.

Referring to (2), it follows that , and . From the fractional Green formula (8) and Figure 2, we can derive Since and , we have Referring to the homogeneous properties of fractional calculus and (10), we can derive that . Then, we have where the sign denotes the inner product. Similar to the definition of fractional-order divergence operator , we consider the -order fractional differential operator to be and the -order fractional divergence operator to be , where both and are the linear operators. In light of Hilbert adjoint operator theory [85] and (12), we can derive where denotes the integral form of the -order fractional inner product. is -order fractional Hilbert adjoint operator of . Then, it follows that where is a fractional Hilbert adjoint operator. When , (12) becomes where denotes the integral form of the first-order inner product, represents the first-order divergence operator, and represents the first-order Hilbert adjoint operator of . As for digital image, we find that Equations (14) and (16) have shown that the first-order Hilbert adjoint operator is the special case of that of fractional order. When , (14) becomes Since the line and line meet at right angle, it has . As for any two-dimensional function , the corresponding and are randomly chosen. According to the fundamental lemma of calculus of variations [83], we know that is required to make (17) established. Since is the positive integer belonging to , we only need to make established. Equation (17) can be rewritten as . When and only when the below equation is satisfied, (17) is established: Equation (18) is the fractional Euler-Lagrange formula, which corresponds to .

Also, if one considers to be the numerical function of vector function and to be the numerical function of differintegrable vector function [1–6], similarly, when , (14) can be written as Equation (19) is established, when and only when Equation (20) is the fractional Euler-Lagrange formula corresponding to .

Since it has no matter what is, we know that the fractional Euler-Lagrange formulas (18) and (20) are irrelevant to the integral order of fractional surface integral . We therefore only adopt the first-order surface integral instead of that of fractional order, when we discuss the energy norm of fractional partial differential equation-based model for texture denoising below.

##### 3.3. The Fractional Partial Differential Equation-Based Denoising Approach for Texture Image

Based on the fractional Euler-Lagrange formula for two-dimensional image, we can implement a fractional partial differential equation-based denoising model for texture image.

represents the gray value of the pixel , where is the image region, that is, . Consider to be the noised image and to represent the desired clean image. Since the noise can be converted to additive noise by log processing when it is multiplicative noise and to additive noise by frequency transform and log processing when it is convolutive noise, we assume that is additive noise, that is, without loss of generality. Consider to represent the additive noise, that is, ; we adopted the fractional extreme to form the energy norm. Similar to the fractional -cover in the Hausdorff measure [96, 97], we consider the fractional total variation of image to be , where is any real number including fractional number and is the hypercube measure. We assume the fractional variation-based fractional total variation as

Consider to be the -order extremal surface of , the test function is the admitting curve surface close to the extremal surface, that is, , we then correlate and merge and by . When , it is the -order extremal surface . Assume that and , where is the cross energy of the noise and clean signal , that is, , and it also the measurement of the similarity between and . We therefore can explain the anisotropic diffusion as energy dissipation process for solving the -order minimum of fractional energy norm , that is, the process for solving minimum of is to obtain the minimum similarity between the noise and the clean signal. Here, plays the role of nonlinear fidelity during denoising, and is regularized parameter. Fractional total variation-based fractional energy norm in surface family can be expressed by

As for , if -order fractional derivative of exists, has the -order fractional minimum or stationary point when . Referring to the linear properties of fractional differential operator, we have where it has as for the vector , and the signal denotes the inner product. Unlike the traditional first-order variation, (23) is the -order fractional extreme of , which aims to nonlinearly preserve the complex texture details as much as possible when denoising by using the special properties of fractional calculus that it can nonlinearly maintain the low-frequency contour feature in the smooth area to the furthest degree and nonlinearly enhance the high-frequency edge information in those areas where gray level changes frequently and also nonlinearly enhance the high-frequency texture details in those areas where gray level does not change obviously [33–38].

Provided that is a fractional number, when , it has . Referring to (11) and Faà de Bruno formula [95], we can derive the rule of fractional calculus of composite function as where and is the constant. is separated from summation item. From (24), we know that the fractional derivative of composite function is the summation of infinite items. Here, satisfies The third signal in (25) denotes the summation of of the combination of that satisfied (25). Recalling (23), (24) and the property of Gamma function, we can derive Without loss of the generality, we consider for simple calculation; it then has and ; thus (26) can be reduced as We know from (25), when , satisfies and ; then it can be furtherly reduced as When takes odd number (, ) and even number (, ), respectively, the expressions of are also different: Put (29) into (28), then it becomes We consider that . Since the corresponding is randomly chosen for any two-dimensional numerical function , − is random. Referring to (20), the corresponding fractional Euler-Lagrange formula of (30) is where it has . It is the corresponding fractional Euler-Lagrange formula of (23).

As for , it has Since the test function is randomly chosen, is also random. According to the fundamental lemma of calculus of variations [83], we know that to make (32) established, it must have Equations (23) and (32), respectively, are the -order minimal value of and . Thus, when we take , 2, and 3, the -order minimal value of (22) can be expressed by where it has . is calculated by the approach of fractional difference. We must compute . If image noise is the white noise, it has . When , it converges to a stable state. Then, we merely multiply at both sides of (34) and integrate by parts over and the left side of (34) vanishes:

Here, we simply note the fractional partial differential equation-based denoising model as FDM (a fractional developmental mathematics-based approach for texture image denoising). When numerical iterating, we need to perform low-pass filtering to completely remove the faint noise in very low-frequency and direct current. We know from (34) and (35) that FDM enhances the nonlinear regulation effects of order by continually multiplying function and power of and enhance the nonlinear regulation effects of order by increasing in the denominator. Also, we know from (34) that FDM is the traditional potential equation or elliptic equation when , the traditional heat conduction equation or parabolic equation when , and the traditional wave equation or hyperbolic equation when . FDM is the continuous interpolation of traditional potential equation and heat conduction equation when and the continuous interpolation of traditional heat conduction equation and wave equation when . FDM has pushed the traditional integer-order partial differential equation-based image processing approach from the anisotropic diffusion of integer-order heat conduction equation to that of fractional partial differential equation in the mathematical and physical sense.

#### 4. Experiments and Theoretical Analyzing

##### 4.1. Numerical Implementation of the Fractional Partial Differential Equation-Based Denoising Model for Texture Image

We know from (34) and (35) that we should obtain the fractional differential operator of two-dimensional digital image before implementing FDM. As for Grümwald-Letnikov definition of fractional calculus in (1), it may remove the limit symbol when is large enough; we then introduce the signal value at nonnode to the definition improving the convergence rate and accuracy, that is, . Using Lagrange 3-point interpolation equation to perform fractional interpolation when , we can obtain the fractional differential operators of YiFeiPU-2, respectively, on the eight symmetric directions [35, 38], which is shown in Figure 3.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

In Figure 3, is the coefficient of operator at noncausal pixel, and is the coefficient of operator at pixel of interest . is often taken odd number. When , we can implement fractional differential operator of , whose coefficients of the operator, respectively, are , , ///, ////////, ////////, ////, and / [35, 38]. Spatial filtering of convolution operator is adopted as the numerical algorithm of fractional differential operator. We take the maximum value of fractional partial differential on the above 8 directions as the pixel’s fractional calculus. As for fractional differential operator in (34) and (35), when , , , , , , , can perform numerical computing referring to Figure 3 and their concerning coefficients [35, 38], and when , we take a special difference as the approximate of first-order differential to maintain calculation stability, that is, + + , + + + .

If time equal interval is , time is noted as when . represents the initial time. The digital image at time is . Thus, we have , where is the original image and is the desired clean image that is a constant, that is, . We can derive the numerical implementation of fractional calculus about time from (1) as where and are the factors for ensuring iterative convergence. Also, the desired clean image is unknown, but the intermediate denoising results is an approximate to the desired clean image, that is, . We consider that for making close to as much as possible. We should compute the fractional calculus on the 8 directions simultaneously in practice in order to improve the calculation accuracy and antirotation capability.

The best image in (36) is unknown in numerical iteration, but the intermediate result is an approximate to , that is, . We consider that to make the iterative result approximate as much as possible. Take in (36) and in (34) and (35), then we can derive the numerical implementation of (34) and (35) as where , , //.

We should pay attention to the following when performing numerical iterative implementation. First, is a small number in (37) to ensure convergence, and here we take . Second, we do not need to know or estimate the variance of noise, but we need to assume to be a small positive number in the first iteration. We therefore assume that in the experiment below. We take to (38) and perform numerical iteration. Each iterative result may be different, but it is the approximate to the variance of noise. Third, it is possible that in numerical iterative calculation, so we take when to ensure that (37) and (38) are meaningful. Fourth, we take when to make have meaning. Fifth, to completely denoise faint noise in very low-frequency and direct current, FDM takes the simple way by reducing the convex in the area where gradient is not changed obviously. We therefore need to perform low-pass filtering for very low-frequency and direct current in numerical iterating. The practices of (37) and (38) are as follows. For one-dimensional signal, we consider that when and and if the noise is not severe in order to ensuring denoising effect, and in the rest conditions we consider that , where , , and if the noise is very strong, we consider that to accelerate the denoising speed. For two-dimensional signal, we consider . When , they are and , when they are and , when they are and , and in the rest occasions, it has