A Study of Adaptive Fractional-Order Total Variational Medical Image Denoising

Abstract

Following the traditional total variational denoising model in removing medical image noise with blurred image texture details, among other problems, an adaptive medical image fractional-order total variational denoising model with an improved sparrow search algorithm is proposed in this study. This algorithm combines the characteristics of fractional-order differential operators and total variational models. The model preserves the weak texture region of the image improvement based on the unique amplitude-frequency characteristics of the fractional-order differential operator. The order of the fractional-order differential operator is adaptively determined by the improved sparrow search algorithm using both the sine search strategy and the diversity variation processing strategy, which can greatly improve the denoising ability of the fractional-order differential operator. The experimental results reveal that the model not only achieves the adaptivity of fractional-order total variable differential order, but also can effectively remove noise, preserve the texture structure of the image to the maximum extent, and improve the peak signal-to-noise ratio of the image; it also displays favorable prospects for applications in medical image denoising.

1. Introduction

In recent years, image denoising has become one of the important steps of image processing in medical diagnosis. Denoised medical images can provide reliable analysis data for clinicians who read large numbers of medical images daily [1,2]. Due to the sophistication of medical imaging, which can be easily disturbed by many factors, such as the complexity of the human body and the uneven density of the target medium, it is difficult for traditional image denoising algorithms to achieve positive denoising effects for medical images [3,4,5]. Among many algorithms for image denoising, this paper adopts the total variational algorithm of interest to improve the effectiveness of medical image denoising.

The Total Variational (TV) image denoising algorithm [6] was proposed by Rudin and Osher in 1992. It is a well-known case of applying partial differential equations to image denoising. The TV model suppresses image noise by minimizing the energy function. Yet, the edge direction of the flat area does not exist, and diffusion only along the edge direction will lead to insufficient noise suppression in the flat area. Currently, research scholars have proposed some improved algorithms for medical images based on the traditional variational method [7,8]. However, these improvements still cannot eliminate the “step effect” that results from excessive image smoothing when noise is removed from the image using the variational method. In the past few years, several advanced fractional image processing methods have demonstrated that fractional-order differentiation can enhance the texture details of images nonlinearly. Researchers such as Yan et al. [9] designed a new edge-preserving filter for image decomposition using fractional-order differential operators and $L_1$ regularization terms. Zhang et al. [10] applied fractional-order differential masks to image fusion preprocessing. The gap between low- and high-frequency signals in the source image was widened. This enabled the preservation of the image edges while effectively reducing unreasonable hole tone inside the image. Many research scholars have successfully applied fractional-order calculus theory, combined with TV theory, to image processing. This allowed researchers to effectively alleviate the problem of excessive smoothing in smooth areas of images by variational methods [11,12,13]. Zhang et al. [14] defined a new space of fractional-order bounded variational functions to partition image textures of different scales into different negative Sobolev space functions, and proposed a class of fractional-order multiscale variational models for image denoising. The model more effectively preserved the texture details with little grayscale variation in the smooth region of the image. However, it was difficult to obtain the right order to generate the best result from the denoising model. Numerous scholars have conducted numerous studies on how to determine the fractional order, but all of them have only considered the local features of the image [15,16,17]. Yu et al. [18] presented an approach to adaptively determine the order of fractional-order differentiation by calculating the local variance to reflect the local texture complexity of the image. Ullah et al. [19] proposed a trial-and-error method to adaptively adjust the fractional order so that the fractional-order total variational denoising model yields better denoising results.

Although the work of the above authors can adaptively determine the order of the fractional-order variational model for different regions in an image, these methods are used for specific metrics for a particular image and still have limitations for use in medical imaging. In recent years, swarm intelligence optimization algorithms are gradually being applied to medical image processing [20,21,22,23]. The simplicity of its implementation and its ability to provide an adjustment to the global situation have attracted the attention of an increasing number of scholars. Among them, the sparrow search algorithm is a new swarm optimization method inspired by the foraging and anti-feeding behavior of sparrow colonies [24]. The algorithm is significantly better than both the gray wolf optimization algorithm and the particle swarm optimization algorithm in terms of accuracy, convergence speed, stability, and robustness. Zhang et al. [25] proposed a semi-supervised integrated classifier based on an improved sparrow search algorithm and applied it to lung computed tomography image detection. Xiong et al. [26] proposed a fractional-order chaotic sparrow search algorithm for the enhancement of long-range red film images. Therefore, using the sparrow search algorithm for the denoising of medical images to retain their more detailed information has potential application value.

Based on previous work, an adaptive fractional-order total variational medical image denoising model based on an improved sparrow search algorithm is proposed to retain more texture details of images while denoising. The main difficulties addressed this paper include how to find the appropriate order of medical image denoising in the model and how to avoid the sparrow search algorithm from falling into the local optimum. The main contributions of this article are summarized as follows:

