A Variable Order Fractional Differential-Based Texture Enhancement Algorithm with Application in Medical Imaging

Texture enhancement is one of the most important techniques in digital image processing and plays an essential role in medical imaging since textures discriminate information. Most image texture enhancement techniques use classical integral order differential mask operators or fractional differential mask operators using fixed fractional order. These masks can produce excessive enhancement of low spatial frequency content, insufficient enhancement of large spatial frequency content, and retention of high spatial frequency noise. To improve upon existing approaches of texture enhancement, we derive an improved Variable Order Fractional Centered Difference (VOFCD) scheme which dynamically adjusts the fractional differential order instead of fixing it. The new VOFCD technique is based on the second order Riesz fractional differential operator using a Lagrange 3-point interpolation formula, for both grey scale and colour image enhancement. We then use this method to enhance photographs and a set of medical images related to patients with stroke and Parkinson’s disease. The experiments show that our improved fractional differential mask has a higher signal to noise ratio value than the other fractional differential mask operators. Based on the corresponding quantitative analysis we conclude that the new method offers a superior texture enhancement over existing methods.


Introduction
Texture plays an important role in the identification of regions of interest in an image, hence texture enhancement is an essential component in digital image processing [1]. The aims of texture enhancement are to improve the overall visual effect of the image through purposefully emphasizing local or whole characteristic features and impair characteristics of little interest. The quality of the image, especially the texture, is critical in the clinical diagnosis and identification of pathology based on medical images. according to the image local statistics and structural features [29]. Results were provided for the Lena image only, and a rigorous quantitative analysis of their findings was not performed. They did show that their method had a higher signal-to-noise ratio and superior image quality than the traditional fractional differential operator and classical integral order differential mask operators. In view of these findings, it appears that variable fractional order calculus provides benefits over fixed order methods. Hence, our work aims to further develop variable fractional order methods in the context of image enhancements.
In our previous work we showed that the use of the Riesz fractional differential operator instead of the Grünwald-Letnikov definition results in higher accuracy due to the positive benefits of using a symmetric second order instead of a one sided first order fractional operator [8,26]. In view of these findings, we changed the YiFeiPU-2 method [25] by replacing the Grünwald-Letnikov definition with the Riesz fractional differential operator. Moreover, we incorporated variable fractional order [29] into the algorithm to estimate the fractal dimension of a local region instead of using fixed fractional order differentiation. The outcome is an image enhancement method applicable across a range of image types that adapts the fractional order of the differential to local features. As such, we achieve larger flexibility by being able to change the weights used in the reconstruction algorithm. The key differences of our method with previous work [8,26] are that we apply a Lagrange 3-point interpolation formula on the Riesz fractional differential operator and incorporate variable fractional order into the algorithm. Key differences of the method with Pu et al. [25] are that we replace the Grünwald-Letnikov definition by the Riesz fractional differential operator and incorporate variable fractional order into the algorithm. Key differences of the method with Huang et al. [29] are that we replace the Grünwald-Letnikov definition by the Riesz fractional differential operator and apply a Lagrange 3-point interpolation formula on the Riesz fractional differential operator to achieve different effects based on the choice of the weights used in the reconstruction. This paper is organized as follows. Firstly, the mathematical preliminaries used throughout the paper are introduced. Secondly, we give the theoretical analysis for our improved fractional differential mask, termed VOFCD, based on the second order Riesz fractional differential operator using a Lagrange 3-point interpolation formula. Then, the VOFCD method is presented. Furthermore, we present an overview of the data acquisition process and give an outline of the algorithm for enhancing the quality of a grey-level image. Finally, we compare our algorithm with Pu et al.'s best algorithm (YiFeiPU-2) [25] and Yu et al.'s algorithm (FCD-1) [8,26] when applied to medical images of the human brain.

Preliminary Knowledge
Here we outline the important definitions used throughout the paper. Definition 1.
The v-order Grünwald-Letnikov based fractional derivative GL D v x sðx; yÞ with respect to x for the finite interval x 2 [a, X] can be expressed by [30,31]: where v is any real number. The v-order Grünwald-Letnikov based fractional derivative GL D v y sðx; yÞ with respect to y can be defined in a similar manner.

Definition 2.
The v-order (0 < v 2) Riesz fractional derivative @ v sðx;yÞ @jxj v with respect to x for the infinite interval −1 < x < +1 is defined as [8]: The Riesz fractional derivative @ v sðx;yÞ @jyj v of order v (0 < v 2) with respect to y can be defined in a similar manner.

