Nonlinear fractional diffusion model for deblurring images with textures

It is a long-standing problem to preserve fine scale features such as texture in the process of deblurring. In order to deal with this challenging but imperative issue, we establish a framework of nonlinear fractional diffusion equations, which performs well in deblurring images with textures. In the new model, a fractional gradient is used for regularization of the diffusion process to preserve texture features and a source term with blurring kernel is used for deblurring. This source term ensures that the model can handle various blurring kernels. The relation between the regularization parameter and the deblurring performance is investigated theoretically, which ensures a satisfactory recovery when the blur type is known. Moreover, we derive a digital fractional diffusion filter that lives on images. Experimental results and comparisons show the effectiveness of the proposed model for texture-preserving deblurring.


1.
Introduction. Image deblurring is a classic image processing problem that has been researched for a long time and yet no clear cut solution exists due to its illposedness. Generally, the image blurring process is modeled as the convolution of a clear image with a shift-invariant kernel plus noise, i.e., where K is a smoothing operator, is noise, and u is to be recovered. The operator K is an integral operator of the form where k(x − x ) is a translation invariant blurring kernel, which is also known as a point spread function (PSF). The image restoration problem is to obtain a reasonable approximation of the original data u(x) from given recorded blurred, noisy data f (x) and known PSF k(x). We note that in some cases the PSF k(x) will also need to be reconstructed, the so-called blind deconvolution problem, but this is beyond the scope of the present analysis. In some applications (such as in artistic painting restoration [36]), it can be considered that there is not noise involved in the process, and the problem is then a pure deconvolution one: (1) f = Ku.
While much progress has been made recently, restoring textured images from blurred ones is still a very challenging problem where available methods tend to output over-smoothed images. In this work, we devote ourselves to the reconstruction of a textured image from the blurred data. Meanwhile, we do not restrict ourselves to a specific kind of blur but view this problem from a more general point of view in order to cover common principles in sharpening various blurs.
To stabilize the restoration process, some a priori knowledge about the unknown image is utilized through the addition of a regularization term, resulting in the model min where Φ reg (u) is the regularizer that enforces the a priori knowledge and the parameter µ is used to balance the two terms. Two classes of regularizers are well known. One is the Tikhonov-like class [44] including Φ reg (u) = Σ j D (j) u 2 , where D (j) 's stand for some finite difference operators and the order j is taken from 1 to the total number of operators. Another well-known class of regularizers is based on total variation (TV), which was first proposed for image denoising by Rudin, Osher and Fatemi in [39], and then extended to image deblurring in [40]. In comparison to Tikhonov-like regularizers, TV regularizers can better preserve sharp edges or object boundaries that are usually the most important features to recover. However, due to the nondifferentiability and nonlinearity of the TV functions, the TV based model is computationally more difficult to implement. To overcome this difficulty, a number of numerical methods have been proposed, such as first order primal dual method [6] and fast total variation minimization (FastTVMM) method [24].
Besides the variational approach, many alternative concepts were developed as well, for instance filtering by smoothing operators (one of the simplest cases is Gaussian convolution) or approximation by appropriate basis functions (the wide field of wavelet approximation falls into this category). For instance, Cai et al. [5] proposed a new wavelet frame based model (WFBM) that explicitly treats images as piecewise smooth functions, which estimate both the image to be restored and its singularity set. However, WFBM will also introduce spurious noise spikes and oscillations (very like pseudo-Gibbs) artifacts in the vicinity of discontinuities of images [10]. As we know, conventional local regularization strategies including TV-based and wavelet-based methods emphasize spatial adaptation from a local smoothness perspective. The limitation of such local approaches is that important structures in images such as edges and textures manifest nonlocal dependency or regularity (e.g., similar patches along the edge orientation or within the texture region) [27].
To the best of our knowledge, there are only a handful of previous work that consider texture preservation in the context of image deblurring. Examples include the work by Honigman et al. [23], where a potential function is incorporated into the diffusion equation which effectively prevents smoothing of specific texture details in an image. Another example by Elad et al. [18] involves handling textures from a separate dictionary by sparseness-based approaches. To date, nonlocal-based methods have been the most promising deblurring approaches for recovering texture.
Two most popular nonlocal regularization techniques in the literature are nonlocal total variation (NLTV) [19,50] and block-matching 3-D (BM3D) [11,12]. In order to restore image details effectively, the NLTV model makes use of the similar image structures by combing the variational framework with the nonlocal self-similarity constraint [19,9,38,29].
Although both NLTV and BM3D have shown excellent performance in image recovery applications, the minimization of the NLTV functional remains as a difficult optimization problem due to the computation complexity and the nondifferentiability [50]. In [13], BM3D frames are defined explicitly and based on a generalized Nash equilibrium approach, the two objective functions for denoising and deblurring parts are balanced. This algorithm is one of the best existing deblurring methods for symmetric blurs (e.g., Gaussian and out-of-focus blurs). Nevertheless, this approach has one major flaw: the prior is imposed on intermediate (patch) results, rather than on the final outcome, and this typically results in visual artifacts. To address this problem, Papyan et al. [31] proposed a multi scale patch based (MSPB) model using learning models of natural image patches. By using a local model on various scales, they managed to keep the simplicity of a low dimensional model while enjoying the non-locality of bigger regions in the image.
As we know, TV based methods rely on the image gradient, which usually models the interactions between pixel pairs. However, it is difficult to model more complex structures of natural images only using adjacent image pixels. As a compromise between the first-order TV regularized models and the high-order derivative based models, fractional partial differential equations based methods were widely used in the fields of image processing. For instance, Bai et al. [3] first introduced fractionalorder derivative into anisotropic diffusion equations for noise removal, where c(·) denotes the diffusion coefficient and D α * x denotes the adjoint operator of D α x , which may be viewed as a generalization of Perona and Malik (PM) model [32]. In addition, Dong et al. [17] proposed 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 variational model was described as follows: where the orders of the fractional derivatives α > 1, β ∈ (0, 1], µ > 0, g(x, y) is an edge detector and λ > 0 are parameters balancing the contribution of the three terms in the objective function. That work suggested that fractional calculus could nonlinearly enhance the complex texture details of fractal like structures. In recent, there have been several works involving discrete forms of an α-order derivative proposed to tackle the image denoising problem [21,22,25,7,8,51]. However, much fewer work employing partial fractional α-order derivatives was applied to the image deblurring problem. In [46], Williams et al. proposed a new method for image deconvolution using fractional derivatives to build a regularisation technique, and then they improved the model by enforcing intensity-range constraints implicitly in the problem functional. Roero et al. [37] presented a method for deconvolution of images by means of a regularizing term with a fractional power of the Laplacian, which ensured the model could recover higher frequencies. Moreover, Yin et al. [49] established a novel class of fractional-order nonlinear anisotropic diffusion equations based image restoration model, which employed the p-Laplace norm of fractional-order gradient of an image intensity function. In that work, the role of the fractional-order gradient was to better accommodate the texture details of an image, and the adaptive factor p was used to diffuse adaptively according to local geometry features, which were fractional-order curvature and fractional-order gradient of an image.
These works demonstrate good performance of the fractional-order derivative in achieving a satisfactory compromise between eliminating staircase effects and preserving important fine-scale features such as edges and textures. Based on the above feature and from the view of system evolution, we implement a fractional diffusion equation based deblurring model.

