Discrete Laplacian Operator and Its Applications in Signal Processing

,


I. INTRODUCTION
Two main operations and building blocks in many engineering disciplines in general [1], [2], image processing and computer vision tasks in particular [3], are first and second order derivatives, otherwise known as the Gradient and Laplacian. The ubiquity of these operators come from the way we model our systems and problems. Mathematical models relying on the language of calculus are at the center. Some of these modeling paradigms include variational methods, partial differential equations (PDE), statistical and linear/non-linear optimization models.
Over the years, many attempts have been dedicated towards the generalization of those operators to different settings. Among the generalizations is the graph Laplacian [4] which found numerous applications in signal and image processing [5].
Of relevance to this work, is the generalization of integer order differential and integral operators to fractional orders. Fractional Calculus (FC) is a 300 years old concept dating back to the days of l'Hôpital and Leibniz [6]. FC has received increased interests over the last 30 years mainly due to their The associate editor coordinating the review of this manuscript and approving it for publication was Peng Wang. long memory property [2]. For a more recent historical survey, the reader is referred to [7].

A. FRACTIONAL DERIVATIVES
The integer differential and integral operators are defined uniquely, and they are local. In other words, they consider the values of very close neighboring points to the point of interest. On the other hand, fractional-order differential operators are non-local; larger neighborhoods are considered in the computation resulting in long-term memory effect. This is one of the main reasons behind their appeal. There is a multitude of definitions of the fractional order derivatives [13], [14]. We are going to adopt the notation D ν to denote the derivative of a fractional order ν. The following is a list of the most common definitions: where ν k is the binomial coefficient and h = |h|e jθ is a complex number with θ ∈ (−π, π]. • Riemann-Liouville (RL) fractional derivatives (2) where n − 1 < ν < n and (.) is Euler's Gamma function.
• Fourier domain fractional derivatives where ω is the Fourier variable, F and F −1 are the Fourier and the inverse Fourier transform respectively. Most of the fractional derivative formulations start from a continuous formula which requires discretization to facilitate implementation.

B. FRACTIONAL LAPLACIAN (FL)
The idea of extending the standard Laplacian operator to fractional order is an old idea (c.f [15] and references therein). The standard Laplacian is a local operator. Local operators utilize the immediate neighborhood only in the calculation where outliers can have big impact on the result. Many attempts have been made to create non-local operators that can alleviate the shortcomings of local operators in image processing applications [16]. Utilizing fractional-order Laplacian is such an attempt.
In a similar vein to the fractional derivatives, a multitude of fractional Laplacian (FL) definitions have been proposed over the past few decades but, a consensus on the most appropriate definition for an application is yet to be reached [17].
Of special interest to our work are two variants: the spectral fractional Laplacian [17] which is defined as follows: where is the continuous Laplacian operator applied on the function f (x). φ i and λ i are the eigenfunctions and eigenvalues of the continuous Laplacian respectively. Secondly, the Fourier-transform based definition (pseudodifferential) [17] It is important to note that both definitions are continuous, and a discretization step is required to make them applicable to digital data. For a comprehensive list of the different formulations, the reader is referred to [12], [15], [17] and references therein.