• For medical imaging, a new fractional-order total variational medical image denoising algorithm is proposed, which combines the amplitude-frequency characteristics of the fractional-order differential operators and the denoising advantages of the TV model, effectively mitigating the “step effect” generated by the TV model when denoising;
• An improved sparrow search algorithm that incorporates a sine search strategy and a diversity variation processing strategy is proposed to avoid the optimization algorithm from falling into local optimal solutions. Additionally, a multi-objective fusion maximization fitness function is proposed to comprehensively evaluate the denoising effect of the optimization model;
• For the feature information of medical images, an optimized fractional-order total variational medical image denoising model based on a fused multi-strategy improved sparrow search algorithm is proposed. The optimization model not only possesses robust adaptivity, but also retains more detailed texture information of medical images while denoising, which can provide more information to help clinicians in diagnosis.

The rest of the paper is structured as follows. In Section 2, the TV model and the defined form of fractional-order differentiation are briefly introduced, and a fractional-order total variational model for medical images is given. In Section 3, the main research of this paper is given. In Section 4, the results of the denoising experiment comparison are given. In Section 5, the conclusion is given.

2. Fractional-Order Total Variational Medical Image Denoising

2.1. Total Variational Denoising Model

The TV model preserves the edge information as much as possible. In general, it is common to use unconstrained extreme value models by introducing Lagrange multipliers $\lambda$, defining the minimized energy generalized form as follows:

$$\min J(u)=\frac{\lambda }{2}{\displaystyle {\int }_{\mathsf{\Omega }}{\left|u-u_0\right|}^{2}dxdy+{\displaystyle {\int }_{\mathsf{\Omega }}\left|\nabla u\right|dxdy}},$$

(1)

where $\lambda >0$ is the Lagrange weight parameter; $u$ is the image after denoising at $x,y\in \mathsf{\Omega }$; $u_0$ is the original image at $x,y\in \mathsf{\Omega }$; $J$ is the energy function of $u$; $\left|\nabla u\right|$ is the image gradient mode with $\left|\nabla u\right|=\sqrt{{\nabla }_x^{2}u+{\nabla }_y^{2}u}$; and ${\nabla }_{x}u$ and ${\nabla }_{y}u$ denote the gradients of $u$ in the $x$ and $y$ directions, respectively. The first term of Equation (1) is the residual, which guarantees that the $u$ retains the main features of the $u_0$. The second term is the regularization term, which serves to eliminate the noise in the image while maintaining the edge information as far as possible.

By the gradient descent method [27], the partial differential equation of Equation (1) is obtained:

$$\frac{\partial u}{\partial t}=\lambda (u-u_0)-\mathrm{div}(\frac{\nabla u}{\left|\nabla u\right|}),$$

(2)

where $t$ is the time step, and div is the scatter of $u$. The local coordinate expression is

$$\frac{\partial u}{\partial t}=\frac{1}{\left|\nabla u\right|}{u}_{\xi \xi },$$

(3)

where $u_{\xi \xi }$ denotes the tangential direction of the $u$. From the analysis in Equation (3), it is obtained that the diffusion operator of this model diffuses only along the image edge direction, which has the characteristic of anisotropic diffusion. When denoising flat areas of an image, the “step effect” is caused.

2.2. Fractional-Order Total Variational Medical Image Denoising Model

Fractional-order calculus is an extension of integer-order calculus, which has been in use for more than 300 years [28,29,30]. So far, there are three main expressions for the definition of fractional order differentiation: Riemann–Liouville definition, Grünwald–Letnikov definition, and Caputo definition [31]. For one-dimensional signals $f(t),t\in [a,t]$, the Riemann–Liouville definition is as follows:

$${}_a{}^R\mathcal{D}_t^{\alpha }=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\left(\frac{d}{dt}\right)}^n{\displaystyle {\int }_a^t\frac{f(\tau )}{{(t-\tau )}^{\alpha -n+1}}d\tau },$$

(4)

where $\mathcal{D}$ denotes fractional-order differentiation; $\alpha$ denotes fractional order differentiation with $n-1<\alpha <n$; and $\mathsf{\Gamma }(\cdot )$ denotes the Gamma function. The Caputo is defined as follows,

$${}_a{}^C\mathcal{D}_t^{\alpha }=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\displaystyle {\int }_a^t\frac{f^{(n)}(\tau )}{{(t-\tau )}^{\alpha -n+1}}d\tau }.$$

(5)

Both the Riemann–Liouville definition and the Caputo definition use the Cauchy integral formula, which is computationally complex. The Grünwald–Letnikov definition can be converted into convolutional form during numerical implementation, and is more adaptable to image signal processing. Therefore, the definition of Grünwald–Letnikov fractional-order differentiation is adopted in this paper. The Grünwald–Letnikov is defined as follows,

$${}_a\mathcal{D}_t^{\alpha }f(t)=\frac{d^{\alpha }f(t)}{d{t}^{\alpha }}=\underset{h\to 0}{\lim }{h}^{-\alpha }{\displaystyle \sum _{j=0}^{\left. [\frac{t-a}{h}\right]}{(-1)}^j\left(\begin{array}{c}\alpha \\ j\end{array}\right)}f(t-jh),$$