1.1.
Motivation and contributions of our approach. Lately, Zhou et al. [53] established a diffusion framework for steering the noise removal progress by the diffusion equation rather than the variational method. Specifically, they proposed a doubly degenerate diffusion equation model which takes into account gray level information by introducing a gray level indicator. Their work indicates that nonlinear diffusion equations are effective tools for image denoising. However, to the best of our knowledge, there are few work based on nonlinear diffusion equations to solve image deblurring problems. An example is the so-called anisotropic filters such as forward-backward-forward diffusion in [20], which can be regarded as an interpolation of the PM model and TV model. This model offered two regularization parameters which makes it possible to harvest the desirable sharpening features of PM models in a controlled fashion. Inspired by these work, we consider the image deblurring problem in the view of anomalous diffusion behavior. In particular, we develop a novel deblurring framework based on a class of fractional diffusion equations. The capability of preserving the edges and texture details of the fractional partial differential equations is an important advantages that leads to the superiority of the proposed approach compared to the traditional integer-order computation-based deblurring algorithms, especially for images with rich texture details.
Our work aims to provide more insight into the effect of texture preservation that fractional diffusion based methods have on the image deblurring process. In addition, we seek to provide a theoretical analysis of the effectiveness of the method by estimating the error between the original image and the resulting one. In Section 4, experimental results on five benchmarks demonstrate that the proposed algorithm based on the fractional diffusion equation performs favorably against the stateof-the-art deblurring methods. We expect the interplay of image processing and fractional diffusion equation is going to be more fruitful in the future.
The contributions of this work are summarized as follows: 1 A novel framework based on fractional nonlinear diffusion equations is developed for the textured image deblurring problem. 2 The approximation properties between the original image and restored image are analyzed in both spatial and frequency domain theoretically, which provides a strong support for the achievement of a satisfactory recovery.  3 An efficient algorithm based on the finite difference method is developed by virtue of the discretization form of fractional derivative operator. 4 The method can be easily adapted to handle various blurs, since the diffusion equation involves a source term containing the information on blur. Based on these four contributions, the proposed deblurring scheme is able to achieve high-quality deblurring results of textured images by taking into account more information of images. As shown in Figure 1, our scheme successfully recovers sharp edges, preserves fine details, and while prevents ringing artifacts.
The rest of this paper is organized as follows. In Section 2, we propose a novel nonlinear fractional diffusion equation based method for solving the image deblurring problem. We shall give detailed arguments for the effect of the regularization parameter on the deblurring capacity. In Section 3, we introduce the corresponding discrete fractional diffusion filter on images. In Section 4, we show the performance of the new filter in comparison with other deblurring methods. Finally, we draw a conclusion in Section 5.
2. Nonlinear fractional diffusion model. We develop the following framework for image deblurring: on Ω, is the bounded domain of the image and has a Lipschitz boundary. Here K is a convolution operator with nonnegative kernel k ∈ L 1 (Ω) satisfying Ω k(x)dx = 1. In addition, we notice the term u r M r is regarded as a gray level indicator. The information about operators div α and ∇ α is present in Appendix A.