C. RELATED WORKS 1) FRACTIONAL DERIVATIVE OPERATORS IN IMAGE PROCESSING
Fractional-order derivatives have found numerous applications in image processing. These applications can be categorized into three main categories based on the framework under which the derivative is used. The first category is linear filtering. In this category, a fractional-order derivative of an image is calculated through a linear filtering process. In some applications, the gradient image is the goal such as in edge detection [18]- [31] which is the earliest image processing application of fractional calculus. Another application is to use the gradient image to enhance the input image through, for example an un-sharp masking scheme [26], [32]- [37]. The third application in this category is contrast enhancement [29], [38]- [40]. The second category is referred to as PDE-based models. Modeling of image restoration problems using integer derivatives dates back to the late 80s and early 90s but utilizing fractional-order derivatives was first presented in 2007 [41] for the image denoising application and was later refined, improved and adopted for a number of applications. Among these applications; image denoising [42]- [49], contrast enhancement [50], [51], image deblurring [52] and image super-resolution [53], [54].
The third category is related to the second and we will refer to it as variational models for image restoration problems. Similar to PDE(s), variational models with integer derivatives in the context of image processing date back to the 80s, but the introduction of fractional-order variational models is recent. Applications in this category include image denoising [55]- [66], in-painting [58], [67], fusion [56], [64], nonrigid registration [68], super-resolution [56] and optical flow estimation [69].
For the sake of completeness, we cover a fourth category that is based on a global optimization technique namely: fractional-order Darwinian particle swarm optimization (FODPSO). The reason for distinguishing this category is that the fractional-order derivative is not applied to the image but to a parameter of the model at hand such as the optimal threshold. The main application in this category is image segmentation [70]- [72].

2) FRACTIONAL LAPLACIAN OPERATORS IN IMAGE PROCESSING
The fractional Laplacian has not seen the same amount of adoption as is the case with the fractional derivative in signal and image processing. In [73], the authors introduced a scale-space model. In [74], the authors proposed a quadratic optimization model for blind image deconvolution involving the fractional Laplacian. In [75], the authors proposed a PDE model for vector field estimation flow-sensitive MRI imaging. In [76], the authors have demonstrated image denoising by solving a fractional diffusion equation (PDE). In [77], the authors proposed a variational model for image denoising based on the fractional Laplacian and later was extended to tomographic reconstruction involving the fractional Laplacian as a regularizer [78].

D. CONTRIBUTIONS
The contributions in this work can be summarized in the following: I) Motivated by the construction of the spectral Laplacian (4), we propose a discrete fractional Laplacian in the form of a matrix operator using the DCT transform. The discrete construction avoids the need for discretization which is typically required by the other fractional Laplacian construction techniques. The discretization step in the case of the pseudo-differential formulation in (5) usually involves solving a filter design problem [79]- [81]. II) Motivated by the reliance of the trend filter on the Laplacian operator, we develop a computationally efficient implementation of traditional trend filtering in the DCT transform domain. III) Utilizing the proposed discrete fractional Laplacian, we extend both the traditional and the 1 trend filters to fractional order and demonstrate their effectiveness at the image denoising task for higher levels of noise. IV) We finally demonstrate the applications of the proposed fractional Laplacian on a number of image processing tasks. However, it is important to stress that we are not claiming state-of-the-art performance in any of these applications. Our aim is to demonstrate the potential of the proposed operator. The remainder of this article is organized as follows. Section II introduces the fractional Laplacian in the 1D and 2D cases. In section III, we provide two groups of applications for the proposed operator. In the first group namely: fast trend filters, we first introduce the trend filtering followed by a DCT-based implementation of the 2 trend filter. We then generalize the trend filtering by proposing a fractional 1 and a fractional 2 trend filters. In the second group, we present a number of image processing applications demonstrating the effectiveness of the proposed fractional Laplacian. Lastly, a discussion and conclusions are presented in section IV.