The Fractional Differential Mask
Utilizing the second order fractional centered difference scheme [8,26,32], we discretize the Riesz fractional derivative @ v s(x, y)/@jxj v (0 < v 2) with respect to x with step h as: Assume that r = x + (vh/2) − kh, then r 2 [x − kh, x + h − kh]. Using a Lagrange 3-point interpolation formula for the three neighboring nodes s(x + h − kh, y), s(x − kh, y), and s(x − h − kh, y), we have sðr; yÞ ffi ðr À x þ khÞðr À x þ h þ khÞ 2h 2 sðx þ h À kh; yÞ Noting that r = x + (vh/2) − kh, we then obtain Compared with s(x − kh, y) in Eq (5), s(x + (vh/2) − kh, y) in Eq (7) is a linear combination of the neighboring nodes, which implies that s(x + (vh/2) − kh, y) contains more information in its neighborhood and leads to richer texture details. Thus, replacing s(x − kh, y) in Eq (5) with s(x + (vh/2) − kh, y), and substituting Eq (7) into Eq (5), we obtain a new Riesz fractional differential with respect to x as By first noting that for 0 < v < 1, GðtÞGð1 À tÞ ¼ p we can rewrite Eq (8) in the form where for k = ±1, ±2, Á Á Á. It can be seen from Eq (10) that the Riesz fractional derivative can be discretized into two parts, one located on the positive x-axis and the other on the negative x-axis. From this formulation the Riesz fractional mask can be obtained.
In the context of medical images, the authors in [8,25,26] discuss the biggest variable of the grey level being limited, and the shortest distance for a change in the grey level image must be at an adjacent pixel. Therefore, the pixel signal is used to measure the duration of a two-dimensional digital image s(x, y) with respect to the two variables x and y. Here, the duration is the dimension of the image matrix assuming that the duration of x and y is For a two-dimensional digital image s(x, y) at pixel signal (x 1 , y 1 ) on the positive x-axis with region [0, x 1 + h], the N + 2 pixels are s N ðx 1 ; y 1 Þ ¼ sð0; y 1 Þ; Á Á Á ; s k ðx 1 ; y 1 Þ ¼ sðx 1 À kh; y 1 Þ; Á Á Á ; s 0 ðx 1 ; y 1 Þ ¼ sðx 1 ; y 1 Þ; s À1 ðx 1 ; y 1 Þ ¼ sðx 1 þ h; y 1 Þ: After truncation, the anterior n + 2 approximate fractional centered difference of the Riesz fractional differential with order 0 < v < 1on the positive x-axis is: Similarly, the anterior n + 2 approximate fractional centered difference of the Riesz fractional differential with order 0 < v < 1on the negative x-axis is: To obtain the fractional differential masks for the eight symmetric directions and make them rotationally invariant, we implement eight fractional differential masks positioned respectively on the negative x-axis, positive x-axis, negative y-axis, positive y-axis, left downward diagonal, right upward diagonal, left upward diagonal, and right downward diagonal. They are correspondingly denoted by W l (l = 1, 2, . . ., 8) (see Fig 1).
In Fig 1, C s 0 is the mask coefficient associated with the pixel of interest. When n = 2m − 1, a (2m + 1) × (2m + 1) fractional differential mask is implemented. To ensure that the fractional differential mask remains symmetric and the centre of the mask aligns with a pixel, in general, n is taken as an odd positive integer.
Digital image processing is based on direct processing for discrete pixels, and the algorithm normally adopts an airspace filtering scheme whose principle is to move the mask pixel by pixel [25]. Therefore, there are separate algorithms for grey and colour image fractional differential masks.
Next, we deduce the following fractional differential algorithm, VOFCD, based on the Riesz fractional differential operator Eq (2). To treat the N x × N y digital grey-level image s(x, y), we perform a convolution filter on the above eight directions in the (2m + 1) × (2m + 1) masks, and propose that the eight fractional differential masks are computed by what we refer to as the VOFCD operator: where l 1 = 1, 2, 3, 4, l 2 = 5, 6, 7, 8, and Thus, we have the digital grey-level images s I (x, y)as where C s k is the mask coefficient given in Eq (17). As for the digital colour image, the algorithm is similar to that for a grey-level image, however the RGB components use the fractional differential respectively. When 0 < v < 1, we implement the fractional differential mask respectively on the eight symmetric directions using what we call the VOFCD operator, having the same structure as YiFeiPU-2 in [25] but with different coefficients. The mask coefficients of the VOFCD operator are given by . . .
which ensures that the fractional differential operator VOFCD produces a sparse matrix having dimension n + 2. Moreover, all the coefficients depend on the fractional differential order v. It can also be proven that the sum of the coefficients is nonzero, which is the main difference between the fractional differential mask and the integral version.

