A fractional-order derivative based variational framework for image denoising

In this paper, we propose a unified variational framework for noise 
removal, which uses a combination of different orders of fractional 
derivatives in the regularization term of the objective function. 
The principle of the combination is taking the order two or higher 
derivatives for smoothing the homogeneous regions, and a fractional 
order less than or equal to one to smooth the locations near the 
edges. We also introduce a novel edge detector to better detect 
edges and textures. A main advantage of this framework is the 
superiority in dealing with textures and repetitive structures as 
well as eliminating the staircase effect. To effectively solve the 
proposed model, we extend the first-order primal dual algorithm to 
minimize a functional involving fractional-order derivatives. A set 
of experiments demonstrates that the proposed method is able to 
avoid the staircase effect and preserve accurately edges and 
structural details of the image while removing the noise.


1.
Introduction. Image denoising is a fundamental problem in image processing and computer vision. It aims at smoothing a noisy image without losing significant features. The difficulty of the problem is in the preservation of image structures and details while removing the noise. Although many methods have been developed to tackle this problem, developing novel methods for better preserving image edges and fine structures, especially for weak edges or textures, is still very important. In this work, we will present a hybrid model by integrating different fractionalorder derivatives with a spatially adaptive weight as the regularization terms for smoothing noisy images.
In this paper, we are interested in restoring an image u from an observed image f corrupted by some additive noise n. Such kind of image degradation process could be mathematically modeled by (1) f = u + n, which is the so-called additive noise model. It is one of the most common image noise models. In the last two decades, variational methods have achieved great success in reducing such noise owing to the use of total variation (TV). For instance, the TV based denoising model, originally introduced by Rudin, Osher and Fetami (ROF) in their pioneering work [19], has been widely applied in various image denoising and reconstruction problems. The unconstrained version of ROF model introduced in [19] reads: where Ω ⊂ R 2 is the image domain, λ > 0 is a parameter balancing the weight of smoothing and data fidelity, and Ω |Du| is the BV semi-norm or TV norm, which is defined as where p L ∞ (Ω) = max{ p 2 1 (x, y) + p 2 2 (x, y)}. A significant advantage of the TV regularization is that it leads to an underlying sparse solution (in the finite difference domain), and allows for sharp discontinuities in the solution for feature preserving. However, both from a theoretical ( [23,17]) and an experimental ( [12,2,5]) point of view, it suffers from the staircase effect, that is to say, the transformation of smooth regions into piecewise constant ones, which may produce an undesirable blocky image. Later, high order partial differential equations (PDEs), especially the fourth order PDE-based methods (see for instance [24,12,5]), were proposed to solve this problem in view of their stronger smoothing ability. But this kind of methods smoothes the image too much to preserve discontinuities. Therefore, the combined models with the TV filter and the fourth-order filter were put forward in [9,2,3,4,13]. As a result, near the image edges, the TV filter plays a dominating role so that the edges are preserved; and in homogeneous regions, the fourth-order filter plays a major role so that the noise can be effectively removed without producing the staircase effect.
Recently, the fractional-order derivative has been paid much more attention in engineering. The fractional-order derivative can be seen as the generalization of the integer-order derivative. It has been studied by many mathematicians (for example Euler, Hardy, Littlewood, and Liouville) [18] for many years. There are several ways to define the fractional-order derivatives. Among them the Grünwald-Letnikov (G-L) fractional-order derivative, Cauchy-integral fractional-order derivative (or Reimann-Liouville (R-L) fractional-order derivative), Caputo fractional-order derivative, and the definition in Fourier domain (frequency domain) [11] are mostly used. But not until Mandelbrot found fractals and applied the R-L fractional-order derivative to the Brownian motion did the fractional-order derivative cause great attention. Up to now, it has been applied to many fields such as noise detection and estimation [10,14], electromagnetic theory [6], wavelets, splines [21,20] and many others.
In very recent years, the applications of fractional order derivatives to image restoration have had a fast development due to its ability to balance the noise removal and edges and textures preservation in images. In [25], Zhang and Wei extended the ROF model to the fractional-order case. They suggested that the fractional order in (1,2) is better for preserving fine scale features such as textures and does not give rise to staircase effect. But the sharp edges will be blurred to a certain degree. Bai and Feng [1] introduced a class of fractional-order anisotropic diffusion models for noise removal under the variational framework, which, in fact, is a generalization of the Perona-Malik model [15]. Their experimental results illustrated that the fractional order greater than one is better than the integer orders. Besides, Guidotti and Longo [7] modified an existing fourth-order PDE model using a fractional Laplacian operator and a fractional gradient, respectively, as edge detectors. Their two fractional high-order diffusion models can eliminate the undesirable staircase effect without giving rise to more blurs when removing the noise. Later, Zhang et al. [27,26] proposed a fractional-order multi-scale variational model by replacing the first-order derivative with the fractional-order derivative in the regularization term, and substituting a kind of multiscale norm in negative Sobolev space for the L2 norm in the fidelity term of the ROF model. The model can better preserve the textural information and eliminate the staircase effect. Various experimental results on the use of fractional-order derivatives indicate that once the appropriate order is chosen, the staircase effect can be avoided, but the use of a single order derivative for smoothing may not well perform denoising, staircase reduction, and sharp edge and texture preservation simultaneously. The order less than one can not remove the noise effectively; and the order greater than one and less than two can not preserve the sharp edges well, although it is able to avoid the staircase effect.
In this paper, based on the fractional-order derivative theory, we propose a novel hybrid model by integrating different fractional-order derivatives in regularizing an image at different regions. In homogeneous regions, we make use of the fractional order equal to or greater than two to reduce the noise and avoid the staircase effect; and near edges, we utilize the fractional order that lies in (0, 1] to preserve them. The homogeneous regions and the edges will be detected by a new edge-texture detector function, which involves two norms of fractional-order derivatives. The orders are less than one, and greater than two, respectively. The norm of the gradient with the fractional order less than one favors identifying weak edges or textures, while the norm of the gradient with the fractional order greater than two helps to detect strong edges or textures. So unlike the traditional edge detector function just associated with the image gradient, this new defined detector function can detect not only the strong edges or textures but also the weak ones. Finally note that the orders of the derivative in the regularization terms may be chosen to be integer or fractional. By our experiments in Section 6, we find that the orders in the regularization terms should be integer for piecewise smooth images and be fractional for textural images. So the proposed model is more flexible to tackle different kind of images.
Moreover, we will adopt the popular G-L fractional order derivative to implement, which is considered as a generalization of the classical finite difference of the integer orders. Besides, in order to effectively solve the proposed model, we extend the first-order primal dual algorithm for solving a convex minimization problem with the integer-order derivatives based regularization to the case with the fractionalorder derivatives.
The rest of this paper is organized as follows. In Section 2, we review some closely related works on reducing the staircase effect, including the integer-order and fractional-order variational methods for image denoising. In Section 3, we analyze the filtering character of fractional-order derivatives based on the definition from the aspect of the Fourier transform. In Section 4, we present a novel variational framework for noise removal based on the property of fractional-order derivatives. The numerical algorithm, an extension of the conventional first-order primal dual method, is introduced in Section 5 to solve the proposed model. In Section 6, we present numerous experiments and comparisons with the closely related existing methods to show the effectiveness of the proposed method in preserving edges and textures as well as avoiding the staircase effect while removing the noise. Finally, a brief conclusion is summarized in Section 7.
2. Related works. In this section, we will briefly review some related models for image denoising: the high order PDE-based methods, typically, the LLT model [12], the models combining TV filter with fourth-order filter [2,3,4,13,9], and a fractional-order derivative based variational model [25]. Those models are effective in tackling the problem of smoothing without producing the staircase effect.
2.1. High order PDE-based methods. High order PDE-based methods are one of the effective ways to avoid the staircase effect because of its stronger smoothing ability than the second-order PDEs. In [24], You and Kaveh proposed a class of fourth order PDEs to denoise, meanwhile eliminate the staircase effect, which are obtained by minimizing the following functional: where ∆ denotes the Laplacian operator, and the function ϕ is an increasing function. The associated Euler-Lagrange equation is a fourth-order PDE.
In [12], Lysaker, Lundervold and Tai (LLT) also proposed a fourth order PDEbased denoising model, and showed the effectiveness in solving the staircase effect problem. The model reads as where |∇ 2 u| = u 2 xx + u 2 yy + u 2 xy + u 2 yx , and λ is a positive parameter. However, the fourth order PDE-based methods suffer from the blurring effect near sharp edges due to the possible over-smoothing. Later, to tackle this problem, the combination of TV filter and fourth-order filter was proposed to respectively preserve edges and alleviate the staircase effect.