II. A DISCRETE FRACTIONAL LAPLACIAN OPERATOR A. DEFINITION FOR 1D CASE
The Laplacian operator of single variable continuous function f (x) is defined as follows: To implement the Laplacian operator, an approximate discretization is used as follows: where L is the discrete Laplacian operator and h is a small constant. This discretization amounts to uniformly sampling the function f (x) with distance h. This operation can be implemented as a convolution as follows: where l = [1, −2, 1] is an FIR filter. The Laplacian in (8) can be formulated as a matrix operator in few different ways based on the boundary condition assumed. In this work we choose a specific formulation with desirable properties as we are going to see later. Specifically, we formulate the Laplacian operator L ∈ R N ×N as follows: where, the rows between [1, N −2] 1 are shifts of the FIR filter discussed earlier. The signals we are dealing with are finite, and as we have seen in (7), the Laplacian filter is centered around x. In other words, to perform the convolution in (8), f (x) needs to exist before and after x which, is not true at the boundaries. As a solution, we extend the signal on both sides [82]. Thus, the first and last rows in (9) are different from the middle rows. This signal extension near the boundary can be performed in several ways [83]. In (9), we assume a symmetric signal extension as we are going to see next.
The operator in (9) has the same first and last rows, this means equal boundary condition on both sides. For a signal u ∈ R N ×1 , the boundary condition in (9) assumes the signal u has zero slope at the boundary (u(0) = 0,u(N − 1) = 0) and is symmetrically extended around the midpoint (u(−1) = u(0), u(N ) = u(N − 1)). For more details about the boundary conditions, the reader is referred to [83].
It is important to note that this particular operator in (9) is diagonalizable by the DCT type-II matrix [83] denoted by M L = M T EM (10) where E is a diagonal matrix of which the diagonal elements denoted {e(k)} are the eigenvalues of L. This spectral decomposition is real which results in efficient computation. Motivated by the spectral Laplacian in (4) and the diagonalization property of the discrete Laplacian in (10), we define the fractional Laplacian L α 2 as follows where E α is diagonal matrix with the diagonal elements raised to the power α.
To consolidate, we started from the discrete Laplacian filter l = [1, −2, 1] which we used to formulate the matrix L in (9). Then we performed the eigen-decomposition in (10) and we raised the diagonal eigenvalues matrix E to the power α to finally get the fractional Laplacian matrix L α in (11).
Recall that when α = 1 in (11) we recover the standard Laplacian in (9) but, for α in the open set α ∈ (0, 1) we get different fractional Laplacian filters. The fractional Laplacians we get are still diagonalizable by the DCT matrix as in (11) and share the same structure as the standard Laplacian 1 Numbering starts from 0 2 We are using the subscript notation for α to avoid the confusion with the exponent as L is a matrix operator in (9) with the difference that smaller values of α result in less sparse operator L α . In other words, smaller α results in Laplacians with more nonzero coefficients hence the memory effect.
As a result of the increase in the number of coefficients, two things become clear: firstly, the dimensions of the fractional Laplacian need to be larger for smaller α to reduce the numerical inaccuracy due to truncation. Secondly, the number of rows in L α that are involved in the boundary condition increases. One way to recover a good filter l α from the matrix L α is to choose the middle row.
It is important to note that the eigenvalues of L α ∈ R N ×N can be obtained by the relation e(k) = 2 − 2 cos k π N [84]. Thus, the eigenvalues of the fractional-order Laplacian L α are defined as follows: (12) and the coefficients of the eigenvectors matrix are the same as DCT-II matrix which are as follows: To get the FIR filter l α from the operator L α , we choose the N −1 2 th row of L α based on our earlier discussion about avoiding the boundary condition. It can be readily demonstrated that the coefficients of the l α have the following form:

B. EXTENSION TO THE 2D CASE
Generalization to the 2D case is straight-forward as the 2D-DCT is usually performed as 1D-DCT calculated once on the rows and once on the columns [83]. The benefits of this formulation of the fractional Laplacian in the DCT domain over the traditional Fourier transform formulation are twofold. First, we can construct the fractional Laplacian filter l α easily allowing for flexibility in applications. Second, it is computationally more efficient as filtering in the DCT domain avoids the need for the image extension as required by the DFT [85].

C. COMPUTATIONAL COMPLEXITY ANALYSIS
From the earlier development, we notice that the proposed fractional Laplacian has two equivalent forms. The first being a matrix operator L α ∈ R N ×N which is constructed using (11). To calculate the fractional Laplacian of a signal using L α , a matrix-vector multiplication is performed, which has O(N 2 ) computational complexity. The second form is a linear FIR filter l α ∈ R N which can be constructed directly using (14). Calculating the fractional Laplacian of a signal of length M using l α is a convolution operation with a computational complexity of O(NM ).