Why use the fractional derivative as a regularization of diffusion?
The motivation to use the fractional derivative is that the model can preserve and even enhance the texture details of original images. Generally, deblurring algorithms attempt to solve the following inverse problem [45]: for a degraded image, f ∈ U, where U ⊂ R N ×N , u ∈ U is a possible solution and k ∈ R n×n is a blurring kernel. The parameter λ is a Lagrange coefficient and g(u) is a regularization function. Partial differential equations (PDE)-based approaches generally use g(u) = G(|∇u| 2 ). A gradient descent minimization then yields the following scheme [45]: where k (η 1 , η 2 ) = k(−η 1 , −η 2 ) and G (|∇u| 2 ) is the diffusivity or edge detection function. This function is chosen so that high gradients are preserved and low gradients are smoothed.
As we know, most of the existing models work well in deblurring for nontextured images, but fail in textured image restoration. Unlike the integer-order derivative operator, the fractional-order derivative operator possesses the "non-local" property because the fractional-order derivative at a point depends upon the characteristics of the entire function and not just the values in the vicinity of the point [33]. This "non-local" property is beneficial to improve the performance of texture preservation. Some progress in the studies on fractional-order image processing not only validates the fractional-order partial differential equations, but also provides interesting and practical suggestions for future research. For instance, a fractional differential provides the flexibility of enhancing the complex texture details of an image in a nonlinear manner [35,48,34], it can maintain the low-frequency contour features in the smooth area of an image in a nonlinear fashion, and creates the possibility of enhancing the high-frequency edges and texture details, in a nonlinear manner, in those areas where the grey level undergoes frequent or unusual variations [35,48,34]. Therefore, to solve the aforementioned common problems of the traditional integer-order contrast deblurring algorithms, it is natural to consider whether it will be possible to hybridize the capabilities of fractional calculus to preserve the edges and texture details with texture image deblurring.
For comparison, we apply two models with fractional and integer gradients in the experiment. The aim of this example is to present the advantage of the nonlocal property of fractional derivative for preserving the texture details of an image. Therefore, we take a signal with rich textures, i.e. high frequency components, in the experiment. The two diffusion based models are described as follows.
The fractional model: The integer model: In the run, we take α = 0.8, r = 0.7, k 1 = 0.01, β = 0.9, ∆t = 0.02, T = 3. As shown in Figure 2, most of high frequency oscillations are preserved by applying the fractional model. Additionally, numerical experiments suggest that the fractional model produces significantly less staircasing than the integer model does. Obviously, the model based on fractional diffusion equations is able to avoid over smoothing in processing texture images. This is most likely due to the fact that the corresponding kernels in fractional derivatives, while still singular, are non-local and thus provide some degree of averaging. In addition, by choosing a proper λ, we can balance the texture enhancement and deblurring abilities of the model to get a better reconstruction. For more details about how to chose parameter λ, we refer to the next subsection.

