# Variable-Order Fractional Diffusion Model-Based Medical Image Denoising

Authors: A. Abirami; P. Prakash; Y.-K. Ma
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Fractional differential models are playing a vital role in many applications such as diffusion, probability potential theory, and scattering theory. In this study, the variable-order space and time fractional diffusion model is employed for denoising the medical images. The finite difference approach is implemented to find the numerical solution of the proposed model. Convergence and stability of the numerical method are presented. The experimental outcomes of the variable-order model are analyzed with those of the fractional and integer-order diffusion models. It was noticed that the peak signal-to-noise ratio (PSNR) value is increased considerably for the proposed model.


---

## Body

## 1. Introduction

The field of fractional differential equations has drawn immense consideration towards theoretical [1–3] and applied research studies [4–6]. Fractional differential equations are found to be an effective tool used in certain mathematical models such as hydrology [7], finance [8], physics [9], and signal and image processing [10–14]. Various theories of fractional integrals and derivatives were developed by many authors, for instance, Riemann-Liouville [2], Gru¨nwald [15], Caputo [16], and Riesz [17].Many researchers are focusing on image denoising based on total variation, wavelet transform, bilateral filter, histogram of gradient, fractionalsinCα, anisotropic diffusion filter, and primal-dual algorithm [18–25]. However, the drawback of using a total variation, second-order, and fourth-order-based image denoising models suffered from the staircase effect, too much of smoothness and preserving in discontinues. The integer-order fractional derivatives are not efficient in specifying some complex diffusion processes. So, our work aims to further development of the variable-order fractional model in the context of image denoising to overcome the above issues.Nowadays, variable-order fractional calculus is especially acknowledged as a helpful and hopeful approach in the modeling of a diffusion process. The pioneering work of variable-order operators can be traced. Chen et al. [26] introduced Caputo-fractional derivatives for removing noisy signal using wavelet transform. In recent years, a finite difference technique for solving multivariable-order fractional integrodifferential equations is introduced using Bernstein basis functions together with the Newton-Cotes collocation points in [27]. Furthermore, they have found a numerical scheme to solve variable-order diffusion-wave and differential equation. Also, in [28], Guo et al. proposed three-dimensional fractional total variation under tensor algebra-based model for 3D image denoising. Recently, in [29], Garrappa et al. investigated the main properties of the emerging variable-order operators and discussed some practical applications of the variable-order Scarpi integral and derivative. Gu et al. [30] proposed an unconditionally stable implicit difference scheme for solving generalized time-space fractional diffusion equations with variable coefficients with numerical scheme that utilizes the L1-type formula for the generalized Caputo-fractional derivative in time discretization and the second-order weighted and shifted Gru¨nwald difference formula in spatial discretization. Also, Fang et al. [31] developed a fast finite difference method for solving a class of variable-order time fractional diffusion equations. In [32], Gu and Wu considered a class of Volterra partial integrodifferential problems and proposed a novel iterative algorithm for parallel-in-time pattern (PinT) computation.In this research, we study the variable-order fractional diffusion model for medical image denoising using the Caputo finite difference scheme for the proposed problem. The experiments demonstrate the advantage of the variable-order fractional model that achieves the highest PSNR value, and it indicates the quality enhancement of the medical images. This research is organized as follows. In Section2, few types of variable-order fractional derivatives are explored. In Section 3, the discrete approximation of the proposed model is developed. In Section 4, the stability and convergence of this model are discussed. In Section 5, the numerical results are compared with fractional and integer-order diffusion models. This study concludes in Section 6 with directions for further research.


## 2. Variable-Order Fractional Operators

The research on fractional calculus is pursued over a long time in different disciplines such as biomedical, computational biology, economics, and control engineering. The classical calculus extends the definition of fractional calculus, where the orders need not be positive integers. On the other hand, the variable-order calculus is a hereditary extension of the integer order calculus. In this regard, the order may use in any variable such as time and space variables or a system of other variables. The basic definitions for variable-order fractional derivatives are provided in the following sequel.Definition 1 (See [16]).
The variable-order Caputo derivative with respect to time is defined as follows:(1)cDtαx,tu=1Γn−αx,t∫0t1t−sαx,t−n+1∂nus,t∂snds,where n−1<αx,t≤n.Definition 2 (See [16]).
The variable-order Caputo derivative with respect to space is defined as follows:(2)cDxβx,tu=1Γn+1−βx,t∫0x1x−sβx,t+n−2∂n+1us,t∂sn+1ds,where n<βx,t≤n+1.