2.2.
Combined models with TV filter and fourth-order filter. Along this idea, in [2], the following model was proposed to perform different types of smoothing in different regions: In this method, the restored image is decomposed into two parts. The discontinuous component of the image u represented by u 1 is regularized by TV norm to preserve the discontinuities; while regions of moderate slopes are allotted to the u 2 component, which is regularized by the second order derivative, so the staircase effect is eliminated. Later, a variant of the model (5) is introduced by Chan, Esedoglu, and Park (CEP) [3] to fast solve the u 2 component with Ω |∇(∇u 2 )| replaced by In order to keep the discontinues of the image, the authors also improved the model (6) by replacing Ω |∆u 2 | 2 with Ω |∆u 2 |, i.e., the CEP2-L2 model [4], On the other hand, as mentioned in subsection 2.1, although the LLT model does not give rise to the staircase effect, the sharp edges may be blurred. In [13], Lysaker and Tai improved the LLT model by combining it with the TV regularization. By alternatively implementing the TV filter to get the solution u and the LLT filter to get the solution v, they generated a new solution by taking the best from each of the two methods, which is achieved by a convex combination of two solutions w = θu + (1 − θ)v , where the weighting function θ ∈ [0, 1] can be found adaptively such that the solution is u along edges and is v in smooth regions. Similar to the idea of combining the LLT model with TV model, Li et al. [9] proposed a new integrated model in the variational framework by a spatially adaptive function to determine when the TV filter or the fourth order filter takes effect. Their variational model is given by where g(x, y) is an edge detector function associated with the gradient of the noisy image which is chosen as and γ > 0 is a small constant to guarantee g < 1, G a is the Gaussian kernel with the standard deviation a. The authors also proved the existence and uniqueness of the minimizer. Their experimental results showed that in homogeneous regions the staircase effect is avoided because the fourth order filter dominates the smoothing, and near the edges the discontinuity are preserved because the TV filter dominates the smoothing.
In this work, we will employ the idea of performing different types of smoothing at different location to propose a novel variational model, based on fractional-order derivatives and a new edge and texture detector.