Variable Fractional Differential Order
The average gradient and information entropy are widely used for quantitative analysis of images [23,25,27]. The average gradient reflects the clarity of the image, and expresses contrast due to small details. It can be used to measure the spatial resolution of the image, i.e., a larger average gradient means a higher spatial resolution [33][34][35]. Generally speaking, a large gradient value is likely to correspond to regions of edges, margins and textures within images, hence a larger fractional order is required to enhance textural details [8,25,26]. On the contrary, a smaller value of the average gradient is likely to correspond to smooth image regions, hence a smaller fractional order is needed to enhance image textures. However, the gradient is quite large for both boundaries and noise. Therefore, a suitable method is needed to identify the margin and noise. Note that the margin of the object is continuously smooth, and it is easy to find another margin point that has a similar magnitude of gradient as the margin of interest. However, the noise is random, and it is difficult to find any noise around the region of interest having a similar magnitude of gradient. Hence, margin or noise can be identified through the following process: where rs ij is the gradient magnitude associated with the pixel of interest s(x i , y j ), rN(s ij ) is the gradient magnitude associated with the pixel around the region of interest s(x i , y j ), and T is the threshold value.
The information entropy evaluates the average information included in the image and reflects the detail information of the image [23,25,27,[36][37][38], that is small in the smooth area and large in the rich texture area of the image. Image entropy is calculated as: where P(x i , y j ) is the probability mass function of s(x i , y j ). The local image roughness is the relative offset measurement of the image pixel grey scale value [36,39,40], that is also small in the smooth area and large in the rich texture area of the image. The local roughness of an image is calculated as where σ is the local variance of the pixel signals of the image. The order of the fractional differential is not only related to the magnitude of the gradient, but also impacted by the image local statistics, namely variance and entropy [23,25,27]. Thus, using a weighted summation of three parameters to reflect local image information, we have: where 0 k 1 , k 2 , k 3 1 and k 1 + k 2 + k 3 = 1.
As suggested by Huang et al. [29], utilizing the exponential function property, we can obtain the variable fractional differential order function as: where α and β are regularization parameters.

Data Acquisition and Algorithm Development
In this section, we outline the implementation of our algorithm to enhance textures in grey level and colour images. We first present the details of the data acquisition. Then, we demonstrate the algorithm development.

Data Acquisition
The research project was approved by the Human Research Ethics Committee of Queensland Health, Brisbane, Australia. Participants provided written consent prior to participating in the study, the process of which was approved by the ethics committee. Furthermore, all data sets were de-identified for the project.
The first data set of the stroke patients was acquired at the Royal Brisbane and Women's Hospital using the 3T Siemens MRI scanner. Patients admitted with a clinical diagnosis of acute ischemic stroke were recruited between May 2011 and April 2012. Patients underwent an MRI examination at admission from which two MRI datasets were randomly selected for this study. Susceptibility weighted magnetic resonance images were acquired on a 3T Siemens Trio human scanner running the Syngo proprietary software housed at the hospital. Data acquisition was performed using the Syngo SWI sequence with the following parameters: matrix size = 224 × 256, repetition time (TR) = 200ms, echo time (TE) = 20ms, flip angle = 15 o , bandwidth = 120Hz per pixel, in-plane resolution, 1mm × 1mm slice thickness and separation = 2mm and number of slices = 72. Magnitude, phase and susceptibility images were saved separately, and only the magnitude images were used in this study. The original image slices 32 to 47 from the first stroke patient show the stroke in the left hemisphere of the brain, as for example, can be seen in Figs 2(a) and 3(a). For the second patient, the stroke was visible in slices 42 to 61 to various extents.
The second data set is for MRI data acquisition-3T in vivo for a patient with Parkinson's disease. The data was obtained from St Andrew's War Memorial Hospital, Brisbane, Australia. The pre-surgery data we use here is from a patient diagnosed with Parkinson's disease. The Fractional Calculus in Medical Imaging brain diffusion tensor MR images are acquired using a GE Medical System (SIGNA 3T) scanner. All image matrix sizes are 256 × 256. Pre surgery data is acquired using an echo time of 93.7 ms and repetition time of 7s. Twenty four interleaved, 5 mm thick slices are acquired in the horizontal plane perpendicular to the coronal and sagittal planes, using the multislice mode. Diffusion sensitization is performed along 35 different diffusion gradient orientations using a diffusion weighting of b = 1000 sec/mm 2 (b is the degree of diffusion sensitization defined by the amplitude and the time course of the magnetic field gradient pulses used to encode molecular diffusion displacements [41]). A reference image without diffusion weighting (b = 0 sec/mm 2 ) not illustrated here was also acquired.

Image Evaluation
The regulation parameters are set to α = 1 and β = 1.7, which gives the variable fractional differential order function as:  where A is the root mean square amplitude. For MRI data, the image background was used to compute A noise . Images are evaluated based on the intensity histogram Entropy (Ent), Standard Deviation (STD), and Mean Absolute Difference Coefficient (MADC) [42]. The higher the value of Ent, the more visual information the image contains, and an increase in Ent means an increase in information and therefore an improvement in image quality. The STD demonstrates how much the image intensities deviate from the expected value, and a variation of STD measures how much corresponding points vary across the source and images. The MADC is used to measure image clarity, such as activity, and is defined as the mean of the sum of difference values over rows and columns of the image pixel intensities [42,43]: where

Simulation Study
The algorithm given in Table 1 provides the process for enhancing the quality of a grey-level image using the symmetric Riesz variable fractional order substitution.

Results and Discussion
In this section, we compare our algorithm to YiFeiPU-2 [25] and FCD-1 [8,26], and apply the method to a set of medical images.

Evaluation of Weights
After testing weights k 1 , k 2 and k 3 ranging from 0 to 1 with step 0.01 subject to the condition k 1 + k 2 + k 3 = 1, Figs 4-6 show the convergence curves of SNR, Ent, STD and MADC with different mask sizes in image slice 35 for the first stroke patient (see Fig 2(a)). Using a 5 × 5 mask, Tables 2 and 3 show different combinations of k 1 , k 2 and k 3 for obtaining the largest and lowest SNR, Ent, STD and MADC. According to the different combinations of k 1 , k 2 and k 3 from Tables 2 and 3, Fig 7 shows the corresponding image slice. Table 1. Algorithm for enhancing the quality of a grey-level image.
Input: Read original grey image, and add Gaussian noise with mean 0 and variance 0.01.

Output: s(x, y)
Choose m = 2, n = 3 and mask 5 × 5 Compute the variable fractional differential order using Eq (23) Compute the mask coefficients C s k using Eq (17)   In a similar manner to the data exhibited for the first stroke patient, Figs 8-10 show the convergence curves for SNR, Ent, STD and MADC with different mask sizes for image slice 51 from the second stroke patient. Figs 11-13 show the corresponding results from clinical data for a patient diagnosed with Parkinson's disease prior to surgery.
Additional experiments were performed on different image slices of stroke patients, the results of which are provided as Supporting Information.
We have shown that the outlined approach can be applied to different data, and the weights can be manipulated to achieve the highest or lowest SNR, Ent, STD or MADC.

Qualitative Analysis
Without loss of generality, we only show results using the 5 × 5 mask applied to the first stroke patient image using two particular combinations of k 1 , k 2 and k 3 . The results of the second stroke patient are provided as Supporting Information. The parameter set (0.45, 0.01, 0.54) obtained the largest SNR in slice 35 depicted in Fig 7, whilst (0.01, 0.01, 0.98) produced the largest MADC. These two choices of parameters allow us to evaluate image SNR versus contrast.     . We now present clinical data used to test our methods on human brain images from first patient diagnosed with stroke. In Fig 2, we compare VOFCD with YiFeiPU-2 and FCD-1. Fig 3 shows the comparison of texture details in the region of interest between the original image slice 35 and its differential using the YiFeiPU-2, FCD-1 and VOFCD methods. It can be seen from Figs 2 and 3 that the image resulting from our method (VOFCD) is qualitatively better than the images constructed using the other methods. We note in particular that the visual effect offered by VOFCD is better than that of YiFeiPU-2 and FCD-1 with fractional order v = 0.5 when using the same mask dimensions. This choice of v allows us to make consistent comparisons Through the use of histograms, Fig 14 shows the comparison of texture details between the original image slices 40, 44, 45, and 46 of first stroke patient, and their fractional differential using VOFCD with a 5 × 5 mask. Again, we can see that VOFCD has enhanced the textural details in these images. Fig 15 provides a comparison of SNR between the original image slices of first stroke patient and the fractional differential with mask 5 × 5 for VOFCD, YiFeiPU-2 and FCD-1 with v = 0.5. From Fig 15, it can be seen that VOFCD produces the largest values of SNR, which implies a superior texture enhancement over the other methods.
We now present clinical data to test our methods on a human brain image from a patient diagnosed with Parkinson's disease prior to surgery.
Applying the fractional differential to the three elements in the HSI colour space, respectively, and then reverting to RGB colour space, one can obtain a colour image without distortion. Fig 12 shows that the parameter set (0.25, 0.74.0.01) obtained the largest SNR for an original fractional anisotropy weighted orientation map from a patient diagnosed with Parkinson's disease prior to surgery with mask 5 × 5, and the parameter set (0.01, 0.01, 0.98) produced Fractional Calculus in Medical Imaging the largest MADC. Fig 16 shows the comparison of texture details between an original fractional anisotropy weighted orientation map and its fractional differential using YiFeiPU-2, FCD-1 and VOFCD using weights (0.25, 0.74.0.01). Table 4 provides the comparison of SNR for a region of interest of a fractional anisotropy weighted colour orientation map between VOFCD using (0.25, 0.74.0.01) and YiFeiPU-2 and FCD-1 with v = 0.5. We can conclude from Table 4 that VOFCD has the largest SNR in comparison to the other fractional order differential methods. The value of this measure is consistent with the VOFCD image shown in Fig 16, wherein features appear to be smoother than in the other images.  Table 5 show the results using weights (0.01, 0.01, 0.98). With this choice of weights, the SNR is reduced however MADC increases. This finding is consistent with a likely increase in image contrast. Hence, the choice of weights can be used to influence or trade-off SNR against contrast. It is unlikely that a fixed choice of weight can be applied to any image, since images not only vary in modality but also have different types of textures and contrast within them.  using VOFCD, the values of STD are larger than the other fractional differential methods. The increase in STD is due to increased signal intensities across the entire image. Furthermore, VOFCD has a smaller value of MADC than the other fractional differential methods, which implies a reduction in noise in VOFCD reconstructed images.
Based on a human brain image from a patient diagnosed with Parkinson's disease prior to surgery, Table 6 provides a comparison of the relevant quantitative analysis for a region of interest of a fractional anisotropy weighted colour orientation map between VOFCD using weights (0.25, 0.74.0.01) and YiFeiPU-2 and FCD-1 with v = 0.5. We conclude from Table 6 that VOFCD has the largest Ent value, and smallest STD and MADC values in comparison to the other fractional order differential methods. Again, the values of these measures are consistent with the VOFCD image shown in Fig 16, wherein features appear to be smoother than in the other images.
In a similar manner to the first set of weights, Fig 20 and Table 7 show the results using weights (0.01, 0.01, 0.98). We have established that the use of the fractional differential not only nonlinearly preserves the contour features in smooth areas, but also maintains a high frequency edge feature in areas where the grey scale has obvious changes. It also preserves high frequency characteristics of texture detail in those areas where the grey scale does not change considerably. Furthermore, VOFCD has a higher SNR than the other fractional differential mask operators, and our quantitative analysis verified that VOFCD leads to superior texture enhancement.  Using our method it is possible to achieve different effects based on the choice of the weights used in the reconstruction (i.e.k 1 , k 2 and k 3 ). For example, we showed that SNR can be maximised through appropriate choices of the weights. Similarly, we found weights that maximised MADC. Therefore, the choices of the weights can be optimised for the application. Nonetheless, it is difficult to establish a set of weights that can be applied robustly across all applications, since particular applications may rely on suppression of noise (i.e. SNR improvements) and others on differentiation of structures within images (i.e. MADC improvements).

Conclusion
For grey scale and colour image enhancement, we derived an improved fractional differential algorithm, VOFCD, based on the second order Riesz fractional differential operator using a Lagrange 3-point interpolation formula. We estimated the fractal dimension of a local region instead of using fixed fractional order differentiation. The experiments showed that the use of VOFCD results in higher SNR than the existing methods of YiFeiPU-2 and FCD-1, which implies superior texture enhancement of medical images. In addition, VOFCD produces qualitatively better results than YiFeiPU-2 and FCD-1 with variable fractional order and same mask dimensions. A quantitative analysis was conducted to verify that VOFCD produces superior texture enhancement in comparison to the other methods considered. This may be helpful in clinical diagnosis and monitoring, and as future work, further analysis of this data will be carried out in conjunction with medical specialists. The optimal choice of reasonable parameters to obtain the variable fractional differential order remains challenging. In future research, we plan to perform an evaluation of how the weighted parameters and regularisation parameters affect the performance of the variable fractional differential order algorithm.    We also thank the reviewers for their constructive comments and suggestions.

Author Contributions
Conceived and designed the experiments: QY FL IT. Performed the experiments: QY. Analyzed the data: QY VV FL IT. Contributed reagents/materials/analysis tools: QY VV FL IT. Wrote the paper: QY VV FL IT.