## 3. Variable-Order Fractional Diffusion Equation-Based Image Denoising

ConsiderΩ is a closed domain in IR2. The format of noisy image would be mathematically modeled as(3)u0x,y=ux,y+gx,y,where u0x,y and ux,y are the observed and clean images, respectively, gx,y is the Gaussian noise, and i,jT takes the location with a rectangular image domain Ω⊆IR2. The proposed work approximates the required clean image ux,y by solving the following variable-order fractional diffusion equation:(4)cDtαx,y,tux,y,t=cDxβx,y,tux,y,t+cDyγx,y,tux,y,t,where xL<x<xR,yL<y<yR, αx,y,t∈0,1, and βx,y,t, γx,y,t∈1,2, where uijn=uxi,yj,tn, αijn=αxi,yj,tn, βijn=βxi,yj,tn, and γijn=γxi,yj,tn. Define tn=nΔt is in the case 0≤tn≤T, Δx=h>0 as the grid size in x-direction, Δy=k>0 as the grid size in y-direction, and x,y∈Ω, where Ω is an image domain with initial and Neumann boundary conditions:(5)ux,y,0=u0x,y,0,∂u∂xxL,y,t=∂u∂xxR,y,t,∂u∂yx,yL,t=∂u∂yx,yR,t.Now, the above Neumann boundary conditions are discretized as(6)ui,Nn=ui,N−1n,uN,jn=uN−1,jn,where M×N is the dimension of the image ux,y. We use these Neumann boundary conditions to preserve continuity of the image boundary. The boundary values of the image can be obtained by reflecting the closest samples inside the image region. The solution to the proposed problem can be obtained by discretization of space and time variables. The discretization of the Caputo-type variable-order space fractional derivative can be stated as follows [33].(7)cDxβx,y,tuxi,yj,tn=h−βijnΓ3−βijn∑s=0i−1ui−s+1,jn−2ui−s,jn+ui−s−1,jns+12−βijn−s2−βijn,where Bijs=s+12−βijn−s2−βijn and Qij=h−βijn/Γ3−βijn.(8)cDyγx,y,tuxi,yj,tn=k−γijnΓ3−γijn∑s=0j−1ui,j−s+1n−2ui,j−sn+ui,j−s−1ns+12−γijn−s2−γijn,where Cijs=s+12−γijn−s2−γijn and Rij=k−γijn/Γ3−γijn.Similarly, the discretization of the Caputo-type variable-order time fractional derivative can be stated as follows [2].(9)cDtαx,y,tuxi,yj,tn=τ−αijnΓ2−αijn∑s=0nuijn+1−s−uijn−s⋅s+11−αijn−s1−αijn,where Aijs=s+11−αijn−s1−αijn and Pij=Γ2−αijn/τ−αij.Discretize equation (4) in the following form:(10)τ−αijn+1Γ2−αijn∑s=0nuijn+1−s−uijn−sAijs,=Qij∑s=0i−1ui−s+1,jn−2ui−s,jn+ui−s−1,jnBijs+Rij∑s=0j−1ui,j−s+1n−2ui,j−sn+ui,j−s−1nCijs.Now, equation (8) can be written as(11)uijn+1=uijn−∑s=1nuijn+1−s−uijn−sAijs+PijQij∑s=0i−1ui−s+1,jn−2ui−s,jn+ui−s−1,jnBijs+PijRij∑s=0j−1ui,j−s+1n−2ui,j−sn+ui,j−s−1nCijs.Equation (9) can be entered at time t=tn in the matrix form:(12)Uijn+1=TijnUijn+Fijn,where Fijn=Aijnfum−1,n−1n,xm−1,n−1,ym−1,n−1,tmn…Aijnfu1,1n,x1,1,y1,1,tmnT, Uijn=um−1,n−1n,um−2,n−2n…u1,1nT, and Tijn=anm is a matrix with the following coefficients:(13)am,n=0,wherem>n+2,PijQijBij0+PijRijCij0,wherem=n+2,1+Aijn+PijQijθBij−2Bij0+PijRijθCij−Cij0,wherem=n+1,PijQijBijn−m−2Bijn−m+1+θBijn−m+2+PijRijCijn−m−2Cijn−m+1+θCijn−m+2,wherem≤nPijnQijnBij0,wherem=1θ=0,form=2,1,otherwise.The computational solutions ofuijn+1 and uijn are described from the identified initial and symmetric boundary conditions (3) and (4). The solution uijn+1 of (8) is meant as a denoised image.The implementation of the proposed method is expressed in the following Algorithm1:Algorithm 1:
(1)
Initialize the iteration process by settingu0x,y(2)
Initializeα∈0,1 and β,γ∈1,2,u0=f, and time step τ=0.1(3)
Computeuijn+1 from (6), where i = 1, 2, …, M, j = 1, 2, …, N, for n = 1,2, …, NT do(4)
Check ifun+1−un/un≤total; then stop.(5)
Setun+1=ux,y(6)
Output displayux,y