2.3.
Fractional-order derivative based variational model. The fractionalorder derivative models under the variational framework have been presented in several references such as [8,27,25]. We will review the one [25] that can be considered to be the generalization of the TV based denoising model, and closely related to our proposed model. It is described as follows where D α u = (D α x u, D α y u) denotes the gradient with order α > 0 that may be fractional or integer number, and λ > 0 is a parameter balancing the contributions between the regularization term and the fidelity term.
Note that when α = 1, the model is precisely the ROF model; when α = 2, the regularization term becomes Ω |D 2 u|, where D 2 u = (D 2 x u, D 2 y u) = (D xx u, D yy u) and |D 2 u| = (D xx u) 2 + (D yy u) 2 , which can be considered to be a simplified version of LLT model without the second-order mixed derivatives. Thus, by choosing different orders of derivatives in (9), the preservation of edges and textures or the elimination of the staircase effect can be achieved when removing the noise.
3. Filtering character of the fractional-order derivative. In order to analyze the filtering character of fractional-order derivatives, we first introduce a definition based on the Fourier transform, which is considered as an extension of the notion of an integer-order derivative using the Fourier transform. More other definitions can be found from [16,11].
3.1. Definition of the fractional-order derivative in Fourier domain. For a real number α ≥ 0, the αth order derivative of a one-dimensional function f : R → R in Fourier domain is defined as follows where i = √ −1 and F (w) is the Fourier transform of f . This definition can be viewed as a direct extension of the following property for an integer-order derivative: Then for any two dimensional function f (x), x = (x 1 , x 2 ) ∈ R 2 , the αth order partial derivative of f with respect to x j is defined as follows where w = (w 1 , w 2 ) and < w, x >= w 1 x 1 + w 2 x 2 .