D. DISCUSSION
We have delayed the discussion about the other formulations of the matrix operator L for a reason that will become clear soon. The alternative to the proposed Laplacian is to formulate it as a circulant matrix, assuming that the signal has a periodic extension as follows [83]: The L matrix in this formulation is diagonalizable by the discrete Fourier matrix (DFT) [83] which could be generalized the same way we did in (11). However, the proposed Laplacian is computationally more efficient. This computational efficiency comes from two aspects. Firstly, the DFT matrix is complex while, the DCT is real. Secondly, the boundary condition in (15) assumes that the signal is periodic, in other words, the signal needs to be padded with a copy of itself on both sides creating discontinuities along the boundaries.
On the other hand, the Laplacian in (9) assumes that the signal is padded with a mirror-image of itself along the boundary. This is a natural boundary condition and it is what MATLAB uses in the imfilter command with the symmetric option.
In [86], the authors have proposed to compute the fractional derivatives of images implicitly by utilizing a discretization of the following result [87]: where ν is the fractional exponent and ω is the continuous frequency variable. In other words, in contrast to what we do in this work, the authors do not provide a direct method for constructing the fractional differential operator (D ν ) in [86]; rather a method to compute its effect on a function. Being able to construct the operators affords more flexibility and applicability in linear algebraic settings.
In [88], the authors propose to approximate the function (16) in the DCT domain by treating it as a FIR filter design problem. This technique produces coefficients of a FIR filter which approximates the ideal fractional derivative operator. The main difference between this technique and ours is that our technique does not use an approximation and it generalizes the Laplacian operator rather than the first order derivative.

III. APPLICATIONS
The Laplacian operator is very ubiquitous in applications and generalizations of which can be examined on such applications. We have split the discussion about applications into two groups. The first group deals with a specific trend estimation technique which relies on the standard Laplacian. We introduce it briefly, propose a new implementation, and VOLUME 8, 2020 generalize it. In the second group, we consider five image processing applications.

A. FAST TREND FILTERS IN THE DCT DOMAIN 1) DEFINITION
Estimating the trend of a signal or time-series is a common problem in many disciplines [89]. More specifically, given a signal b ∈ R N ×1 , it is assumed that the signal is composed of two components: a slowly varying component u known as the trend and a rapidly varying white Gaussian noise component denoted as follows: Trend estimation is the process of producing an estimateû for the underlying trend u. The literature on trend estimation is very rich with various parametric and non-parametric techniques [90]. In this article however, we are interested in a specific trend estimation technique which we will refer to as trend filter.
Concretely, we use capital letters to represent matrices and small letters to represent column vectors. For example, the nth column vector of the matrix M is denoted M n , while the nth element of a vector u is denoted u n . The p −norm of u is ||u|| p . The trend filter [89], [91] is formally defined as follows: A special case of the model (18) is the generalized Wiener filter [92], [93] when p = 2, which is otherwise known as the Hodrick-Prescott trend filter in the statistics community [91]. By utilizing the fact that the Laplacian matrix L is diagonalized by the DCT matrix (10), the cost function can equivalently be defined in the DCT domain: whereū = Mu andb = Mb. By minimizing the cost function, we haveū This result shows that the trend filter can be very efficiently implemented in the DCT domain as an element-wise division since the eigenvalues e(k) can be pre-calculated. This result also allows us to interpret the trend filter as a low-pass filter in the DCT domain. To this end, let G = M T , the trend filter can thus be represented as follows: On the other hand, the observed vector b can be represented as Comparing (23) with (25), we can clearly see an interpretation of the trend filter as a shrinkage operation in the DCT domain. The level of shrinkage is controlled by the regularization parameter λ. Since each column vector G(k) can be regarded as a frequency component with k = 0 corresponding to the DC and k = N − 1 corresponding to the highest frequency, the shrinkage is also frequency dependent. It can be shown that the eigenvalues have the property e(0) = 0 and e(i) < e(k) for i < k. As a result, a higher frequency component will be shrunk more. Similar to a linear low-pass filter which attenuates the high frequency components in the Fourier transform domain, the trend filter is a low-pass filter which attenuates the high frequency components in the DCT domain.