## 4. Stability and Convergence Analysis

### 4.1. Stability Analysis

ConsiderWn+1 and Un+1 be the numerical solutions of (6) with initial values given by Wij0 and Uij0, respectively.Theorem 1.
The explicit approximation method defined by (9) to variable-order space-time diffusion equation (2) is stable, that is,(14)Wn+1−Un+1≤CW0−U0for anyn,where C is a positive number independent of h, k, and τ.Proof.
DefineWn+1−Un+1=ϵn+1. From (6),(15)ϵijn+1=Tijnϵijn+Fϵn,where(16)Fϵn=Am−1,n−1nfU1,m−1n,x1,m−1n,y1,n−1n,tm,n−Am−1,n−1nfW1,m−1n,x1,m−1n,y1,n−1n,tm,n,⋯,Am−1,n−1nfU1,1,x1,1,y1,1,tm,n−Am−1,n−1nfW1,1n,x1,1,y1,1,tm,nT≤Am−1,n−1nLm−1,n−1nϵm−1,n−1n…,A1,1nL1,1nϵ1,1n=ΔFnϵn,ΔFn=diagAm−1,n−1nLm−1,n−1n…,A1,1nL1,1nT.
We note thatLijn≤L for any i,j.
LetA¯=maxAm−1,n−1n….,A1,1n. From equation (15), Tn+ΔFnm,n≤2+A¯L. Then,(17)ϵijn+1∞≤Tn+ΔFn∞ϵijn∞.
The stability via mathematical induction is analyzed in [16]. From (1), ϵijn∞≤Cϵij0∞, where C is a constant.
From (11), ϵijn+1∞≤C2+A¯Lϵij0∞≤C1ϵij0∞, which proves that the explicit scheme is stable. If there is a perturbation in uij0, then the small change would not cause large error in numerical solution.


### 4.2. Convergence Analysis

Theorem 2.
Letun+1xi,yj,tk be the numerical solution of equations (4)–(6) at mesh point xi,yj,tk and computed by equations(5)-(6), and (8). Then, there is a positive constant C, such that |uijn−uxi,yj,tn|≤Cτ+h+k, where i=1,…,M−1,j=1,…,N−1, and k=1,…,n.Proof.
Defineeijk=uxi,yj,tk−uijk and ek=e1k,…,en−1kT. Using e0=0 and uijk=uxi,yj,tk−eijk, substitution into equation (16) leads to(18)εijn+1=Tijnεijn+Fijn,εij0=0,and Fijn≤Cτατ+h+k, where C is a positive constant C>0. Using mathematical induction, we provide the convergence analysis as follows. For n=1,(19)letεij1∞=eij1=max1≤i≤M−1,1≤j≤N−1eij1.
We haveεl1=Fl1≤Cτατ+h+k,(20)y1∞≤CB0−1τατ+h+k.
Suppose thatεij1∞≤CB−1τατ+h+k, and let(21)eK+1≤el0−∑s=1kuijn+1−s−uijn−sAijs+PijQij∑s=0i−1ui−s+1,jn−2ui−s,jn+ui−s−1,jnBijs+PijRij∑s=0j−1ui,j−s+1n−2ui,j−sn+ui,j−s−1nCijs.≤1+∑s=1kAijs+PijQij∑s=0i−1Bijs+PijRij∑s=0j−1Cijs,≤CBk−1τατ+h+k,elk+1∞≤CBk−1τατ+h+k,because(22)ltk⟶∞Bk−1kα=ltk⟶∞k−αk+11−α−k1−α=ltk⟶∞k−11+1/k1−α−1=ltk⟶∞k−11−αk−1=11−α.
Hence, there is a constant C,ek∞≤Ckατ1+α+ταh+ταk=Cτ+h+k. If kτ≤T is finite.
Thus, we see that for anyx,y, and t as h,k, and τ approach 0 in such a way that ih,jk,τ tends to x,y,t. This proves that uijn+1 converges to uijn as h,k, and τ tends to zero [34]. Hence, the conclusion follows.