3.2.
Filtering character of the fractional-order derivative. To explore the idea of how to apply fractional-order derivatives for better image denoising with edge and texture preserving, in this subsection we will illustrate the filter character of fractional-order derivatives. To this end, we will resort to the Fourier transform of the αth-order derivative. For simplicity, we use the definition (10) for the one-dimensional function, and the same properties can be obtained for the two-dimensional function. By (10), we have where f is the Fourier transform of f , i.e., F (w). From (10) and (13), we can view the fractional-order derivative operator as a filter in frequency domain, and the frequency response is (iw) α . More specifically, for α > 0, D α can be viewed as a high-pass filter, i.e., the high frequency components of an image can pass by this operator. In fact, the amplitude of D α f in frequency domain is (14) | that is to say, in frequency domain, the amplitude of f is changed by multiplying |(iw) α | through the operator D α . So by analyzing the function |(iw) α |, we can understand the filtering character of the fractional-order derivative D α . Let The curve of H(w) with different α is shown in Fig.1. From the Fig.1, we can see that when α > 0, the high frequency (|w| > 1) is enhanced, and the low frequency (|w| < 1 ) is restrained. Hence, the fractional-order derivative operators can be viewed as high-pass filters. Moreover, we observed that the high frequency components of f are relatively less enhanced for 0 < α < 1 than that for α > 1, but the low-frequency components of f are restrained more by applying the fractional-order derivative D α with α > 1. Based on these observations, by minimizing |D α u| in the model (9), the edges corresponding to the high frequency components can be better preserved if the fractional order α is less than or equal to one, while the homogeneous regions corresponding to the low frequency can be better smoothed, if the fractional order α is greater than one.
Besides, from Fig. 1, we observe that for 0 < α < 1 the low frequency of f are not restrained more than that for α > 1, and for α > 1 the high frequency of f are enhanced more than that for 0 < α < 1. Generally, the edges and textures of an image have different frequencies. If the textures have higher frequency than edges, then the derivatives with order α 2 > 1 would be better to detect textures, while the weakly high frequency will be restrained, so derivatives with order 0 < α 1 ≤ 1 would be required to detect edges. On the contrary, if the edges have higher frequency than textures, the derivatives with order α 2 > 1 would be better to detect edges, and derivatives with order 0 < α 1 ≤ 1 would be required to detect textures. By combining those two different order derivatives into the edge detector, the ability of detecting textures and edges can be enhanced. Therefore, we propose a new edge detector for better edge and texture detection, which uses the derivative with order larger than one to detect higher frequency components, and the derivative with order less than one to identify weakly high frequency components.  4. Our proposed model. In this section, we present a new variational model for selective smoothing with edge and texture preserving. In this model the regularization is carried out by minimizing the magnitudes of two derivatives with different fractional (or integer) orders, aiming at performing different type of smoothing at different regions. Those regions (homogeneous, texture or edges) are identified by a new edge detector, which also involves a combination of different fractional-order derivatives.
For any real number α > 0, denote D α u (D α x u, D α y u). The proposed variational model can be described as follows: where the orders of the fractional derivatives α > 1, β ∈ (0, 1], µ > 0 and λ > 0 are parameters balancing the contribution of the three terms in the objective function. Here g(x, y) is an edge detector defined using the fractional order derivatives and given by (16) g(x, y) where α 1 ∈ (0, 1], α 2 > 1, k 1 and k 2 are two positive constants, and G a is a Gaussian kernel with the standard deviation a. The edge detector defined in [7] only uses the magnitude of fractional order derivative |D α1 u| (0 < α 1 ≤ 1), which is a special case of (16) with k 2 = 0. Although their experimental results showed the effectiveness of this edge detector in edge-preserved smoothing, we observed that by adding |D α2 u| with α 2 ≥ 2 the edge detector can better identify the strong edges or textures. Since |D α2 u| with α 2 ≥ 2 is larger in the regions corresponding to the high frequency components (i.e. strong edges, textures) than that corresponding to relatively less higher frequency components (i.e. weak edges, textures). By using a combination of |D α2 u| with α 2 ≥ 2 and |D α1 u| with 0 < α 1 ≤ 1, the detector is able to identify not only strong edges or textures but also weak ones. Consequently, the proposed model can have better performance in preserving the edges and textures.
In the following, we will analyze the denoising effect of the proposed model. In homogeneous regions, the magnitudes of the fractional-order gradients |D α1 u| and |D α2 u| are all small, thus the edge-texture detector function g is close to one, so that the regularization is dominated by the first term in (15), where the order of the derivative is greater than one to avoid the staircase effect. Near the location of edges and textures, the magnitude of the fractional order derivatives |D α1 u| and |D α2 u| are all large, then the edge detector function g takes smaller value (close to zero). Hence, the regularization in the model is dominated by the second term, which favors the preservation of the edges and textures in view of the derivatives with order β ∈ (0, 1]. However, the amount of smoothing near strong or weak edges and textures could be different due to the different value of g there. The proposed model (15) generalized several existing models by choosing various values of α, β, and g. For instance, when g = 0 or g = 1, the model (15) reduces to the fractional-order derivative based variational model (9). When g = 0 and β = 1, the model (15) is exactly the ROF/TV denoising model (2). When g = 1 and α = 2, the model becomes a second-order variational based model similar to the LLT model without the term of D xy u in D 2 u. Furthermore, if α = 2, β = 1, α 1 = 1 and k 2 = 0 in g, the model (15) is similar to the model (8) with the little difference in the second-order regularization term. 5. Numerical algorithm.