2) FRACTIONAL 2 TREND FILTER
As an application of the proposed operator in (11), we can generalize the trend filter (19) by replacing the L operator with L α operator leading to a new model. We call it the fractional trend filter which is the minimization of the following cost function: The solution to this cost function, can be performed in DCT domain as follows: Compared with (20), we have more control over the magnitude response by changing α. To demonstrate the impact of α on the shape of the filter, we plot the filter function (27) at various values of α in Fig. 1. To facilitate the comparison, we normalize the x-axis because, at each α we generate a new Laplacian L α that has its diagonal elements e(k) ∈ [0, 4 α ]. It is important to note that a similar model was proposed in [74] with the difference that the solution was done in the Fourier domain.
From Fig. 1 we notice that lower values of α allow more high frequencies to pass than the higher values of α.

3) FRACTIONAL 1 TREND FILTERING
We extend the 1 trend filtering [89] to fractional 1 trend filtering of the form: min u ||u − b|| 2 2 + λ|L α u| 1 To solve this problem, we adopt the ADMM algorithm [94] as follows: The shrinkage effect of the fractional 2 trend filter (first term in (27)) using various values of α.

Algorithm 1 ADMM for Fractional 1 Trend Filter Using Matrix Form
Result: u k+1

Algorithm 2 ADMM for Fractional 1 Trend Filter Using Convolutions
Result: u k+1 In algorithm 2, l α is one of the middle rows of L α . It is a FIR filter. In 1D, l α and its transpose l T α are convolved with signals. In 2D, we generate two kernels using l α for the x and y directions as shown Fig. 2. l T α is identical to l α as they are symmetric around the center. For images, we use an  an-isotropic extension of the 1D model as follows

4) DISCUSSION
We note that the fractional Laplacian in (14), being purely discrete, is a FIR filter of order N . To determine the order of the filter, we need to exploit the structure of the Laplacian. The Laplacian filter is even symmetric and has a positive and relatively high middle coefficient surrounded on both sides by negative coefficients which decay towards zero approaching both ends. Smaller values of α correspond to filters with sharper roll-off in the frequency domain and as a result the filters are generally longer in the discrete domain.
To determine the length of the filter l α , we propose an empirical procedure that iteratively increases the length VOLUME 8, 2020  of the filter until the first and last coefficients become lower than a user specified threshold as can be seen in algorithm 4. Frequency and impulse response of fractional Laplacian of various α's is presented in Fig. 3

5) NUMERICAL EXAMPLES
We demonstrate the use of the fractional trend filter in 1D and 2D settings. In the 1D case, we generate a synthetic time-series data with a piece-wise linear trend and additive white Gaussian noise. In Fig. 5 and Fig. 6 synthesized data, synthesized trend and filters' results are reported for three values of α ∈ {1, 0.6, 0.2}. From Fig. 5 and Fig. 6 we notice that different fractional orders of trend filtering result in trend estimates that belong to different function families. In the case of α = 1 and 1 fractional trend filter, which corresponds to the original 1 trend filter, the trend estimates are close to piece-wise linear but that is not the case for α = 0.6 or 0.2.
In the 2D case, we test fractional trend filter on the image denoising task. We conducted an experiment on six gray-scale images (cameraman, house, lena, peppers, pirate and blonde woman) shown in Fig. 7. We start with a ground truth imageũ then we add noise to it with three noise levels σ ∈ {15, 25, 50} forming an image b and we experiment with various values of α ∈ {0.1 . . . 1} with the goal of recovering an imageû as close as possible to ground truth imageũ. Each experiment was run for 10 times (noise instances) and we are reporting the average performance across the 10 runs. Performance is measured in the form of mean squared error: In the case of fractional 1 trend filtering, the results can be found in Table (1). We notice a pattern in the results that for higher noise levels, better reconstruction is achievable with smaller α. It is important to note here that the best regularization parameter λ was chosen using an exhaustive search.
In the case of fractional 2 trend filtering, experimental results can be found in Table (2). Similar procedure to 1 case was conducted. The results here again show that there is a benefit in using a fractional Laplacian for higher levels of noise.
In some cases, such as the ''Cameraman'' image and the ''House'' image, it turns out that the setting α = 1 (corresponding to the Laplacian operator) leads to the best results. We believe such results should not be regarded as a weakness of the fractional Laplacian operator. Instead, we believe that this is an advantage of the fractional Laplacian operator which includes the Laplacian operator as a special case. Compared with the Laplacian operator, the fractional Laplacian operator permits the user to ''tune'' the parameter α to achieve the desired result. As such, in other cases we have achieved better results by tuning α.
To get an idea about the difference in performance between the presented fractional trend filters and the state-of-the-art in image denoising, we have compared with BM3D [95] provided with the underlying noise standard deviation (not an estimate) in Table (3).