(6)

where $\left(\begin{array}{c}\alpha \\ j\end{array}\right)=\frac{\alpha (\alpha +1)\cdots (\alpha +j-1)}{j!}$ is the combination parameter.

According to Equation (6), if the period is divided by $h=1$, $n=\left. [t-a/h\right]=[t-a]$. Thus, the fractional differential expression is obtained [32]:

$$\frac{d^{\alpha }f(t)}{d{t}^{\alpha }}\approx f(t)+(-\alpha )f(t-1)+\frac{(-\alpha )(-\alpha +1)}{2}f(t-2)+\cdots +\frac{\mathsf{\Gamma }(-\alpha +1)}{n!\mathsf{\Gamma }(-\alpha +n+1)}f(t-n).$$

(7)

Only the case of $0<\alpha <2$ is considered in this paper. Under the Grünwald–Letnikov definition, the amplitude-frequency characteristic curves with fractional-order differential orders of 0.2, 0.4, 0.6, 0.8, 1, 1.2, and 1.4 are shown in Figure 1. The amplitude-frequency characteristic is the curve of the amplitude of the system frequency response as a function of frequency. As shown in Figure 1, the frequency response of the fractional-order differential operator is equivalent to a nonlinear filter. Noise with high signal strength in the image will be filtered out. Noise with weak signal strength is also suppressed by its nonlinearity. The fractional-order differential operator enhances the texture details of images by denoising.

As observed in Figure 1, there is some enhancement of the lowfrequency signal as $\alpha <0.4$ and $1<\alpha <2$. However, the enhancement is relatively low compared to the case of $0.4<\alpha <1$. Thus, different fractional orders have different effects on image denoising. By changing the order of the fractional-order according to the different features of different regions, the effective denoising of different regions can be achieved.

For medical images, the structural details and edges of the image are critical, and the loss of important details may cause the doctor to make an incorrect diagnosis. Therefore, distinguishing the details and noise of an image is crucial to medical image denoising algorithms. The complete medical image is divided into detail regions and edge regions. These regional features are represented by different signal values. Due to the fractional-order differential operator enhancing the detail region of the signal with low frequency effectively, the fractional-order differential operator is introduced in this paper and a new fractional-order total variational model is proposed. To achieve better denoising result, replace the first-order differential operator in the regular term of Equation (1) with a fractional-order differential operator, as shown in the following equation:

$${\left|u\right|}_{\mathrm{TV}}={\displaystyle {\int }_{\mathsf{\Omega }}\left|{\mathcal{D}}^{\alpha }u\right|dxdy}.$$

(8)

The new energy function of the TV model for medical image denoising is obtained by introducing the fractional-order differential operator ${\mathcal{D}}^{\alpha }u$ [33], which is the Fractional-order Total Variational (FTV) model:

$$\min J(u)=\frac{\lambda }{2}{\displaystyle {\int }_{\mathsf{\Omega }}{\left|u-u_0\right|}^{2}dxdy+{\displaystyle {\int }_{\mathsf{\Omega }}\left|{\mathcal{D}}^{\alpha }u\right|}}dxdy,$$

(9)

where $u$ is the denoised medical image; $u_0$ is the original medical image; $\alpha$ is the $X\times Y$ dimensional order matrix; ${\mathcal{D}}^{\alpha }u=({\mathcal{D}}_x^{\alpha }u,{\mathcal{D}}_y^{\alpha }u)$; $\left|{\mathcal{D}}^{\alpha }u\right|=\sqrt{{({\mathcal{D}}_x^{\alpha }u)}^2+{({\mathcal{D}}_y^{\alpha }u)}^2}$; and ${\mathcal{D}}_x^{\alpha }u$ and ${\mathcal{D}}_y^{\alpha }u$ denote the $\alpha$-order differentiation of $u$ in the $x$, $y$ direction.

The fractional-order partial differential equation is obtained by generalizing the partial differential Equation (2) as follows:

$$\frac{\partial u}{\partial t}=\lambda (u-u_0)-\mathrm{div}\left(\frac{{\mathcal{D}}^{\alpha }u}{\left|{\mathcal{D}}^{\alpha }u\right|}\right).$$

(10)

From the fractional-order partial differential Equation (10), the model degenerates to the TV model as $\alpha =1$. Therefore, to solve the model, there is need only to solve for ${\mathcal{D}}^{\alpha }u=({\mathcal{D}}_x^{\alpha }u,{\mathcal{D}}_y^{\alpha }u)$. Ref. [10] defines fractional-order partial differentiation of digital images $u(x,y)$:

$$\begin{array}{l}{\mathcal{D}}_x^{\alpha }u(x,y)=u(x,y)+(-\alpha )u(x-1,y)+\frac{(-\alpha )(-\alpha +1)}{2}u(x-2,y)\\ +\cdots +\frac{\mathsf{\Gamma }(-\alpha +1)}{n!\mathsf{\Gamma }(-\alpha +n+1)}u(x-n,y)+\cdots \end{array}$$