5.1.
Discrete forms of the fractional-order operators. In this subsection, we introduce a discrete form of the fractional order derivative proposed by Grünwald and Letnikov. It is an extension of the discrete form (using finite difference) of the integral-order derivatives. It is more convenient for numerical computation.
For a real function u : Ω → R, where Ω ⊆ R 2 is a bounded open set, as given in section 2.3, the fractional-order derivative is represented as (17) D α u = (D α x u, D α y u), α ∈ R + The G-L discrete form of the fractional order derivatives D α x and D α y introduced by Grünwald-Letnikov is given as follows where C k α is the generalized binomial coefficient, i.e., , and the gamma function Γ(·) is given by the following formula Let D α * x and D α * y be the adjoint operators of D α x and D α y , then they are defined respectively by For a vector function p(x, y) = (p 1 (x, y), p 2 (x, y)), we give the fractional divergence as follows Then the adjoint operators and the fractional divergence operators have the following relationship Remark 1. According to the above notations, when α = 1, C k 1 = 0 for k > 1, then the gradient and divergence operators become which are consistent with the classical backward and forward differences.

Remark 2.
In general, the parameter K in the G-L fractional-order derivatives (18) is chosen to be larger than 3, so the fractional-order derivative is related with more pixels than the integer-order derivative, which is related with only its two or four neighbor pixels. Thus, the fractional-order derivative will be more favorable to structure preservation in image restoration.