B. IMAGE PROCESSING APPLICATIONS 1) 1D FILTERING
To get a better idea about the effect of using the proposed fractional Laplacian, we conduct a comparison on a test signal from the Wavelet Toolbox in MATLAB named blocks. The choice of this signal was made because; it is a piecewise constant signal with sharp edges, allowing for clearer demonstration of the impact of the filters on edges. Results of linear filtering in Fig. 8 demonstrate the long memory effect. The fractional Laplacian with smaller values of α have longer memory effect than the standard Laplacian.

2) IMAGE SHARPENING
One use case for fractional Laplacian is to increase the image sharpness as follows: where γ is an amplification factor and H α is the 2D isotropic fractional Laplacian kernel in Fig. 9 formed by adding two kernels similar to the ones in Fig. 2.   Fig. 10 is a comparison of image sharpening performance between the standard Laplacian (α = 1), the fractional Laplacian (α = 0.5) and the guided filter [96]. The results in Fig. 10 clearly demonstrate that a sharper image was achieved using a fractional Laplacian with α = 0.5 than  the standard Laplacian (α = 1), this results is expected as the fractional-order Laplacian captures more information than the standard Laplacian as can be seen in the impulse responses of Laplacian filters in Fig. 3. Sharpening using the guided filter is presented as a comparison with a 89700 VOLUME 8, 2020  state-of-the-art edge-aware filter. Recall that the fractional Laplacian is a linear operator while the guided filter is nonlinear.

3) EDGE DETECTION
The literature on edge detection is rich and it is beyond the scope of this study to list and compare with all techniques in the literature however, we are examining the potential use of the proposed fractional Laplacian in this task. We begin with the traditional Marr-Hildreth edge detector [97]. A grayscale image is first smoothed with a Gaussian filter to reduce the impact of noise and make the detection more robust. This is followed by Laplacian filtering step to find the edges. Filtering an image with the Laplacian kernel results in zerocrossings where edges are potentially located. Next, slopes at the zero-crossings are computed and finally a threshold T is applied to keep significant edges only. The algorithm is summarized in Fig. 11.  We generalize Marr-Hildreth edge detector by replacing the Laplacian by its fractional-order generalization. Fig. 12 illustrates the results that can be achieved using various values of α. We notice in Fig. 12 that more robust VOLUME 8, 2020  Edge detectors based on the second order derivative are known to be sensitive to noise [98] which begs the question: does the fractional Laplacian suffer from sensitivity to noise as well? To this end, we used a synthetic image with various forms of edges that commonly exist in natural images, then we added Gaussian noise to it and finally processed it with the fractional Marr-Hildreth edge detector in Fig. 13. From Fig. 13 we see that; the fractional Laplacian is less sensitive to noise as expected because it has more coefficients (memory property) than the standard Laplacian.