2.2.
Influence of the regularization parameter λ. As indicated by many reseachers, the difficulty with the restoration method of considering the TV regularization is the determination of the important but unknown regularization parameter, which controls the trade off between a good fit and an oscillatory solution. The work of [24] illustrates that real images cannot be efficiently represented by homogeneous models, since they are made of various textures, sharp edges and constant areas. So the prior parameter should adapt to the local structure of the image in order to provide a better reconstruction (i.e. a less noisy result in homogeneous regions, and sharper details in other areas). However, an adaptive regularizing model based on statistical control parameters [28] seems to be too computationally demanding, both for estimation and reconstruction. We perceive that most existing work empirically provides a range of the regularization parameter without any theoretical guidance, even though the results are acceptable in some situations. Therefore, we shall present a theoretical analysis of the new digital filter and prove convergence for suitable regularization parameter λ balancing the diffusion and the source term. The theoretical analysis is performed both in physical and frequency domain.
In order to discuss the problem in spatial domain, we review a function space Q s p (Ω) (refer to [52]). For any positive integer p ∈ N + , let For more details, we refer to [52]. We assume that the exact solution u 0 (x, y) of (1) satisfies u 0 (x, y) ∈ Q α 2 (Ω). To the best of our knowledge, there is lack of theoretical support for analyzing the well-posedness of the proposed equation. Some existing analyses of the doubly degenerate diffusions with integer order, i.e. α = 1, can be found in [4,14,15,16]. Ideally, we assume the steady state of the proposed equation can be attained after a long time of diffusion. In other word, we assume that there is a function v(x, y) ∈ Q α 2 (Ω) satisfying the steady state which leads to the estimation In this sense, if the blurring operator K is injective, i.e., there is a positive constant In conclusion, the solution of (2) approaches that of (1) when λ is close to infinity.
On the other hand, in order to analyze the influence of λ for the model in frequency domain, we replace derivative operators in classical image processing approaches with operators of the type D α (please refer to Appendix B for more details) and investigate the properties of the resulting filter method. As shown in [3], assume the domain Ω has certain regularities, then the extension theorem [1] ensures there exists a linear operator E which maps H 2α (Ω) onto H 2α (R 2 ). So we can prolongate u ∈ H 2α (Ω) to u ∈ H 2α (R 2 ) by the linear extension operator E. Therefore, we consider the equation in the whole space R 2 . To simplify, we analyze the model by taking r = 0 and k 1 = 0. Then the equation in the model becomes the following parabolic equation involving the fractional power of the Laplacian: Next, we demonstrate the influence of λ for the proposed model. More specifically, we will point out the effect of λ on the asymptotic behavior of the solution u(x, t), i.e. how the solution converges to the wanted result u 0 . Furthermore, the effect of the parameter λ on the restoration will be deduced, which provides a selection criterion of λ in practice.
Particularly, if we choose h = 1 and employ the central fractional derivative (9), then we can obtain the finite fractional-order difference as follows: Before we introduce the finite difference discretization of the fractional order derivative, we define a spatial partition (x l1 , y l2 ) (for all l 1 , l 2 = 0, 1, ..., N − 1) of image domain Ω. We consider the discretization of the α-order fractional derivative at all points of Ω along the x-direction and the y-direction by using the approach where ω α i = (−1) i C α i , i = 0, 1, ..., N − 1, and ω α j = (−1) j C α j , j = 0, 1, ..., N − 1, and Proposition 1. (see [30]) For any 0 < α < 1, the coefficients {ω α l } ∞ l=0 have the following properties: In practice, all N equations of fractional derivatives and their adjoint operators along the x direction in (4) can be written simultaneously in the matrix form Let U ∈ R N ×N denote the solution matrix at all nodes (l 1 h; l 2 h), l 1 , l 2 = 0, ..., N −1, corresponding to x-direction and y-direction spatial discretization nodes. Therefore, we have ∆ αx U = B α N U . Similarly, the αth order y-direction derivative of u(x; y), In addition, the adjoint operators can be rewritten as the matric form: ∆ * αx U = (B α N ) T U and ∆ * αy U = U B α N . Remark 2. As shown in [52], the standard definitions for fractional derivatives require a function to have zero Dirichlet boundary conditions due to end singularity, but for imaging applications such conditions are unrealistic and too restrictive. Since the fractional derivative is nonlocal, any inappropriate boundary condition will influence the whole image. Therefore, we calculate the fractional-order difference for all points in the image based on (4) and (5). Experimental results show that such a treatment could enable the model to get satisfactory results to some extent. Future work will focus on finding an effective treatment of the boundary conditions for such problems involving fractional derivatives.