## 4.1. Stability Analysis

ConsiderWn+1 and Un+1 be the numerical solutions of (6) with initial values given by Wij0 and Uij0, respectively.Theorem 1.
The explicit approximation method defined by (9) to variable-order space-time diffusion equation (2) is stable, that is,(14)Wn+1−Un+1≤CW0−U0for anyn,where C is a positive number independent of h, k, and τ.Proof.
DefineWn+1−Un+1=ϵn+1. From (6),(15)ϵijn+1=Tijnϵijn+Fϵn,where(16)Fϵn=Am−1,n−1nfU1,m−1n,x1,m−1n,y1,n−1n,tm,n−Am−1,n−1nfW1,m−1n,x1,m−1n,y1,n−1n,tm,n,⋯,Am−1,n−1nfU1,1,x1,1,y1,1,tm,n−Am−1,n−1nfW1,1n,x1,1,y1,1,tm,nT≤Am−1,n−1nLm−1,n−1nϵm−1,n−1n…,A1,1nL1,1nϵ1,1n=ΔFnϵn,ΔFn=diagAm−1,n−1nLm−1,n−1n…,A1,1nL1,1nT.
We note thatLijn≤L for any i,j.
LetA¯=maxAm−1,n−1n….,A1,1n. From equation (15), Tn+ΔFnm,n≤2+A¯L. Then,(17)ϵijn+1∞≤Tn+ΔFn∞ϵijn∞.
The stability via mathematical induction is analyzed in [16]. From (1), ϵijn∞≤Cϵij0∞, where C is a constant.
From (11), ϵijn+1∞≤C2+A¯Lϵij0∞≤C1ϵij0∞, which proves that the explicit scheme is stable. If there is a perturbation in uij0, then the small change would not cause large error in numerical solution.


## 4.2. Convergence Analysis

Theorem 2.
Letun+1xi,yj,tk be the numerical solution of equations (4)–(6) at mesh point xi,yj,tk and computed by equations(5)-(6), and (8). Then, there is a positive constant C, such that |uijn−uxi,yj,tn|≤Cτ+h+k, where i=1,…,M−1,j=1,…,N−1, and k=1,…,n.Proof.
Defineeijk=uxi,yj,tk−uijk and ek=e1k,…,en−1kT. Using e0=0 and uijk=uxi,yj,tk−eijk, substitution into equation (16) leads to(18)εijn+1=Tijnεijn+Fijn,εij0=0,and Fijn≤Cτατ+h+k, where C is a positive constant C>0. Using mathematical induction, we provide the convergence analysis as follows. For n=1,(19)letεij1∞=eij1=max1≤i≤M−1,1≤j≤N−1eij1.
We haveεl1=Fl1≤Cτατ+h+k,(20)y1∞≤CB0−1τατ+h+k.
Suppose thatεij1∞≤CB−1τατ+h+k, and let(21)eK+1≤el0−∑s=1kuijn+1−s−uijn−sAijs+PijQij∑s=0i−1ui−s+1,jn−2ui−s,jn+ui−s−1,jnBijs+PijRij∑s=0j−1ui,j−s+1n−2ui,j−sn+ui,j−s−1nCijs.≤1+∑s=1kAijs+PijQij∑s=0i−1Bijs+PijRij∑s=0j−1Cijs,≤CBk−1τατ+h+k,elk+1∞≤CBk−1τατ+h+k,because(22)ltk⟶∞Bk−1kα=ltk⟶∞k−αk+11−α−k1−α=ltk⟶∞k−11+1/k1−α−1=ltk⟶∞k−11−αk−1=11−α.
Hence, there is a constant C,ek∞≤Ckατ1+α+ταh+ταk=Cτ+h+k. If kτ≤T is finite.
Thus, we see that for anyx,y, and t as h,k, and τ approach 0 in such a way that ih,jk,τ tends to x,y,t. This proves that uijn+1 converges to uijn as h,k, and τ tends to zero [34]. Hence, the conclusion follows.