5.2.
Algorithm. In this subsection, we provide the corresponding algorithm for solving the proposed model, which is based on the dual formulation of the norm of the fractional-order derivatives.
First, we derive the dual form for the functional where h(x, y) is a positive valued function in R 2 .
To discretize the equations (30)-(32), we assume that the discrete image u is of N × N pixels, and let u ij = u(i, j) for i, j = 1, · · · , N . Thus we have In the numerical experiments, we choose K = 15.
Then the discrete versions of the equations (30)-(32) are represented by To automatically determine the parameter λ, we adopt the following strategy. For the saddle point problem (26), fixed p and q, we have the minimization problem So the stationary point can be obtained by solving By multiplying u − f in (37) and integrating it over Ω, we have where and σ is the standard deviation of image noise. Consequently, we can update λ by Finally, we update g(x, y) by As a result, the primal dual algorithm for solving our proposed model can be summarized as follows: Step1. Update p n+1 by (33); Step2. Update q n+1 by (34); Step3. Update u n+1 by (35); Step4. Update λ n+1 by (39); Step5. Update g by (40) for every 50 iterations; Step6. Check whether the stop criterion is satisfied. If not, return to step 1.
Two stopping criterions are given to terminate the algorithm. The first one is that if the condition (41) u n+1 − u n 2 / u n 2 < ε is satisfied, then end the algorithm. The other criterion is that the number of iterations does not exceed a prefixed iterations M . When either criterion is satisfied, we terminate the experiment. In our experiments, we choose ε = 1e − 8 and M = 500.
Remark 3. If the standard deviation of the noise σ is unknown, we can choose an appropriate parameter λ by the trial and error method as usual to obtain the desirable result.
For the quantitative comparison, we use the measures including the peak signal to noise ratio (PSNR), the signal-to-noise ratio (SNR), the mean square error (MSE) and structural similarity (SSIM) metric [22] of the restored image to evaluate the restoration results. They are defined by P SN R = 10 log 10 where u * , u and N 2 are the restored image, the true image and the size of the image, respectively.ū andū * represent the mean values of u and u * in image domain Ω, σ u and σ u * denote the variances of u and u * , σ uu * is the covariance of u and u * , and C 1 and C 2 are the constants chosen considering the human visual system perception and to ensure numerical stability of the division.
In all experiments, we fix τ = 0.1, γ = 1, k 1 = 0.002, k 2 = 0.02 and the Gaussian kernel G a with zero mean and variance a = 0.8. We choose the orders α = 2 and β = 1 for piecewise smooth images and fractional orders for textural images. For the integer-order derivatives based models, we still use the proposed numerical algorithm, because when the orders are integer, the G-L fractional-order derivative becomes the classical finite difference. Other parameters such as α 1 , α 2 and µ need to be chosen according to the observed images. In our experience, we observed that α 1 is better to be less than 1 but close to 1, and α 2 is better to be greater than 2 but close to 2. The parameter µ should be smaller for images with much textures than those with little textures. For each experiment, we will give the specific values for those parameters.   The aim of the first experiment is to preliminarily examine the performance of the proposed model in image restoration by comparison with the following models: ROF, Second-order, CEP2-L2, TV-LLT. In this experiment, we applied those models to the noisy images shown in Fig.2(b) and Fig.3(b), respectively. Both noisy images are generated by adding Gaussian noise with zero mean and the standard deviation σ = 20 to the original clean images shown in Fig.2(a) and Fig.3(a), respectively. We choose α 1 = 0.8, α 2 = 2.2 and µ = 1/5 for the proposed model. As shown in Figs. 2(c)-2(l) and 3(c)-3(n), the images reconstructed from the proposed model ( Fig. 2(k), Fig. 3(k) and 3(m)) are visually better than those obtained from the ROF, Second-order, CEP2-L2, TV-LLT models. On one hand, comparing with the ROF model, the proposed model avoided the staircase effect; on the other hand, the proposed model reduced blurring on edges that may occur due to the over-smoothing by the second-order derivative based models, and preserved more details than CEP2-L2, TV-LLT models. The ability of the edge preservation can be further seen from the difference images f − u * . In all difference images, the residual structures are less by the proposed method. Moreover, from Tables 1, we see that the proposed method has higher PSNR, SNR and SSIM and less MSE than other models in comparison. Those results initially illustrate that the proposed method has better performance in image restoration.
To further show the advantages of the proposed model, we compare it with some first and second order derivatives combined regularization based models, including CEP2-L2 model and TV-LLT model. In addition, by comparing M15 21Reg with g 1 , M15 21Reg with g 2 , M15 21Reg with g and M15 FracReg with g, we also evaluate the contribution of the fractional-order derivative based g and regularization to the efficiency of the denoising process. The Lena image, which is a piecewise smooth image, is used for comparing the performance of the above models. The denoising results are displayed in Fig. 4. Firstly, from the difference image f − u * , the CEP2-L2 model has more residual structures (Fig. 4(d)) than other models. Secondly, comparing the TV-LLT model (Fig. 4(e)) with the proposed model with g = g 1 (Fig. 4(g)), we find that the TV-LLT model showed some over-smoothing, especially for the fringes of the hat, which indicates that the proposed second-order derivative based regularization favors preserving image details. Moreover, from the Table 3, we observe that the use of the new defined edge detector g in the model (15) with α 1 = 1, α 2 = 2, i.e., g 2 (shown in the 4th column) results better PSNR, SNR, MSE and SSIM than the conventional one, i.e. g 1 = 1 1 + k 1 |∇G a * u| 2 (shown in the 3rd column). This demonstrates that by incorporating the higher-order derivative into the edge detector, the edge and texture are better preserved. Although the restored images by M15 21Reg with g and M15 FracReg with g look similar (see Fig. 4(k) and 4(m)), the last column of Table 3 shows the best PSNR, SNR, MSE and SSIM than all other columns. This is the result obtained by using the fractional-order derivatives in the new edge detector with α 1 = 0.8, α 2 = 2.2 and the regularization with α = 2.1, β = 1.0. So the new fractional-order derivative based edge-texture detector function g and regularization are better than other integer ones.  The next experiment is to validate the effectiveness of the proposed framework in restoring textured images. The Barbara image is used in this experiment. Firstly, by comparing integer order based edge detectors with the proposed one, we observe the ability of the proposed fractional-order based edge detector g in discriminating the homogeneous region from edges and textures. Fig. 5(c)-5(f) shows the magnitudes of the integer and fractional order derivatives of the Barbara image. Comparing Fig. 5(c) and 5(d), it shows that the magnitude of fractional order less than 1 (Fig.5(d)) can detect more edge and texture structures, especially the weak edges and textures, than that of the first order one (Fig.5(c)). From Fig. 5(e) and 5(f), the magnitudes of fractional order larger than 2 and the second order do not show big difference visually, however, in the textural regions the value of |D 2.3 (G a * u)| is larger than |D 2 (G a * u)|, which will make the edge detector function (16) smaller so that the weight of the second regularization term is large, and then the textures are preserved better. Fig. 5(g)-5(i) display the images of edge detectors g 1 , g 2 and g, from which, we see that three edge detectors can discriminate the homogeneous region from the textures and edges, but the proposed edge detector g contains more textures and edges (shown as low intensities) than the classic edge detector g 1 and integer orders based g 2 . It illuminates that the proposed edge detector can better discriminate different regions. Secondly, we will further compare the restored results of some combined regularization based models. To show the effectiveness of the fractional-orders based regularization and g for textured images, we also examine the performance of M15 21Reg with g 1 , M15 21Reg with g 2 , M15 21Reg with g and M15 FracReg with g on the Barbara image. The denoising results are displayed in Fig. 6. Comparing with the images resulted from CEP2-L2, TV-LLT models shown in Fig. 6(a) and 6(c), the proposed model can better avoid the staircase effect in the smooth regions such as the hand region. We also compare the adaptive fractional-order multi-scale method in [27] with the proposed one. In [27], the idea of adaptively choosing the differential order α in the regularization term, such that different regularization plays a role in different regions, is very nice. However, the restored result ( Fig. 6(e) still has the slightly staircase effect, especially, in the face region. This indicates that their differential order in the non-textured region limited from 1.0 to 1.2 might be too small to avoid the staircase effect. In addition, from the difference images f − u * , we see that the proposed model has less residual, so it can better preserve image structures, such as edges and textures. In order to have a better visual comparison, we also displayed two corresponding zoom-in areas in Fig. 7. The first zoom-in area contains both smooth parts and textural parts of the image. From Fig. 7(f1)-7(i1), we can see that M15 21Reg with g = g 1 is less effective in texture preservation than M15 21Reg with g = g 2 and the proposed model with fractional-order derivatives based g. The second zoom-in part of all restored images in Fig.7(c2)-7(i2) validated the same conclusion as the first zoom-in area. Besides, from the aspect of preserving textures in the testing image, it may not be easy to visually compare the results of those models from Fig. 7(g1)-7(i1) and 7(g2)-7(i2). But the Table 4 shows that the PSNR, SNR and SSIM are higher and the MSE is smaller by using M15 FracReg with g. It indicates that the fractional-order derivatives based regularization is more suitable for textured images.  Finally, in order to exam the influence of fractional orders α and β in the regularization terms on the results of image recovery, we select the noisy Lena image (Fig. 4(b)) and Barbara image (Fig. 5(b)) as the test images. They can be viewed as the representative images with less and more textures, respectively. We apply the proposed model with α ranging from 2 to 2.3 and β ranging from 0.7 to 1 to those two images, and then observe the changes of PSNR and SSIM. From Fig. 8, we see that for Lena image, both PSNR and SSIM achieve the largest values when α = 2.1 and β = 1; and for Barbara image, both PSNR and SSIM achieve the largest values when α = 2.3 and β = 0.7. The results indicate that for different images, the fractional orders of the regularization term should be chosen differently. The fractional orders α and β could be close to the integers 2 and 1, respectively, for the images with less textures, and α larger than 2 and β less than 1 for the images with more textures. We may set α = 2 + and β = 1 − with a small for the latter case. 7. Conclusion. In this paper, we proposed a variational model, which uses a combination of two different fractional-order derivatives in the regularization and edge detector to perform different type of smoothing at different regions (homogeneous, near edges, textures). The idea of the proposed framework is based on the analysis of the filter character of the fractional-order derivatives. The experimental results with the comparison of several related models showed that the proposed model has better performance in preserving edges and textures, and avoiding the staircase effect when removing the noise, either from the visual effect or from the quantitative analysis. But in this paper, the selection of the fractional orders is not automatic. In our further work, we will try to discuss the adaptive parameter selection methods for the fractional orders in the proposed variational model.