4) SHOCK FILTERING
Shock filtering was initially proposed by Osher and Rudin [99] for image enhancement but, the technique received a lot of interest from researchers. Shock filters are formulated as PDEs that are evolved over time to come up with the filtered image which is characterized to be piecewise constant. The basic formulation of a shock filter is as follows: A more robust version utilizes a smoothing operator such as a Gaussian kernel, and is formulated as follows: G σ represents a two-dimensional Gaussian filter with σ being a smoothness parameter and is set to 1 in all experiments in this work.
We extend the proceeding model to fractional order which can be solved using the explicit scheme 89702 VOLUME 8, 2020 This extra parameter α gives more control over the filtering effect. In Fig. 14, we compare the fractional (α = 0.6) and the standard Laplacian (α = 1). The fractional Laplacian produces better object segmentation effect with sharper edges than the case with the standard Laplacian. To give a sense of what is happening, and because the shock filter results in piece-wise constant images, we compute a normalized anisotropic TV-L2 cost at every iteration. The TV-L2 is: where N is the number of pixels in the image u. Fig. 15 demonstrates the measure TV-L2 for different runs of the shock filter at different values of α. It is clear from Fig. 15 that the shock filter might not converge. To avoid the shock filtered result from diverging, we use TV-L2 as a stopping criterion as in Algorithm 5.

Algorithm 5 Fractional Shock Filter With Stopping Criterion
Input: α, I input Output: u n+1 1 threshold = 10 −3 ; 2 u 0 = I input ; 3 while TV-L2 (37) > threshold do 4 u n+1 using (36) 5 n ← n + 1 6 end In Fig. 16, the impact of using different values of α in the proposed shock filter is presented for different images where, the number of iterations is determined by the stopping criterion.
One thing to notice in Fig. 16 is that the fractional shock filtering results in more abstract images and sharper edges. To further validate this observation, we plot in Fig. 17 a scan line from the three channels (RGB) for the four images in Fig. 16. In each subplot of Fig. 16, we have scan-lines of two cases to facilitate visual comparison.

5) α SCALE-SPACE
The fractional Laplacian could be used as replacement for the standard Laplacian in the heat equation, this is known as α scale-space model [100]: This leads to different diffusion effects and rates. The implementation of this diffusion can be carried out efficiently in the DCT domain. First, we write the explicit scheme [3]: Consequently, the filtered signal at iteration n becomes: A demonstration of the effect of α scale-space is presented in Fig. 18. It is important to stress that this is an iterative linear filtering process. The filters are characterized as having lowpass response. The authors of [73] have explored combining multiple fractions of the Laplacian implemented in the Fourier transform domain.

IV. SUMMARY AND CONCLUSION
In this article, we have presented a new technique for constructing a fractional Laplacian using the DCT transform.
The proposed operator is a matrix operator which avoids the need for discretization, a necessary step for implementation, typically done in DSP-based constructions. We have also studied the trend filter and provided a new DCT-based implementation for it, which we used later in generalizing trend filters to the fractional-order. The proposed fractional Laplacian allowed us to make the generalization of another version of the trend filter namely the 1 trend filter to fractional-order.
To test the efficiency of the new operator, we have incorporated it in five applications that traditionally relied on the Laplacian operator. Firstly, we used the proposed fractional trend filters for image denoising and showed that for higher noise levels, the fractional-order Laplacians tend to produce better results than the standard Laplacian. Secondly, we used the proposed operator in image sharpening and stronger sharpening effect was achieved compared to the standard Laplacian. Thirdly, Marr-Hildreth scheme for edge detection was generalized and more robust edge detection was demonstrated. Fourthly, we showed that the use of fractional-order Laplacian in shock filtering resulted in better segmented images. Finally, we incorporated the proposed fractional Laplacian in the α scale-space scheme and showed that the fractional Laplacian resulted in faster diffusion.
Finally, we reiterate; the aim of this paper is to present an alternative way to define the fractional Laplacian operator, which, to our best knowledge, has not been published before. We demonstrate its potential applications in solving a wide range of problems from data modeling to image processing. We observe from the simulations, the performance is sometimes sub-optimal. This is a result of applying the fractional Laplacian operator in a direct and non-sophisticated fashion. Further improvement can be achieved by using the fractional Laplacian operator as a building block in some successful algorithms.