(11)

$$\begin{array}{l}{\mathcal{D}}_y^{\alpha }u(x,y)=u(x,y)+(-\alpha )u(x,y-1)+ \frac{(-\alpha )(-\alpha +1)}{2}u(x,y-2)\\ +\cdots + \frac{\mathsf{\Gamma }(-\alpha +1)}{n!\mathsf{\Gamma }(-\alpha +n+1)}u(x,y-n)+\cdots \end{array}$$

(12)

Equation (10) is solved using the difference method. Assuming that the image pixels can be represented as an $X\times Y$ dimensional matrix, the time step is $\Delta t$, the spatial grid size is $h$, and time and space can be discrete as

$$t=n\Delta t,i=xh,j=yh,$$

(13)

where $n=0,1,2,\cdots$, $x=1,2,\cdots ,X$, $y=1,2,\cdots ,Y$. Then, the numerical solution of the nonlinear FTV model (10) at pixel $(x,y)$ after $n+1$ iterations can be calculated by Equation (14):

$$u_{x,y}^{n+1}=u_{x,y}^n+\Delta t\left(\lambda (u_{x,y}^n-u_{x,y}^0)-{\left({\mathcal{D}}_x^{\alpha }\left(\frac{{\mathcal{D}}_x^{\alpha }u}{\left|{\mathcal{D}}^{\alpha }u\right|}\right)\right)}_{x,y}-{\left({\mathcal{D}}_y^{\alpha }\left(\frac{{\mathcal{D}}_y^{\alpha }u}{\left|{\mathcal{D}}^{\alpha }u\right|}\right)\right)}_{x,y}\right),$$

(14)

where the initial condition is $u_{x,y}^0=u_0(xh,yh)$ with $x=1,2,\cdots ,X$ and $y=1,2,\cdots ,Y.$

3. Fractional-Order Total Variational Model Optimization

In this paper, the framework of the proposed model is presented in Figure 2. An Improved Sparrow Search Algorithm (ISSA) is first designed, which integrates a sinusoidal search strategy and diversity variation processing strategy to avoid falling into local optima during the individual search. Then, the algorithm is introduced into the FTV model, and an adaptive FTV medical image denoising algorithm is proposed to address the problem that the order of the FTV denoising model cannot be determined, and to further improve the adaptivity of the FTV model and to preserve the edges and smoothed areas of the images to a greater extent.

3.1. Improve Sparrow Search Algorithm

The Sparrow Search Algorithm (SSA) was proposed by Xue and Shen in 2020 for solving optimization problems. The algorithm has a high merit-seeking ability and fast convergence [24]. In the standard SSA, a set of sparrows that may represent the optimal solution in the solvable space is first initialized, and the ratio of producers and predators is set. The fitness value is calculated and ranked according to the fitness function. Then, the producers and predators update their positions according to the predation rules. Each time the sparrow’s position is updated, the adaptation value is calculated at that time. Iteration is continued until the optimal value is found or the maximum number of iterations is reached.

The updates of producers’ and predators’ positions in the standard SSA adopt the same update formula with the problem of slow search speed. Therefore, the sinusoidal search strategy is introduced in the optimization algorithm of this paper. It enables the individuals with better adaptation in the original population to search near their original positions, reinforcing the local search ability of the algorithm; the individuals with poor adaptation in the original population can explore far from their positions, strengthening the global search capability of the algorithm. The formula for changing the weight of the sine search strategy is as follows [34]:

$$w=w\min +(w\max -w\min )\times \left(\sin \left(\left(\frac{f_i^t-f_{\mathrm{best}}^t}{f_{\mathrm{worst}}^t-f_{\mathrm{best}}^t}+1\right)\times \frac{\pi }{2}+\pi \right)+1\right)$$

(15)

where $w\min$ is the minimum value of the weight variation range; $w\max$ is the maximum value of the weight variation range; $f_i^t$ is the fitness value of the i-th sparrow in the t-th iteration of the population; $f_{\mathrm{best}}^t$ is the optimal fitness value of the population in the t-th iteration; and $f_{\mathrm{worst}}^t$ is the worst fitness value of the population in the t-th iteration.

Applying the values $w$ in the sinusoidal search algorithm strategy to the SSA, assuming that the number of sparrows is $n$ and the dimensionality of the medical image denoising model of order $\alpha$ to be optimized is $d$, the update formulae for producers and predators of the improved algorithm are obtained as follows:

$$x_{i,j}^{t+1}=\left\{\begin{array}{cc}{x}_{i,j}^t\cdot \exp (-\frac{1}{\alpha \cdot it})& R_2<ST\\ x_{i,j}^t+w\cdot Q\cdot L& R_2\ge ST\end{array}\right.$$

(16)

$$x_{i,j}^{t+1}=\left\{\begin{array}{cc}Q\cdot \exp \left(\frac{x_{\mathrm{worst}}-x_{i,j}^t}{i^2}\right)& i>n/2\\ x_P^{t+1}+w\cdot \left|x_{i,j}^t-x_P^{t+1}\right|\cdot A^{+}\cdot L& i\le n/2\end{array}\right.$$

(17)

$$x_{i,j}^{t+1}=\left\{\begin{array}{cc}{x}_{\mathrm{best}}^t+w\cdot \beta \left|x_{i,j}^t-x_{\mathrm{best}}^t\right|& f_i>f_g\\ x_{i,j}^t+w\cdot K\frac{\left|x_{i,j}^t-x_{\mathrm{best}}^t\right|}{(f_i-f_w)+\epsilon }& f_i=f_g\end{array}\right.$$

(18)

where $t$ denotes the number of iterations and $x_{i,j}^{t+1}$ denotes the individual position information of the i-th sparrow in the j-th dimension.

Equation (16) represents the location update of the explorer, where $it$ is the maximum number of iterations, $R_2$ is the warning value issued after the predator is detected, $ST$ is the preset safety threshold, $\alpha \in [0,1]$ is a random number, $Q$ is a random number obeying a normal distribution, $L$ is a $1\times d$ matrix, and each element is $1$.

Equation (17) is the location update of the predator, where $x_{\mathrm{worst}}$ denotes the current location with the lowest adaptation, $x_p$ is the location of the predator, $A$ is a $1\times d$ matrix with each element of the matrix being assigned a random value of $1$ or $-1$, and $A^{+}=A^T{(A{A}^T)}^{-1}$.

Equation (18) expresses the position update of a sparrow aware of the danger, where $x_{\mathrm{best}}^t$ is the optimal position of the population; $\beta$ is the random number of step control parameters obeying a mean of 0 and a variance of 1 in normal steps; $K\in [-1,1]$ is the stochastic number; $f_i$ denotes the current individual fitness value; $f_g$ and $f_w$ denote the current global optimal fitness value and the global worst fitness value, respectively; and $\epsilon$ is the smallest constant.

To solve the problem that the optimization algorithm is prone to fall into local optimum, the index $A$, which indicates population aggregation in biology, is introduced along with the sinusoidal search strategy. The expression of an indicator $A$ is as follows [35]:

$$A=\frac{\delta -\overline{x}}{{\overline{x}}^2},$$

(19)

where $\delta$ denotes the variance of sparrow population fitness and $\overline{x}$ denotes the mean of sparrow population fitness. The population exhibits an aggregated state as $A\gg 0$; when $A$ tends to 0, the population exhibits a random state. To avoid the initial appearance of iterations in the aggregated state, the population is treated using the Cauchy variation.

When the value $A$ is greater than the preset threshold with $t\le it/2$, the global optimal solution is mutated using Equation (20):

$$x=x_{\mathrm{best}}+x_{\mathrm{best}}\cdot \mathrm{Cauchy}(0,0.5).$$

(20)

3.2. Fitness Function

To consider both the quality of the medical image and the similarity with the original image, differentiating from the previous single-objective fitness function, a multi-objective maximization equation is designed as the fitness function of the ISSA algorithm in this paper. In the ISSA algorithm, the value of the fitness function for each sparrow is used to evaluate the denoising effect of medical images. The medical images processed by the ISSA algorithm should filter out most of the noise while retaining the texture structure of the original image. The fitness function is as follows:

$$Fit=w_{1}P+w_{2}S,$$

(21)

where $w_1$ and $w_2$ are constants with $w_1+w_2=1$ and are used to represent the weights of the relative importance of the objective function. Take $w_1=w_2=1/2$ to ensure a balanced denoising result in this paper.

$P$ in Equation (21) represents the peak signal-to-noise ratio (PSNR) [36], and a higher value of PSNR indicates better denoising of the image. The expression formula is as follows:

$$\mathrm{PSNR}=10\times {\log }_{10}\frac{d^2}{\frac{1}{X\times Y}{\displaystyle \sum _{x=1}^x{\displaystyle \sum _{y=1}^y{(u_0(x,y)-u(x,y))}^2}}},$$

(22)

where $d$ denotes the maximum value of the image grayscale; $x$ and $y$ denote the length and width of the medical image, respectively; $X\times Y$ denotes the number of pixel points in the medical image; $u_0(x,y)$ denotes the original medical image; and $u(x,y)$ denotes the denoised medical image.

$S$ denotes structural similarity (SSIM) [37], and the higher the SSIM value, the more similar the structure is. The SSIM expression is as follows:

$$\mathrm{SSIM}=L(u,u_0)\times C(u,u_0)\times H(u,u_0),$$

(23)

where $L$, $C$, and $H$ are defined as:

$$L(u,u_0)=\frac{2{\mu }_u{\mu }_{u_0}+C_1}{{\mu }_u{}^2+{\mu }_{u_0}{}^2+C_1},$$

(24)

$$C(u,u_0)=\frac{2{\sigma }_u{\sigma }_{u_0}+C_2}{{\sigma }_u{}^2+{\sigma }_{u_0}{}^2+C_2},$$

(25)

$$H(u,u_0)=\frac{{\sigma }_{u_{0}u}+C_3}{{\sigma }_u{\sigma }_{u_0}+C_3},$$

(26)

where ${\mu }_u$ and ${\mu }_{u_0}$ denote the mean values of images $u$ and $u_0$, respectively; ${\sigma }_u$ and ${\sigma }_{u_0}$ denote the standard deviation of images $u$ and $u_0$, respectively; ${\sigma }_u^2$ and ${\sigma }_{u_0}^2$ denote the variance of images $u$ and $u_0$, respectively; ${\sigma }_{u_{0}u}$ denotes the covariance of images $u$ and $u_0$; $C_1$, $C_2$, $C_3$ are constants; and $C_3=C_2/2$.

3.3. Steps of the Proposed Denoising Algorithm

A FTV medical image denoising method based on an ISSA is proposed to better filter out medical image noise. Unlike the FTV model with a constant order matrix, the adaptive FTV model proposed in this paper adaptively adjusts the order of each pixel point according to the edge region and detail region of medical images.

The first step is to input the medical image; $n$ sparrows are randomly selected and divided into producers and predators in a certain proportion. The fitness value of each sparrow is calculated and ranked according to the objective function (21). Secondly, the location information of each producer and predator is updated according to Equations (16)–(18), and the fitness value of each sparrow is updated according to Equation (21). At the end of the iteration, the order matrix $\alpha$ suitable for denoising different features in different regions of the input image is obtained. Finally, the order matrix $\alpha$ is substituted into the FTV model (14) to obtain the denoised medical image $u$.

To express the proposed algorithm more clearly, the specific steps of the proposed algorithm are described using a grayscale image of size $70\times 70$ in Figure 3a as an example. The steps of the ISSA-based Fractional-order Total Variational (ISAFTV) medical image algorithm are as follows:

Step 1: Initialize the FTV model. Input noisy image $u_0$ as shown in Figure 3a. After several tests, the ideal experimental parameters were obtained, where the maximum number of model iterations $K=50$, step size $\Delta t=0.2$, and weight coefficient $\lambda =5$.

Step 2: Initialize the ISSA algorithm. The pixel size of the noisy image is $X\times Y$, i.e., the dimension of the optimization problem is $X\times Y$ and the order matrix to be solved is $X\times Y$. The number of sparrows in the population is $n=200$; the search range is $[0,2]$; the maximum iteration is set as $it=100$; the weights in the sinusoidal search strategy are set as $w\max =3$ and $w\min =0$; the threshold $ST=0.8$ for sparrows to raise alarm; and the producers and predators are classified according to the $20\%$ of the total proportion of producers.

Step 3: Calculate the PSNR and SSIM of the noisy image. Obtain the fitness values of each sparrow according to the maximum objective function (21) and rank them.

Step 4: Update the location formulae. For producers and predators, update the location information according to Equations (16)–(18).

Step 5: Introduce the diversification variation processing operation and perturb the optimal position using Equation (20) when the aggregation of individuals reaches a specified value.

Step 6: Output ISSA algorithm results. At the end of the loop, the maximum fitness function value is obtained, which is used to find the fractional order $\alpha$ for medical image denoising; the visualization of the order matrix $\alpha$ is shown in Figure 3b.

Step 7: Output the processed image. The obtained order matrix $\alpha$ is substituted into the FTV model (14) to obtain the denoised image $u$ as in Figure 3c.

The pseudo-code of the ISAFTV algorithm is shown in Algorithm 1.

Algorithm 1. Pseudocode for ISAFTV algorithm, ISAFTV

Input: Medical image $u_0$; FTV model parameters: $K$, $\Delta t$, $\lambda$; ISSA algorithm parameters: order limit range, $ST$, $it$, $n$, number of explorers, $w\max$ and $w\min$. Output: Model order matrix $\alpha$; denoised image $u$. 1: Read denoised image information;

2: Initialize the population;

3: for $i=1$ to $it$ do

4:   Calculate the PSNR and SSIM of the denoised image according to Equations (22) and (23);

5:   Calculate the fitness value for each sparrow;

6:   Sort by fitness value;

7:   Update explorer and predator locations with Equations (16)–(18);

8:   Calculate population aggregation index $A$;

9:   if $A\gg 0$ then

10:    Perturb of the optimal position according to Equation (20);

11:   end

12:   Obtain the current position;

13:   If the new position is better than the previous one, update the optimal position;

14: end

15: Obtain the model order matrix $\alpha$; 16: for $n=1$ to $K$ do 17:   Denoising the image according to the fractional order total variational model (14);18: end

19: Output denoised image $u$.

4. Simulation Experiments and Results Analysis

In this section, Experiment 1 is designed to verify the effect of fractional order on medical image denoising by comparing different orders of the FTV model. Experiment 2 demonstrates that for medical images, the ISAFTV model proposed in the paper retains more detailed information during denoising by comparing it with various denoising models. Experiments were performed with Matlab R2018b on a Windows 10 (64-bit) desktop computer with an Intel Core i5, a 2.90 GHz processor, and 4.0 GB of RAM (Dell, Xiamen, China).

Real images of ophthalmic diseases (vitreous hemorrhage and vitreous mechanization) provided by the hospitals cooperating with the project are selected as the experimental objects. Figure 4a shows the image of vitreous hemorrhage. Due to blood vessel rupture, blood pools in the vitreous cavity, causing blood to escape and collect in the vitreous cavity. Figure 4b shows the vitreous mechanization image. The blood in the vitreous is not absorbed for a long time, and the fibrous membrane is formed by intravitreal mechanization. The experimental images are all $128\times 128$ in size. Due to the aging of the machine and equipment, the images contain a large quantity of noise, making it difficult for clinicians to form judgments.

4.1. Different Order Comparison Experiment

Different orders of FTV models bring different results for medical image denoising. Fractional-order differentiation can nonlinearly enhance the detailed regions of medical images. Compared to $\alpha <0.4$ and $1<\alpha <2$, the fractional order enhances the detail areas of the image more significantly as $0.4<\alpha <1$. In Experiment 1, the results of typical 1-order, 0.4-order, 0.8-order, 1.2-order, and 1.6-order models after denoising the vitreous hemorrhage and vitreous mechanization images are presented in Figure 5 and Figure 6.

Figure 5 and Figure 6 are similar, rendering it difficult to distinguish the denoising effect of different orders, except that the images become blurred after the 1-order model processing. Thus, a quantitative analysis of

Table 1. Evaluation index of the image after denoising of different orders $\alpha$.

Order	Vitreous Hemorrhage		Vitreous Mechanization
	PSNR	SSIM	PSNR	SSIM
0.2	30.2744	0.9806	32.7779	0.9655
0.4	31.0549	0.9904	33.5478	0.9934
0.6	31.0656	0.9841	33.6727	0.9838
0.8	31.6578	0.9978	34.9359	0.9920
1.0	25.1980	0.7818	26.6687	0.7953
1.2	31.5042	0.9967	34.6754	0.9943
1.4	31.4595	0.9935	34.5923	0.9904
1.6	31.4362	0.9941	34.1148	0.9871
1.8	31.1416	0.9976	33.7817	0.9715

is also needed. Table 1 details the PSNR and SSIM values of vitreous hemorrhage and vitreous mechanization images after processing at different orders $\alpha$ in the range of 0 to 2.

From Table 1, it is evident that the results are often not the most desirable with $\alpha =1$. Overall, the PSNR and SSIM values first increase and then decrease. The optimal denoising effect is $\alpha =0.8$. It is further shown that a suitable order helps to improve the denoising effect of the FTV model, as well as to better preserve the textile features of the image.

4.2. Different Model Comparison Experiments

In Experiment 2, the model proposed in this paper is analyzed in a series of comparisons with the TV mode [6], the 0.8-FTV model, the AFOTV model [16], the FNM model [3], and the F-DSG-NLM model [38]. The comparison results reveal that this model retains more details in medical images and provides better assistance to clinicians.

Figure 7, Figure 8, Figure 9 and Figure 10 show the comparative effects of the six denoising models on the denoising of ophthalmic images with the addition of Gaussian noise, with a standard deviation of 15. Figure 7a–g are images of vitreous hemorrhage and Figure 8a–g are magnified views of the portion of the lesion in the red box of Figure 7a, which is a detailed area of the medical image. Figure 9a–g show the vitreous mechanization images and Figure 10a–g show the portion of the lesion within the red box in Figure 9a, which belongs to the marginal area of the medical image.

Table 2. Objective evaluation criteria of different models for the denoising effect of vitreous hemorrhage images.

Standard Deviation	10		15		20
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
Noise image	28.1391	0.6574	24.5870	0.4819	22.1206	0.3670
TV model	25.6914	0.7989	26.1665	0.7995	25.9783	0.7579
FTV model	29.0345	0.8836	27.8864	0.7957	26.0211	0.6750
AFOTV model	29.7141	0.8131	27.8282	0.8051	26.6172	0.7236
F-DSG-NLM model	29.7963	0.8376	29.3433	0.8160	28.9241	0.7897
FNM model	29.8856	0.8388	29.2419	0.8163	28.8017	0.7898
Our model	30.4862	0.8857	29.7705	0.8331	28.9800	0.7995

and

Table 3. Objective evaluation criteria of different models for the denoising effect of vitreous mechanization images.

Standard Deviation	10		15		20
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
Noise image	28.1928	0.6497	24.6366	0.4825	22.1240	0.3658
TV model	25.9562	0.8034	26.4000	0.8007	26.2366	0.7586
FTV model	28.9670	0.8582	27.4736	0.7478	26.6295	0.7266
AFOTV model	29.8894	0.8290	28.0774	0.8212	26.9852	0.7449
F-DSG-NLM model	32.7222	0.8790	29.0509	0.8594	27.3167	0.7348
FNM model	32.7744	0.8796	29.0377	0.8595	27.2980	0.7349
Our model	33.8371	0.9532	30.7497	0.8614	27.7744	0.7635

show the performance indexes of the images of vitreous hemorrhage and vitreous mechanization corrupted by different degrees of noise after being processed by six denoising models.

Figure 7, Figure 8, Figure 9 and Figure 10 show that all six denoising models can remove the noise from the images to a certain extent. The contours of the lesion portion of the image became sharper after the TV model denoised the vitreous hemorrhage image, as shown in Figure 8b. The image became overall more blurred, as shown in Figure 7b. This is due to the diffusion operator of the TV model only diffusing along the orthogonal direction of the gradient. In processing the vitreous hemorrhage image, there is no edge in the flat area of the image (the hemorrhage portion of the red area in Figure 7a). The TV model over-smoothed the image, resulting in false edges in the blood contour, also known as the “step effect.” The denoising effect of vitreous mechanization images located at the edge is better than that of vitreous hematopoietic images. The PSNR values of vitreous mechanization after the TV model treatment in Table 3 are also higher than those of the vitreous hemorrhage images in Table 2.

The denoising results of the 0.8-FTV model are shown in Figure 7c and Figure 9c. Compared with the TV model, the 0.8-order differentiation preserves the detailed areas of the image nonlinearly, which alleviates the “step effect” caused by the integer differentiation in the TV model by over-smoothing the flat areas of the image. Therefore, it better preserves the detailed information of the blood accumulation in the part of the lesion in the vitreous hemorrhage image, as in Figure 8c. Table 2 and Table 3 show that the PSNR values as well as the SSIM values of the images processed by the 0.8-FTV model are higher than those of the TV model, but the order of the fractional order often requires many experiments to obtain, and the artificially selected order is only a constant, which cannot adapt to different feature points of different images, resulting in incomplete image noise removal and incomplete detail retention.

The AFOTV model achieves order adaptation by Bregman iteration according to the features of the image compared to the 0.8-FTV model. However, Figure 7d and Figure 9d show that the denoising effect of this model is not significantly different from that of the 0.8-FTV model.

Figure 7e and Figure 9e show the denoising results of the F-DSG-NLM model for two diseases, and Figure 7f and Figure 9f show the denoising results of the FNM model. Both models are improvements to the non-local denoising algorithm. It can be observed from the denoising result plots that the images become more blurred after the denoising of the ophthalmic ultrasound images for both types of images, which destroys the detailed information inside the vitreous in the images. It is proven that traditional image denoising methods do not apply to medical images.

The denoising results of the model in this paper for two disease images are shown in Figure 7g and Figure 9g. The noise in the image is largely eliminated in both the detailed area of vitreous hemorrhage and the edge area of vitreous mechanization. The detailed information in the image is also retained. This is because the adaptive search for the optimal order of fractional-order for different regions of the whole image using the ISSA algorithm improves the denoising ability of the fractional-order differential operator and overcomes the loss of image detail information, thus obtaining better visual results and helping doctors to diagnose diseases more efficiently and correctly.

In the numerical analysis of two disease images with Table 2 and Table 3, six denoising models outperform medical image denoising for medical images with low noise intensity than those with high noise intensity. The PSNR values of the proposed model for medical image denoising are better than other models regardless of the variation of noise intensity. Experiments prove that the proposed model in this paper has a positive denoising effect. Meanwhile, the SSIM values of the TV model and the FTV model for the two disease images after denoising are generally kept above 0.7, and the SSIM values of the AFOTV model, FNM model, and F-DSG-NLM model for image denoising generally remain around 0.8, while the SSIM values of this model for the medical images with a standard deviation of 10 added remain above 0.85 overall, and the SSIM values of the images with standard deviations of 15 and 20 added are always kept around 0.8, indicating that the details of the image structure after denoising of this model remain intact. The combined Table 2 and Table 3 demonstrate the effectiveness of the proposed model for ophthalmic ultrasound image denoising.

5. Conclusions

The adaptive FTV medical image denoising model better protects the sharp edges of the image, highlights the focal part of the patient, has a better balance between removing noise and protecting image features, and compensates for the shortcomings of the TV model in processing images. Meanwhile, the introduction of the ISSA algorithm avoids the artificial selection of a fractional-order differential order, improves the denoising ability of a fractional-order differential operator, and realizes the adaptiveness of a fractional-order differential operator. Finally, this paper verifies the effectiveness of this model by comparing the denoising effect of ultrasound images of two ophthalmic diseases.

In future research, we will focus on Ref. [9] to accurately distinguish the detailed regions as well as the edge regions of medical images by use of a fractional-order differential Sobel edge operator so that more detailed information can be retained in the process of the denoising of medical images to provide more help for clinicians. In addition, the additional term added to the regular term in Ref. [9] is able to attenuate the effect of the infrared background and provide us with new research ideas. Later, we will also consider the application of the proposed model with the additional term in Ref. [9] to the denoising of complex infrared images, and thus seek to improve the utility of the proposed model.