## 5. Experimental Results and Discussion

The effect of this variable-order model is estimated by calculating the numerical indicator PSNR value, which is frequently used to find the quality of the processed image and also compute high quality of restoration in image firmness. The PSNR is described via the mean square error (MSE) for two images, that is to say,u0 (corrupted image) and u (denoised image), respectively.(23)MSE=1MN∑i=1M∑j=1Nu0i,j−ui,j2,PSNR=10log10maxu0,u2MSE.The PSNR and MSE values of the denoised image are determined which act as a computable quality for comparison of the proposed model, fractional-order, and integer-order diffusion models (Table1). Figure 1 shows the relationship between β,γ, and PSNR values. Here, NT  = 250, that is, after 250 iterations, we obtained the resulting image which is very close to the input image. The computational time (CPU) is given in Table 1. To this end, we performed on the benchmark image databases where a noisy image is evolved by using the fractional-order and integer-order diffusion models and the proposed model with α (0<α≤1) fractional-order derivative would attain the spotting of the inflection point without any relocation. The parameter σ shows Gaussian noise, and each pixel in the image will be changed from its original value by a small amount. The terminating formula for the iteration is un+1−un/un≤ε. We studied ϵ=10−5 in the numerical experiment. In experimental and numerical analyses, we have used the version of MATLAB R2012a. The proposed model has been studied with four types of medical images, namely, OCT (optical coherence tomography) image, MRI (magnetic resonance image), CT (computed tomography) scan, and X-ray image. From Figures 2(e), 3(e), 4(e), and 5(e), the proposed model provides affirmative evidence for better performance in maintaining the details and edge information than the fractional and integer-order models.Table 1
The PSNR and MSE values obtained by applying different test images.
Modelsσ=10σ=15σ=20σ=25Avg. PSNRMSEIterations (sec)OCT-retina image (α=0.5 and β=γ=1.5)Integer-order diffusion31.9736.534.0531.4533.490.990.69Fractional-order diffusion38.1439.9535.2734.8137.040.670.67Proposed model41.0444.6344.0541.2942.750.440.53MRI-brain image (α=0.6 and β=γ=1.6)Integer-order diffusion24.1924.1624.1924.1524.170.750.68Fractional-order diffusion36.1436.1836.1936.1536.170.250.66Proposed model43.5043.7343.8444.1343.80.150.50CT-lung image (α=0.7 and β=γ=1.7)Integer-order diffusion26.7225.6924.9725.7325.780.270.66Fractional-order diffusion36.4937.1937.1136.5536.840.210.64Proposed model44.0944.5744.5044.2944.360.190.54X-ray-spider image (α=0.8 and β=γ=1.8)Integer-order diffusion30.8932.7232.7532.2632.150.350.66Fractional-order diffusion40.3041.0641.0741.0140.860.310.62Proposed model44.1344.0344.0143.9944.040.240.56Figure 1
Comparison of three models betweenβ, γ, and PSNR values.Figure 2
The denoising results of OCT-retina image withσ=20, α = 0.5, and β=γ=1.5. (a) Input image. (b) Corrupted image. (c) Noise revmoved image by integer order diffusion. (d) Noise revmoved image by fractional order diffusion. (e) Noise revmoved image by fractional variable order diffusion.
(a)(b)(c)(d)(e)Figure 3
The denoising results of MRI-brain image withσ=15, α = 0.6, and β=γ=1.6. (a) Input image. (b) Corrupted image. (c) Noise revmoved image by integer order diffusion. (d) Noise revmoved image by fractional order diffusion. (e) Noise revmoved image by fractional variable order diffusion.
(a)(b)(c)(d)(e)Figure 4
The denoising results of CT-lung image withσ=10, α = 0.7, and β=γ=1.7. (a) Input image. (b) Corrupted image. (c) Noise revmoved image by integer order diffusion. (d) Noise revmovedDenoised image by fractional order diffusion. (e) Noise revmoved image by fractional variable order diffusion.
(a)(b)(c)(d)(e)Figure 5
The denoising results of X-ray-spider image withσ=25, α = 0.8, and β=γ=1.8. (a) Input image. (b) Corrupted image. (c) Noise revmoved image by integer order diffusion. (d) Noise revmoved image by fractional order diffusion. (e) Noise revmoved image by fractional variable order diffusion.
(a)(b)(c)(d)(e)Figure2 (OCT image) shows us beneath the surface of the retina. The OCT is the most valuable advancement in the retinal diagnostic image, and we can greater understand the very fine changes which can be indicated abnormally. It constructs a cross-sectional view of ocular structures accurate to less than 10 microns. So, denoising these types of images is still challenging. Hence, the retinal diagnostic image (Figure 2(a)) is used for the first experiment. Figure 2(a) is a corrupted image by an additive Gaussian noise at four various intensities such as 10, 15, 20, and 25, as shown in Figure 2(b). The denoised image 1(e) by the presented model is added clear, and it has a larger PSNR value than Figure 2(c) and Figure 2(d). Also, Figure 2(d) shows the information of retina that is critical in picking up the most subtle changes same from Figure 2(a). Therefore, the proposed model can be used for enhanced eye examination and in addition to all the normal tests.Figure3 (MRI image) shows a brain tumor image. Recently, biomedical images are taking an important role in identifying anatomy and the physiological process of the body in both health and disease. This MRI image has proven successful for the diagnosis of all parts of the body including cancer, stroke, heart and vascular disease, breast disease, and musculoskeletal disorders. Detection of tumors from MRI data is tedious for physicians and challenging for computers in all medical disciplines. MRI scans detect a large number of lesions and define the location more readily, and they are also the best at detecting spread to the meninges, the lining around the brain (and spinal cord). Figure 3(a) is the brain image affected by the tumor. Figure 3(b) is affected by additive Gaussian noise at four different values such as 10, 15, 20, and 25. Figures 3(c) and 3(d) have lost the information as well as the appearance of the tumor is not clear. Figure 3(e) shows original information without any loss and can be used for diagnosis and surgical purposes.Figure4 (CT image) shows a lung cancer image. This CT image is a diagnostic image that is used to test and construct detailed images of internal organs, soft tissue, bones, and blood vessels. A chest CT scan can assist determine the cause of lung symptoms such as shortness of breathing or chest pain or check lung problems such as a tumor and excess fluid around the lung. As shown in Figure 4(a), lung image affected by cancer and CT lung screening is capable of detecting lung nodules as small as 2 or 3 millimeters. Figure 4(b) is corrupted by additive Gaussian noise at four different levels such as 10, 15, 20, and 25. In Figure 4(c) and Figure 4(d), small nodules and affected areas are not visible. Figure 4(e) shows malignant tumors when they are still small and can be removed before the disease spreads to other areas of the body.Figure5 (X-ray image) shows a spider image. X-ray image has enhanced an important method for visualizing cellular and histological structures in a broad range of biological and medical studies. Figure 5(a) is a high-resolution X-ray image of a spider. Figure 5(b) has been perturbed by additive Gaussian noise at 10, 15, 20, and 25. Figures 5(c) and 5(d) are not making possible to see smaller details. But, Figure 5(e) is making it possible to see similar details of Figure 5(a).


## 6. Conclusion

In this work, the convergence and stability of the numerical method are presented. The remarkable difference between integer, fractional, and variable-order fractional models of diffusion equation can be seen easily. The variable-order fractional model describes the characteristics of denoising with more accuracy compared to integer and fractional-order. Thus, the fundamental goal of this work to construct an image denoising algorithm1 for the variable-order fractional diffusion equation in space and time by using finite difference approximation is achieved. The present work shows the validity and great potential of variable-order fractional diffusion equation for denoising the noised images. The experimental results show that the quality of denoising images and the highest PSNR and least MSE values are obtained by a fractional variable-order diffusion model.



---

*Source: 8050017-2021-12-22.xml*