3:
for i, j = 1, ..., N − 2 do 4: if u n1+1 doesn't satisfy the given stopping condition then end if 12: end for 4. Experimental results on textured image deblurring. We present experimental evaluations of the proposed algorithm against the state-of-the-art deblurring methods for textured images. In experiments, we consider five deblurring scenarios used as the benchmarks in many publications. The PSF for each scenario are summarized in Table 1. The results are compared with those obtained by some state-of-art algorithms from past research such as the FastTVMM method [24], the MSPB method [31], the WFBM method [5] and the BM3D method [11]. The source codes of these methods are either publicly available on websites or requested from the authors.  Figure  3. In the proposed algorithm, we set ∆t = 0.0005, α = 0.9, k 1 = 0.0005, λ = 4000, r = 0.9, β = 0.8. The iteration termination condition of the proposed algorithm is where mean(|u n1 − u n1−1 |) is the average of the absolute values of u n1 − u n1−1 . In all experiments, the choice of tol only depends on the noise level σ. Specifically, we choose tol = 0.001, 0.06, 0.35 in the recovery of noisy blurred images with noise levels σ = 0, 2, 5, respectively.
For the FastTVMM algorithm, the MSPB algorithm, the WFBM algorithm and the BM3D algorithm, the free parameters are set as suggested in the reference papers [24], [31], [5] and [11]. The restoration quality can be assessed by visually comparing the restored image with the original image (which is provided here as a reference only, and is never used in the restoration procedure). PSNR (Peak Signal to Noise Ratio, unit: dB) is used to evaluate the objective image quality. We recommend the reader to carefully examine the different parts of the images (constant areas, sharp edges, small details and oriented textures) to observe how these parts are restored.
Firstly, two sets of experiments are conducted to verify the performance of the proposed method for image deblurring. To be more precise, the blurred image k * u is generated in MATLAB by "imfilter(u, k, 'circular')". In Figures 4-7, we show the results of our method for four test images and for five blurring kernels. We see that our method works well for a wide range of images and blurring kernels. The reason is that the diffusion term in the model enables our scheme to be fully independent from any assumption on the spatial distribution of kernels.
Next, we compare the proposed model with a closely related model-fractional variational model. In fact, in the case of r = 0, β = 1, the proposed model is equivalent to the evolution equation solving the total fractional-order variation model: which was studied in [52]. In order to check the equivalence, we apply both two models to restore Boat image degraded by scenario I. In this experiment, we set α = 0.9 and λ = 10 for both models. As shown in Figure 8, the recovery of the proposed model is almost the same as that of the fractional variational model [52] in term of PSNR.  Table 2. Values of PSNR for the restoration by using the proposed model with α = 0.9, k 1 = 1, β = 1, λ = 5.
Images r = 0 r = 0.5 r = 1 r = 1.5 r = 2 In addition, since we take r > 0 in the model, the proposed model is totally different from the variational model. This is due to the fact that the model can not be regarded as the evolution equation of any fractional variational model. In this sense, we remark that the model is established in the framework of fractional diffusion equations rather than fractional variational models. For comparison, we applied our model with different parameter r to restore all six test images degraded by scenario I. Here, we regard the result of our model with r = 0 as that of fractional variational model, which is reasonable due to the equivalence we mentioned above. Results in Table 2 illustrate that the proposed model (r = 0) performs better than the fractional variational model (r = 0) in most cases.
Next, we compare our method with the closest competitors FastTVMM [24] and WFBM [5]. In the experiment, we restored Barbara image and Aerial image which are degraded by application of scenarios II and III. The visual quality of some of the restored images can be evaluated from Figures 9-12. As we know, FastTVMM [24] is  a TV-based deblurring approach that can reconstruct the piecewise smooth regions well but often fails to recover fine image details. WFBM [5] is a new wavelet frame based image restoration model that explicitly treats images as piecewise smooth functions. For Barbara image (see Figures 9-10), our method noticeably outperform the other two. Between our method and others, the difference mainly comes from the more faithful reconstruction of some low-contrast texture areas such as the scarf areas with shading and in the shadow. For the Aerial image (see Figures 11-12),  the subjective quality of deblurred image by our method is apparently superior to that of other competing schemes, ours does not contain any artifact in the lower left region like others do. For a detailed comparison, we give in Figure 13 the zoomed versions of two parts of Barbara and one part of Aerial. Comparing the texture details of the deblurred result and the test image, we can clearly find that our method restores the image successfully while preserving significantly more texture details (see .   (a) (b) (c) (d) Figure 9. Deblurring of the Barbara image, scenario II. From left to right: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 10. Deblurring of the Barbara image, scenario III. From left to right: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 11. Deblurring of the Aerial image, scenario II. From left to right: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 12. Deblurring of the Aerial image, scenario III. From left to right: blurred, reconstructed by FastTVMM [24], WFBM [5], our method.
(a) (b) Figure 13. The original Barbara image and Aerial image, the marked rectangles are marked for zooming.
Then we consider the noisy blurred case, where blur and noise are assumed to be known. The noise is considered to be Gaussian white noise with variance σ 2 . In Figure 21, two types of blur were used: scenario V and scenario II with noise level set at σ = 2 and σ = 5, respectively. In Figure 22, scenario I was used with the noise level set at σ = 2, 5, respectively. The restoration in Figures 21-22 shows that the proposed model can preserve textures as well as remove the noise. Peak signal-to-noise ratio (PSNR) values for the experiments shown in Figures 21-22 are Figure 14. Zooming regions. Figure 15. Deblurring of the Barbara image, scenario II. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 16. Deblurring of the Barbara image, scenario II. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method.
given in Table 3. Results presented in Table 3 illustrate that the proposed algorithm performs well in terms of PSNR. In Figure 23 another comparison is offered that emphasizes the efficiency of the proposed model as a recovery tool for noisy blurred images. Table 4 offers a side-byside comparison of the corresponding PSNR values. The proposed model delivers better results in a shorter amount of time compared to the reference algorithms. This proves that a PDE-based method can indeed be very efficient. Figure 17. Deblurring of the Barbara image, scenario III. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 18. Deblurring of the Barbara image, scenario III. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 19. Deblurring of the Aerial image, scenario II. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 24 shows that the proposed method not only recovers a better image, but also computes its output (on the right) in only 2 seconds. This is to be compared with the 220 seconds it took MSPB to deliver its output image.

5.
Conclusion. In this paper we try to solve the textured image deblurring by means of a fractional diffusion equation. Our contributions are three-fold. Firstly, we develop the diffusion based method by introducing the fractional derivative and Figure 20. Deblurring of the Aerial image, scenario III. From left to right, zoomed fragments of the following images are presented: blurred, reconstructed by FastTVMM [24], WFBM [5], our method. Figure 21. Top row: Noisy blurred images. Bottom row: Recovered images corresponding to the above images. Left two columns: scenario V and noise levels σ = 2, 5, respectively. Right two columns: scenario II and noise levels σ = 2, 5, respectively.  Figure 22. Top row: Noisy blurred images with scenario I and noise levels σ = 2 (left two columns), σ = 5 (right two columns). Bottom row: Recovered images corresponding to the above images. Figure 23. Left to right: Noisy blurred image (scenario I and σ = 2), recovered image from FastTVMM [24], MSPB [31], WFBM [5], BM3D [11], and recovered image with proposed model. the gray level indicator, which performs well in both deblurring and texture preservation. Secondly, we propose to understand the restoration capacity of the method by analyzing the approximation properties between the result solution and the original one theoretically. From the theoretical perspective, our deblurring technique is essentially a new diffusion process with a proper source term, which makes it possible to obtain a satisfactory result when the diffusion time and the parameter  There are several open issues that remain to be addressed in the future. Among them, the choice of relaxation parameter α has been shown critical to the performance but the rule of choosing α has seldom been studied. There are also several possible extensions along this line of research, e.g., in this work we have exclusively considered non-blind scenarios; potential implications into blind image deconvolution [26] remains unexplored. We note that our approach could be extended for recovering images degraded by double (even multiple) blurs. In those cases, the source term needs to be modified to cater for processing such problems.
Finally, as the proposed method essentially solves the pure deblurring problem, we wonder whether the output of the algorithm could be injected into a denoising algorithm to process the blurred and noisy problem. In this case, it could be possible to build an iterative algorithm. The stability of such a method is not certain and convergence studies have to be carried out.

Appendices
Appendix A. In this appendix, we give information on some fractional derivatives, see [52] for more details.
Let Γ(·) denote the Gamma function. For any positive integer n, and a real number s between n − 1 and n, i.e. n − 1 < s < n, the GL derivative, the Caputo derivative and the Riemann-Liouville (RL) derivative of order s are respectively defined as follows.
The right-sided GL derivative: The left-sided Caputo derivative: The right-sided Caputo derivative: The right-sided RL derivative: The Riesz-RL (central) fractional derivative Remark 3. One notes that when a function is n−1 times continuously differentiable and its nth derivative is integrable, the fractional order derivatives by the above definitions are equivalent subject to d i v(a) dx i = d i v(b) dx i = 0, i = 0, 1, ..., n − 1 [33].  [52] implies the operator (−1) nC div s is the adjoint of operator R ∇ s . The superscript C is used for Caputo derivative based quantities such as C div s and C ∇ s , while superscript R means that a quantity is based on the RL derivative.
Based on Remark 4 and two kinds of Riesz fractional derivative in (8) and (9), we apply the Riesz-Caputo derivative for the operator div s and the Riesz-RL derivative for the operator ∇ s throughout this paper.
Appendix B. In this appendix, we introduce some basic notions related to fractional powers of the Laplacian. We consider the Fourier transform of a function f ∈ L 1 (R d ) pointwise defined by It is clear that We note that the Schwartz space of rapidly decreasing functions is defined as follows The set of all continuous linear functionals on S(R d ) is denoted S (R d ) and is called the space of tempered distributions. It is easy to verify that F : S(R d ) → S(R d ).
We define F −1 by We now give the relation between the Fourier transform and convolutions. The convolution u * v of two functions on R d is defined by (10) u Computing the Fourier transform of (10) leads immediately to the formula F(u * v)(ξ) = (2π) d/2û (ξ)v(ξ).
Following the notation in [43], we define the Sobolev space From [47] we deduce the spectral decomposition of the Laplacian as multiplication operators in the Fourier domain. For a detailed overview, see [41,42].