An image super-resolution reconstruction model based on fractional-order anisotropic diffusion equation

1.   Introduction

Image resolution is an important indicator of the ability to present the details of an image, which describes the number of pixels contained in the image. However, in many cases, due to the limitations of some hardware devices in the imaging system, people are often unable to obtain high-resolution images. Image super-resolution reconstruction refers to the technology of using a single or sequence of low-resolution images to reconstruct a high-resolution image^[1]. Image super-resolution reconstruction is a typical ill-posed inverse problem^[2] from a mathematical point of view, and its ill-conditioned features determine that image reconstruction is not a simple problem. Therefore, it is necessary to transform the morbid problem into a non-morbid problem in image reconstruction.

In recent years, the image processing technology based on variational method^[3,4,5] and partial differential equation^[6,7] has made great progress and development, and achieved good results in removing noise and preserving boundary information. For image denoising, a image denoising and decomposition model based on total variational regularization was proposed by Rudin, Osher and Fatemi^[8] (also known as ROF model). This algorithm is effective in preserving edge information. However, there is an obvious "staircase effect" ^[9]after image reconstruction due to the model allows the function to have discontinuous information such as jump and discontinuity in the BV space (bounded variation function space). Perona and Malik^[10] proposed a classical image denoising model based on nonlinear diffusion equation (also known as PM model), which is a milestone in the use of nonlinear partial differential equations to process and study images. Different denoising strategies are adopted in the homogeneous region and edge, not only the effect of noise removal is better, but also because it includes backward diffusion, it can also enhance the image boundary. However, the model is easy to oscillate in the diffusion process, resulting in a "staircase effect", and the details of the reconstructed image will be lost. In order to reduce the "staircase effect" in image processing, many improved methods have been proposed.

In order to solve the problem that the total variational regularization model is easy to produce "staircase effect", Zemin Ren^[11] proposed a fractional total variational regularized image super-resolution reconstruction model, which effectively suppresses the "staircase effect" and produces clear edges by introducing fractional derivatives. An image denoising model based on fractional-order anisotropic diffusion equation proposed by Jian Bai^[12], this model can not only eliminate noise, but also preserve edge information, and effectively reduce the "staircase effect". Yao^[13] proposed a fractional nonlinear diffusion model, which not only removes the multiplicative noise of the texture image, but also preserves the structural details of the texture image, but it will produce some isolated white spots when removing the noise, and it will also cause the "staircase effect" of the image because of the excessive smoothing when removing the noise. Yin^[14,15] proposed a model of image denoising and image super-resolution reconstruction based on fractional nonlinear diffusion driven by differential curvature, which effectively avoids speckle effect and gradient effect, but the convergence speed of the system is too slow. Abirami^[16] effectively applied the difference format of CN-GL and proposed a denoising model based on the spatial fractional diffusion equation, but the denoising effect is not outstanding when the image is high in noise. Making use of the advantages of ROF model and fractional total variational model, Fariba^[17] proposed a total variational denoising model based on the mixture of integer order and fractional order. Based on the characteristics of local entropy and gradient of image, Yu^[18] proposed an image denoising algorithm based on adaptive fractional calculus, which can effectively retain the edge and texture information of the image, but the "staircase effect" of the image is still obvious in the case of high noise.

Although good results have been achieved in denoising and reducing the "staircase effect" through different improved methods, the current image denoising and image reconstruction models based on nonlinear partial differential equations still have some shortcomings that cannot be ignored. First, the ill-conditioned problems caused by nonlinear partial differential equations cannot be solved efficiently, that is, non-local information, such as texture information, cannot be effectively processed. Second, it is impossible to efficiently avoid the "staircase effect" and fuzzy edge information generated in the image reconstruction process. In fact, when using anisotropic partial differential equations to smoothly eliminate noise for many times, the image will inevitably transition to a block image, causing a staircase effect of the image. Even though the algorithm based on partial differential equations at the present stage has achieved good results in suppressing the staircase effect, it will also cause the image edge and detail information to be incompletely retained. It is in the two problems of suppressing the staircase effect and effectively retaining the edge information. An effective neutralization method cannot be obtained. Third, the image reconstruction algorithm related to partial differential equations at this stage cannot effectively suppress the noise in the image in the case of high noise. At the same time, in the case of severe noise, the convergence speed of the system is too slow.

In this study, for images with complex structure and texture information, the research focus shifted to how to find a balanced algorithm that can effectively suppress the image staircase effect while retaining the edge information of the image. And to solve the problem that the traditional super-resolution reconstruction algorithm can not effectively suppress the noise in the image under the condition of blur and serious noise, and accelerate the convergence speed of the system, we try to propose a new fractional order adaptive fidelity image super-resolution reconstruction and denoising model based on fractional nonlinear diffusion equation. First, a gradient adaptive operator is proposed to balance the smooth image of the diffusion equation and the gradient fidelity term and the ability to retain edge information. Secondly, the washout filter is introduced into the model as the control term of the system to improve the stability of the system. In addition, the original fractional diffusion model is improved, and the order of the fractional order is adaptively changed based on the local variance of the image, which effectively reduces or cuts down the "staircase effect" in the image reconstruction process. Finally, numerical calculation and experimental comparative analysis of the algorithm proposed in this paper are carried out. Simulation results show that the algorithm proposed model can not only better reduce the "staircase effect" in the process of image super-resolution reconstruction and image denoising, but also has a good anti-noise performance, and can effectively retain the edge information of the image. The peak signal to noise ratio of the image reconstruction result and the visual effect of the image have be improved. Moreover, the introduction of washout filter can improve the stability of the system and accelerate the rapid convergence of the system to a stable state. This paper made the following contributions:

$\bullet$ $\quad$The gradient fidelity term is combined with the anisotropic diffusion system, and a new adaptive fidelity term and its coefficient calculation method is proposed. The coefficient of the fidelity term can be adaptively changed according to the gradient of the image. The use of gradient adaptive operators can balance the diffusion equation and the retention of the smooth image and edge information of the gradient fidelity term, and effectively suppress the image staircase effect.

$\bullet$ $\quad$By introducing the washout filter into the model as a control item of the system, the stability of the system is improved, and the rapid convergence of the system to a stable state is accelerated, which effectively solves the problem of the slow convergence speed of the anisotropic diffusion system.

$\bullet$ $\quad$Based on the original fractional anisotropic diffusion model, a new adaptive fractional differential operator is constructed. In this paper, we are committed to seeking the mapping relationship between the order of the fractional order differential operator and the local variance of the image, so that the order of the fractional order will change adaptively with the magnitude of the local variance of the image, avoiding large oscillations at the edges of the image, which can better preserve the texture information of the image.

Moreover, other similar image super-resolution reconstruction algorithms based on regularization or diffusion^[19,20,21] always have some shortcomings. Using some regularization techniques based on adaptive diffusion and L2 norm, in the flat area, the important image features can be well preserved, but in the edge area, the staircase effect cannot be better suppressed. In the research of this article, in addition to the realization that the image characteristics of the flat area and the edge area are well preserved in the image reconstruction process, and in the case of high image noise, the algorithm can still effectively reduce the noise. At the same time, the proposed algorithm has fast convergence and stability.

The research content of this paper is arranged as follows. We will briefly introduce the fractional-order anisotropic diffusion equation and the washout filter in section 2. In section 3, the new model of coupled fractional-order anisotropic diffusion equation and the proposed adaptive fractional gradient fidelity term are analyzed and how to introduce the washout filter to act as the control is presented. In section 4, the numerical calculation method of the new model is given, and the simulation experiments and experimental results are analyzed. The conclusion of this paper will be given in section 5.

2.   Related work

In this section, we will briefly introduce the theory of fractional calculus, fractional-order anisotropic diffusion equation and washout filter.

2.1.   Related theory of fractional calculus

In recent years, the application of fractional calculus in nonlinear fields (chaos, image, pattern dynamics and complex systems) has become increasingly prominent, and many achievements have been made in engineering applications. However, it is a pity that there is not a unified definition of fractional calculus, and its mathematical and physical meaning is not clear. At present, the definition of fractional calculus is mainly divided into four definitions in spatial domain: Caputo definition, Riemann-Liouville (R-L) definition, Grunwald-Letnikov (G-L) definition, the Marchaud definition^[22] and two definitions in frequency domain: the definition of Fourier transform domain and the definition of Laplace transform domain. Here we mainly introduce the definition in the Fourier transform domain^[11,23].

In terms of image processing, Fourier transform is to transform the image from the spatial domain to the frequency domain. As an extension of integer order calculus, fractional calculus and fractional calculus have something in common in Fourier transform. For a given arbitrary real function $f(x, y)\in L^2 (R)$, the definition of Fourier transform of integer order in frequency domain is as Eq (2.1).

According to the differential property of Fourier transform, the definition of fractional calculus in Fourier transform domain can be extended from integer order to fractional order. For any function $f(x, y)\in L^2 (R)$, the definition of the fractional partial derivative in the Fourier transform domain can be obtained as Eqs (2.2) and (2.3).

where $F^{-1}$ is a two-dimensional continuous Fourier transform operator, which means the inverse Fourier transform.

2.2.   Fractional-order anisotropic diffusion equation

The anisotropic diffusion equation was first applied to image denoising. The purpose of image denoising is to recover the information of the original image $u$ as much as possible from the known observed image $f$ and some noise statistical information (such as mean and variance), which is a typical inverse problem. Therefore, the problem of image denoising is ill-posed, and the regularization method is generally used to solve the problem. The more classical regularization method is the method of partial differential equation. In 1990, Perona and Malik^[10] proposed the anisotropic diffusion equation (2.4).

where $div$ is the divergence symbol, $C(s)$ is the diffusion function and $\nabla u$ represents the gradient of the image. And $C(s)$ satisfies the following basic properties: $C(s)$ is a monotonously decreasing function and $C(0) = 1, \lim\limits_{x\to +\infty}C(s) = 0$. According to the fluid diffusion system, in the uniform region, the corresponding $|\nabla u|$ is smaller, $C(s) \to 1$, the diffusion is similar to the thermal equation, which can well remove noise; in the edge region, the corresponding $|\nabla u|$ is larger, $C(s) \to 0$, the diffusion basically stops, which can well protect the boundary information. However, this algorithm will produce a "staircase effect" in the diffusion process, resulting in some false boundaries in the grayscale gradient region, which will affect the visual effect of the image. Perona and Malik give the two diffusion coefficients in their model.

where $K$ is a given parameter that controls the degree of diffusion. Through the improvement of the diffusion coefficient $C(s)$, many new models have emerged. Zhai^[24] proposed two gradient adaptive diffusion coefficients. It is pointed out that fractional order instead of integer order diffusion coefficient can be more effectively used for edge detection and noise elimination. Zhang^[25] introduced a new method for calculating the diffusion function, and the coefficient of the function is constructed by the hyperbolic tangent function. V. B. S. Prasath^[26] proposed a fuzzy diffusion coefficient, which takes into account the variability of local pixels in order to better remove edge noise and selectively smooth the image.

Considering that the fractional derivative is not the local information of the image, the texture information can be well processed, and the introduction of fractional order can avoid large oscillations at the edge. Bai Jian^[12] combined with the concept of fractional order, proposed an image denoising model based on fractional-order anisotropic diffusion, improved the mechanism of anisotropic diffusion, achieved a good image visual effect, and effectively reduced the gradient effect. The proposed model is as Eq (2.7).

where $t$ denotes time, $C(s)$ is a decreasing continuous function. Eq (2.7) can be obtained by minimizing equation (2.8). We consider defining the energy function in continuous image space.

In Eq (2.8), $\Omega$ is the image region, $f(\cdot)\ge0$ is an increasing function related to the diffusion coefficient, and the relationship between the function $f(\cdot)$ and $C(\cdot)$ is as Eq (2.9).

where $D_\alpha$ denotes the fractional derivative operator defined by $D_\alpha u = (D_{\alpha x} u, D_{\alpha y} u)$, and $|D_\alpha u| = \sqrt{D_{\alpha x}^2+D_{\alpha y}^2}$. By using the variational method of Functionals, the Euler-Lagrange equation of the energy function is obtained, and the gradient descending flow equation can be obtained as Eq (2.10).

The observed image is taken as the initial condition, and then the finite difference method is used to solve Eq (2.10). When the iterative step $t \to +\infty$, the final solution of the equation is obtained, but the time evolution can be stopped earlier to achieve the best effect between noise removal and boundary preservation.

2.3.   Washout filter

Washout filter^[27] is a kind of high-pass filter which is widely used in engineering. It has the advantage of not changing the equilibrium point of the original system. It has been used in many industrial fields. Its transfer function is as Eq (2.11).

where $p$ is the time constant. When $p > 0$, the washout filter works in steady state, and when $p < 0$, it works in unstable state. The equation of state is expressed as.

where, $\omega$ is the state variable of washout filter, $x$ is the input variable and $\upsilon$ is the output variable. When the system is stable, $\dot{\omega} = 0, \upsilon = 0$, which means that the input signal $x$ is washed out, so it is called washout filter.

In recent years, the washout filter chaos control technology has been applied in the control of chaotic systems and nonlinear systems due to its advantages of simple implementation, low control cost, certain robustness and suitability for engineering applications. Du^[28] applied this filter to the bifurcation control of chaotic systems, and delayed the bifurcation of balanced chaotic systems by adjusting the parameters of the controller. Zhou^[29] introduced washout filter control in the fractional-order chaotic synchronization system. Consider the following fractional-order response system 2.13 and drive system 2.14.

and the coupled system with washout filter controlling is constructed as Eq (2.15).

The research results show that the range of coupling parameter $k$ in fractional order chaotic synchronization system can be reduced by washout filter control, and the system trajectory can be converged quickly. Therefore, it is of a great significance to study the specific application of washout filter in various chaotic systems. Generally speaking, considering the advantages that washout filter can improve the anti-jamming ability of the system and accelerate the convergence speed of the system, the washout filter control technology is introduced into the anisotropic diffusion system as the control item of the system in order to improve the stability of the system.

3.   The proposed model

In this section, we will introduce our new model and its analysis. In order to solve the problem that the fractional anisotropic diffusion equation is easy to produce staircase effect and sometimes lose some fine texture in image denoising at present, and the traditional super-resolution reconstruction algorithm can not effectively suppress the noise in the image in the case of blur and serious noise, we propose a new fractional fidelity term model based on gradient adaptive. It is realized that the diffusion speed is slowed down by adaptively changing the fidelity coefficient, effectively removing noise and realizing image reconstruction, while alleviating the staircase effect. At the same time, in order to solve the problem of slow convergence of the anisotropic diffusion equation, the washout filter is introduced as the control term of the system, which makes the system converge to the stable state quickly and enhance the stability of the system. The following new model of image super-resolution reconstruction and image denoising are proposed.

where, ${div}\left(c\left|\nabla^{\alpha} u\right| \nabla^{\alpha} u\right)$ is the fractional diffusion term, $\left(\nabla^{\alpha} u-G_{\sigma} \otimes \nabla^{\alpha} u_{0}\right)$ is the adaptive gradient fidelity term, and $\lambda$ is the coefficient of the fidelity term. $c$ is the fractional anisotropic diffusion coefficient, $w$ is the state variable of washout filter, and $\beta \frac{\partial w}{\partial t}$ is the washout filter, which can be used as the control term of the model. $\beta$ is the weight control constant, which is used to balance the ratio between the gradient fidelity term and the control term, so as to maximize the effect of the fidelity term and the control term. Usually $\gamma$ is constant and $\gamma = m⁄n$, where $m, n$ is a positive integer and satisfies $(m, n) = 1$.

In the process of image super-resolution reconstruction, in order to obtain the true value of the HR image and test the effectiveness of the algorithm proposed in this paper, firstly, we will downsample the original HR image to obtain the corresponding LR image, and then use the proposed algorithm The HR image is recovered from the LR image. At the same time, in order to verify the anti-noise ability of the proposed algorithm, we will also add Gaussian noise to the down-sampled LR image, and perform the detection of the anti-noise ability of the proposed algorithm while performing image super-resolution reconstruction.

3.1.   Analysis of the new model

Considering that the gradient estimation of noisy image is an ill-posed inverse problem, we first use Gaussian function to regularize the preprocessed image, and then calculate the gradient of the image. At the same time, considering that the fractional derivative can better deal with the texture details of the image, combined with the fractional-order anisotropic diffusion system, the exponential adaptive fractional gradient fidelity term model is proposed and the energy functional is as Eq (3.2).

where $\Omega$ is the definition domain of the image, $\sigma$ is the variance of the Gaussian function, $\otimes$ represents the convolution operation, $D_{\alpha} = \left(D_{\alpha x}, D_{\alpha y}\right) = \left(\nabla_{x}^{\alpha}, \nabla_{y}^{\alpha}\right) = \nabla^{\alpha}, \left|D_{\alpha} u\right| = \sqrt{D_{\alpha x}^{2}+D_{\alpha y}^{2}}$, and $D^{\alpha}$ and $\nabla^{\alpha}$ are fractional differential operators. The expression $f\left(\left|D^{\alpha} u\right|\right)$ in the Eq (3.2) is the fractional-order anisotropic diffusion equation, and $\lambda \cdot\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right)^{2}$ is the adaptive fractional gradient fidelity term. The fidelity term $\lambda \cdot\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right)^{2}$ makes the gradient of the reconstructed image closer to that of the original image. In addition, the introduction of fractional derivative can better deal with the texture information, which can make the reconstructed and denoised image piecewise smooth and effectively prevent the staircase effect. $\lambda$ is the fractional gradient fidelity term coefficient, Xie^[30] proposed the coefficient of fidelity term which can be adjusted adaptively according to the characteristics of the image, in which the expression of the coefficient of fidelity term is as Eq (3.3).

It can be known from the Eq (3.3) that in the smooth region of the image, the gradient is smaller and the coefficient of the gradient fidelity term is larger, so the role of the fidelity term will be enhanced, making the denoising image gradient closer to the original image gradient. In the edge region of the image, the gradient is larger and the coefficient of the value of gradient fidelity term is smaller, which can better remove the noise. This method has a good effect on denoising and fidelity at the image boundary, but in some areas with non-edge or small gradient, the effect of gradient fidelity term will be over-enhanced, and the image will lead to mosaic and incomplete denoising due to excessive fidelity.

In order to solve the problems of over-enhancement and incomplete denoising caused by the gradient fidelity term coefficient in the region with small gradient in (3.3), we propose an adaptive adjusted coefficient as (3.4) for fractional gradient fidelity term.

The gradient fidelity coefficient $\lambda$ is a function with respect to $\varphi\left(\left\|\nabla^{\alpha} u\right\|\right)$, and the expression of $\varphi\left(\left\|\nabla^{\alpha} u\right\|\right)$ is as Eq (3.5).

where $T$ is a given parameter, tunable 0.4 in this article. The coefficient $\lambda$ can adaptively adjust by the image gradient $\nabla^{\alpha} u$. In the (3.4) and (3.5), we introduce the concept of image gradient of fractional order, which can be adaptively adjusted according to the local variance of the image to avoid large oscillations at large gradients such as image edges. $\varphi\left(\left\|\nabla^{\alpha} u\right\|\right)$ can be adjusted adaptively according to the gradient information of the image. In the region with large gradient, $\varphi\left(\left\|\nabla^{\alpha} u\right\|\right)$ is close to 2, and in the region with smale gradient, the value of $\varphi\left(\left\|\nabla^{\alpha} u\right\|\right)$ function is relatively large. In this way, the gradient fidelity term converges in the region with small gradient, which appropriately enhances the diffusion degree of the system; in the region with large gradient, the gradient fidelity term keeps the original value as far as possible, and alleviates the diffusion degree in the edge region of the image. It is realized that the diffusion speed is slowed down and the "staircase effect" is alleviated by adaptively changing the size of the index term.

According to Eq (3.2), taking into account the complexity of partial differential equations or irregular boundary conditions, it is difficult to obtain an analytical solution to solve the equation, and it has to be converted to an approximate solution that meets the requirements of the degree of approximation. In this paper, the variational method is used as an approximate method to solve the energy function. The so-called variational method is the extreme value problem of the functional.

The process of solving the energy functional is the process of finding its minimum value. It is easy to prove that Eq (3.2) has a minimum value, so that the minimum value is $u$. Generally, the variational method can be used to solve the energy functional and obtain the corresponding Euler-Lagrange equation. For any test function $\eta \in C^{+\infty}(\Omega)$, we can define the Eq (3.6).

let $\varepsilon = 0$, and satisfy $f^{\prime}(\sqrt{s}) = c(s) \cdot \sqrt{s}$, then

where $D_{\alpha x}^{*}$ and $D_{\alpha y}^{*}$ denote the conjugate operators of $D_{\alpha x}$ and $D_{\alpha y}$ respectively. Let $\Phi^{\prime}(0) = 0$, the Euler-Lagrange equation (3.8) can be obtained.

then the above Euler-Lagrange equation is solved by the gradient descending flow method, and the nonlinear differential equation can be obtained as Eq (3.9).

The above shows that the range of coupling parameter k in the fractional-order chaotic synchronous system can be reduced by the washout filter control, and the trajectory of the system converges rapidly. Considering the advantages that the washout filter can improve the anti-disturbance ability of the system and accelerate the convergence speed of the system, this paper introduces the washout filter control technique into the anisotropic diffusion system and uses it as the control term of the system to achieve the purpose of improving the stability of the system.

Next, coupling the washout filter in the above equation nonlinear equation, and the Eq (3.10) can be obtained.

At the same time, we consider that the fractional derivative can enhance the low frequency information of the image, and has a better suppression effect on the high frequency information of the image than the integer order, which avoids the large oscillation at the edge of the image, and can better retain the texture information of the image. and it also has a good denoising ability. Based on aforementioned analysis, we introduce the concept of fractional derivative into the proposed exponential adaptive gradient fidelity term, and construct an adaptive fractional differential operator, and the order of fractional order will change adaptively with the size of the local variance of the image. The expression of its adaptive fractional order $\alpha$ is as Eq (3.11).

where $\sigma_{x, y}$ is the variance of the image $u(x, y)$, $\min \left(\sigma_{x, y}\right)$ represents the minimum variance of the image, $\max \left(\sigma_{x, y}\right)$ represents the maximum variance of the image, $k_1$ and $k_2$ are constants with values of 0.5 and 0.693 respectively, and $k_2$ denotes the maximum value of $\vert \nabla I(x, y) \vert$, so that the local variance can be controlled between [0.5, 0.5+exp($k_2$)], making the effect of the local variance on the diffusion coefficient comparable to the effect of the gradient on the diffusion function.. We set N to be a matrix with $(x, y)$ as the center and a size of $n_1\times n_2$, then the mathematical formula of $\sigma_{x, y}$ is as Eq (3.12).

where, $\mu(x, y)$ is expressed as the mean value of the image $u(x, y)$, and its calculation formula is as follows.

As shown in the Eq (3.13), the order $\alpha$ can be adjusted adaptively according to the size of the local variance of the image. In the high-frequency parts such as complex texture or edges, the local variance of the image is larger, and the order of the fractional order is larger. The details and edge information of the image are highlighted by strengthening the high-frequency components of the image. In the flat area of the image, the local variance of the image is smaller, and the order of the fractional order is smaller, which effectively suppresses the noise in the flat region and enhances the sharpening effect.

3.2.   Theoretical analysis of gradient fidelity terms

In view of the fractional gradient fidelity term introduced in our proposed model, it can be proved theoretically that the constraint term $\lambda \cdot\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right)^{2}$ will not change the existence and uniqueness of the global optimal solution of the original model. In the process of proof, the fidelity term coefficient $\lambda$ can be regarded as a constant.

Lemma 1: if the functional is convex, the function has a unique global optimal solution.

Lemma 2: for $\forall{a, b} > 0$, when $a+b = 1$, there exists functions $u_1$ and $u_2$ such that the functional $E$ satisfies $E(a \times u_1+b \times u_2)\le a\times E(u_1)+b\times E(u_2)$, then the functional $E$ is a convex function.

Proposition 1: the fractional gradient fidelity term does not change the existence and uniqueness of the global optimal solution of the original model.

Proof: for fractional gradient fidelity terms, the Eq (3.14) holds.

that is

from Eq (3.15), it can be concluded that the fractional gradient fidelity term $\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right)$ is a convex function, so there is a unique global optimal solution for the gradient fidelity term according to Lemma 1. And the gradient fidelity term introduced in our proposed model will not change the existence and uniqueness of the global optimal solution of the original model. However, in the aspect of equation theory of classical PM model, the existence and uniqueness of the solution of the enhanced model have not been proved, and in the aspect of numerical realization, the corresponding variational form of PM equation is to solve the minimum value of the functional $\int_{\Omega} \ln \left(1+(|\nabla u| / K)^{2}\right)$, which is a non-convex function with linear growth, so the well-posedness of the solution can not be obtained. Therefore, it is difficult to prove the well-posedness, existence and uniqueness of the global solution of the hybrid model which combines the PM model and the gradient fidelity term.

In fact, the ill-posed of the equation leads to instability in the discrete problem, and there are two main ways to solve the appeal problem. The first improved method is to use Gaussian convolution to obtain the regularized equation $\frac{\partial u}{\partial t} = \operatorname{div}\left(c\left(\left|\nabla G_{\sigma}\right| * u\right) \nabla u\right)$. The existence and uniqueness of the solution of this equation have been proved. Another method is to use higher-order PM equations. A new class of super-resolution and denoising models based on fractional-order anisotropic diffusion is proposed in reference^[31,32]. These models can be regarded as generalizations of second-order and fourth-order anisotropic equations. Therefore, these two methods and their promotion methods will also be the focus of our next study, this paper will not do too much application and explanation.

At the same time, we can prove that the fractional gradient fidelity term introduced in our proposed model is a bounded function in a bounded variation space (BV space) and the function is integrable.

Proposition 2: in the bounded variation space (BV), the fractional gradient fidelity term $\int_{\Omega}\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right)$ is a bounded function.

Proof: for the fractional gradient fidelity term $\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|$, the Eq (3.16) holds.

by integrating both sides of the Eq (3.16), the Eq (3.17) can be obtained.

in bounded variation space, there is $\int_{\Omega}\left\|\nabla^{\alpha} u\right\| d x d y < {+\infty}$, $\int_{\Omega}\left\|\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\| d x d y < {\infty}$according to the property of bounded variation space (BV space). Therefore, according to the properties of bounded function, it can be concluded that the fractional gradient fidelity term function $\int_{\Omega}\left(\left\|\nabla^{\alpha} u-\nabla^{\alpha}\left(G_{\sigma} \otimes u_{0}\right)\right\|\right) d x d y$ is also bounded function. The bounded function must be integrable, so it is integrable in BV space, and the introduced fractional step fidelity term will not change the denoising and edge preserving ability of the original model.

4.   Numerical experiments

In this section, we will construct the numerical calculation scheme for the proposed image processing model, and illustrate numerical simulation experiments and analysis of experimental results.

4.1.   Numerical calculations

We adopt the finite difference method to calculate the numerical value of Eq (3.10). In the case of a two-dimensional image, $u_{i, j}$ is defined to represent the pixel value of the image at $(ih, jh)$, and $h$ is the spatial step size, usually $h = 1$. $\Delta t$ is the time step, this article sets its value to 0.01. Let $u_{i, j}^n$ and $w_{i, j}^n$ be the approximate values of $u(n\Delta t, ih, jh)$ and $w(n\Delta t, ih, jh)$ respectively. The size of the image is $Ih\times Jh$, then the partial derivatives of space variables $\partial u \verb|/|\partial x$ and $\partial u \verb|/|\partial y$ can be approximated by the central difference and expressed as Eq (4.1).

The partial derivatives of the time variable $\partial u \verb|/|\partial t$ and $\partial w \verb|/|\partial t$ can be obtained by forward difference.

Then the difference calculation format of the new model 3.10 can be obtained as 4.3.

According to the above calculation format, $u_{i, j}^0$ and $w_{i, j}^0$ are used to calculate $u_{i, j}^1$ and $w_{i, j}^1$, and then continue this iteration to obtain a series of evolutionary processes until the iteration is terminated at the best result (PSNR maximum). In Eq (4.3), $D_{\alpha x} u$ and $D_{\alpha y} u$ u are expressed as fractional partial derivatives of central difference, and $D_{\alpha x}^*$ and $D_{\alpha y}^*$ are conjugate operators of respectively. In this paper, the discrete form of fractional differential in Fourier domain is adopted, and the calculation formulas are as Eqs (4.4) and (4.5).

The conjugates corresponding to Eqs (4.4) and (4.5) are defined as Eqs (4.6) and (4.7).

Based on the appeal analysis, it can be concluded that the flow of the proposed super-resolution reconstruction algorithm is shown in algorithm 1.

4.2.   Simulation experiment and analysis

In order to verify the image super-resolution reconstruction effect and anti-noise performance of our proposed model. In the experiment, we would like to compare it to some methods which include image super-resolution by adaptive sparse domain selection^[33](ASDS), adjusted anchored neighborhood regression for fast super-resolution^[34] (A+), anchored neighborhood regression for fast example-based super-resolution^[35] (ANR), image super-resolution using deep convolutional networks^[36] (SRCNN), multimodal image super-resolution via joint sparse representations induced by coupled dictionaries^[37] (CDLSR), super-resolution via image-adapted denoising CNNs: incorporating external and internal learning^[38] (IDBP), robust single-image super-resolution via CNNs and TV-TV minimization^[39] (CNNS-TV). The experimental procedures are implemented with MATLAB R2018b. The experimental environment is: the operating system is Windows 10, the processor is i5, and the memory size is 8GB.

In the experiment, we select 13 kind of images, such as "lena", "butterfly", "head", "parrot", and "bird", to be test images. Before the experiment, the degraded low-resolution image is obtained after the original image is sampled at equal intervals at 2:1 and 3:1 as input, and then the image super-resolution reconstruction was carried out using different methods to restore the estimated high-resolution image, and finally compare the original image with the estimated high-resolution image obtained.We take peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) as the main evaluation indexes for image super-resolution reconstruction and image noise immunity quality. Among them, PSNR is used to evaluate the image quality, and a larger value of PSNR indicates a better quality of the reconstructed image; SSIM focuses on evaluating the integrity of the structural and detail information of the image, and a larger value of SSIM indicates a smaller difference between the reconstructed image and the original image. Assuming that the size of the image is $m\times n$, the original input image is $u_0 (x, y)$, and the reconstructed output image is $u(x, y)$, the calculation formulas of PSNR and SSIM are as Eqs (4.8) and (4.9).

where $\mu_{x}$, $\mu_{y}$, $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{x y}$ are the local mean values, standard deviations and cross-covariance along the $x$ and $y$ directions of the image, respectively. Generally, the images with better quality have higher PSNR and SSIM values.

4.2.1.   Image super-resolution reconstruction

First of all, the super-resolution reconstruction experiments are carried out to enlarge the image by two times and three times. After double and triple sampling, the degraded low-resolution image is obtained as input, and then the image is reconstructed by different methods, and the results are compared. In the experiment, the parameters are set to $\beta = 0.1$, $\gamma = 2/5$, $T = 0.4$, and the iteration step size $\Delta t = 0.01$. In the target data set of the image, we choose the image "butterfly" and the image "bird" to carry out the comparison experiment of the image magnification by two and three times of various algorithms respectively. The "butterfly" image and "bird" image have complex structures, contain a lot of structural information, and rich texture information, which can be used as test images for image comparison. The comparison experiment results of magnifying the image twice and three times without super-resolution reconstruction algorithm are as Figures 1 and 2.

As shown in Figures 1 and 2, where Figure 1(a) is the original test image, and Figures 1(b)-(j) are respectively the images after image reconstruction using 8 different methods. At the same time, we will enlarge the wing edge of the reconstructed image "butterfly" and the lower beak of the bird in the reconstructed image "bird" locally. From Figure 1(b) and Figure 2(d), we can see that the sawtooth effect in the image reconstructed by ASDS method is more obvious, and the ANR algorithm will produce obvious ringing phenomenon. In Figure 1(g) and Figure 2(h), although there is no obvious ringing phenomenon in IDBP algorithm and CNNS-TV algorithm, it can be seen from the enlarged image of local details that both methods have aliasing effect and artifacts appear near the edge of the wing. For Figure 1(c) and Figure 2(c), we can see that the reconstruction method based on A+ algorithm makes the reconstructed image partially blurred. It can be seen from Figure 1(e) and Figure 2(e) that the reconstruction effect of SRCNN algorithm is improved, but the whole reconstructed image is too smooth and some details are lost. Through observation, the reconstruction effect of the algorithm based on CDLSR is similar to that of the proposed algorithm, the sawtooth effect and ringing effect are less in the algorithm results, and the reconstructed image has no obvious staircase effect, so it can restore better edge information. It can be seen from the overall enlarged image and the partial detailed enlarged image that the reconstruction error of the contour of the wing edge in the "butterfly" image and the edge of the mouth in the "bird" image of the proposed algorithm is smaller than that of the CDLSR algorithm, and is closer the original image, and the details of the image are reconstructed more clearly and more realistically. To sum up, the proposed reconstruction algorithm has a better effect in the preservation of image structure texture information and the restoration of edge information. And the reconstructed image will not appear obvious sawtooth effect and ringing phenomenon, which the staircase effect is suppressed effectively, and the overall performance of reconstruction is better than that of other methods.

Because some algorithms have little difference in image details, it is not easy to observe the reconstruction effect of each algorithm from a subjective point of view, and there is a certain subjectivity. Therefore, we give the comparison of the objective evaluation indicators of PSNR value and SSIM value of each algorithm for multiple image reconstruction, and the results are shown in the Tables 2 and 3.

Tables 2 and 3 show the image magnification performance of different reconstruction methods when the target image data set is magnified twice and magnified three times, respectively. Among them, the bold value in the table is expressed as the optimal value of the relevant evaluation index under the same conditions. From the data presented in the table, under different magnification, the objective evaluation index of our algorithm has achieved most of the lead, and the proposed algorithm has the highest average objective evaluation value, and has strong competitiveness compared with other algorithms. furthermore, it shows the effectiveness of the proposed algorithm.

In order to more intuitively reflect the superiority of image reconstruction of the new algorithm proposed in this article, we draw a line chart between the PSNR value and SSIM value obtained by each comparison algorithm and the proposed algorithm and the test image. As shown in Figures 3 and 4, the abscissa represents different test images, and the ordinate represents the PSNR and SSIM values of each image. Each curve represents a different algorithm, and the red curve represents the algorithm in this article. Through the line graph, it can be more intuitively shown that the proposed new algorithm has a better image reconstruction effect.

In order to facilitate further observation and analysis, we obtain the average of the PSNR and SSIM values of each algorithm for all test images. At the same time, the average PSNR value of each algorithm is used as the X axis, and the average SSIM value is used as the Y axis, and then the algorithm of this paper is used in turn Compared with the other 7 algorithms, a scatter plot is drawn according to the corresponding coordinates. As shown in Figure 5, the farther the scattered points are from the origin of the coordinates, the better the overall performance of the algorithm. From the comparison in Figure 5, it can be seen that the comprehensive performance of the A+ algorithm, the ANR algorithm and the ASDS algorithm is poor, and the comprehensive performance of the algorithm in this paper is better than other algorithms.

4.2.2.   Image anti-noise detection

In order to verify the robustness of the anti-noise ability of the proposed algorithm, in this section, we will illustrate the anti-noise simulation experimental results on "lena" and "peppers" images. First, a Gaussian white noise with a mean of 0 and a variance of 15 is added to the low-resolution image, and then the noisy low-resolution image is magnified two and three times. As shown in the Figures 6 and 7, where Figure 6, 7(a) is the original image, Figure 6, 7(b) is the image after adding noise with a variance of 15, and the rest are represented as the results of denoising and magnification by other reconstruction algorithms. From the reconstruction results, we can see that there is a lot of noise in the reconstruction results of other algorithms, except for the denoising effect of CNNS-TV algorithm, ANR algorithm and this algorithm. In all the experimental results, although there is some noise in the output of the proposed algorithm, the denoising effect is obviously better than that of other algorithms, so in the case of adding noise and image magnification, the denoising ability of the proposed algorithm is the best, and has a good anti-noise ability.

Finally, we evaluate the difference between the original image and the denoised reconstructed image objectively by calculating the two objective evaluation indexes of the noisy and enlarged image: PSNR and SSIM. As can be seen from Table 4, under the condition of adding Gaussian white noise with a mean value of 0 and a variance of 15, the PSNR and SSIM values of this algorithm are much larger than those of other algorithms, regardless of the magnification of two times or three times. From an objective point of view, it shows that this algorithm is superior to other algorithms in denoising.

In order to verify the superiority of the proposed algorithm in terms of algorithm running time when performing image denoising, we first added Gaussian noise with variance of 15 to the test image, and then performed image denoising experiments with each comparison algorithm and the proposed algorithm, and recorded the CPU running time of each algorithm, and obtained the following table information.

As can be seen from the Table 5, the proposed algorithm in this paper achieves most of the lead in running time compared with other compared algorithms, only longer than the running time of IDBP algorithm and CNN-TV algorithm. However, the IDBP algorithm and the CNN-TV algorithm have poor denoising effect, and both of them are deep neural network-based image super-resolution reconstruction models that enjoy the advantages of deep learning, but it takes a lot of time to train the data from the training database.

Meanwhile, to better visualize the adaptivity of the proposed fidelity term coefficients in this paper, we plot the waveforms of the adaptive coefficients $\lambda$ as the number of iterations $k$ increases with the number of iterations $k$, as shown in the Figure 8, using the number of iterations of the system as the horizontal coordinate and the adaptive fidelity term coefficients as the vertical coordinate. Here, we choose four test plots to verify and show that the proposed fidelity term coefficients will achieve adaptivity with the gradient of the image.

At the same time, in order to more intuitively reflect the superiority of this algorithm for image denoising while achieving image super-resolution reconstruction, we have drawn the PSNR of different algorithms for denoising "lena" and "peppers" images. And a histogram of SSIM values. As shown in Figures 9 and 10, the abscissa represents different algorithms, the ordinate represents the PSNR and SSIM values after denoising, each bar represents a different algorithm, and the red represents the algorithm in this article. It can be clearly seen from the histogram that the image reconstruction effect of this algorithm is better than other algorithms in the case of high noise.

In addition, this paper proposes to introduce the washout filter control into the anisotropic diffusion system as a control term. In order to verify whether the washout control term can improve the stability and accelerate the convergence of the system, we use the image "lena" magnification three times as a reference. The waveform of PSNR value and iteration number $k$ after image super-resolution reconstruction using washout filter as control term $\beta = 0.1$, $\gamma = 2⁄5$ and the waveform diagram of PSNR value and iteration number $k$ without washout filter $\beta = 0$, $\gamma = 0$ are shown in the Figure 11. As can be seen from Figure 11(a), after introducing washout filter control, the system tends to be stable when the number of iterations is about $k = 90$. However, when the washout filter is not introduced, the system does not tend to be stable at $k = 500$. To sum up, this paper introduces the washout filter into the anisotropic diffusion system, which can improve the stability and accelerate the convergence speed of the system. The practicability of introducing washout filter as a control term into anisotropic diffusion system is further demonstrated, which lays a foundation for the subsequent introduction of washout filter to other systems.

5.   Conclutions

In our model, based on the fractional anisotropic diffusion equation, a new method for calculating the coefficients of adaptive fidelity terms is proposed. And the proposed method is based on the relationship between image gradient and diffusion function. It is proved mathematically that the proposed fractional fidelity term will not change the existence and uniqueness of the global optimal solution of the original model. Simulation results show that the addition of the adaptive fidelity term can effectively suppress the staircase effect and aliasing effect, better reconstruct the edge and texture detail information of the image, and make the reconstructed image more clearly and naturally. And in the case of serious noise, the reconstruction algorithm effectively suppresses the noise and has a strong denoising ability. Besides, by introducing the washout filter to act as the control term of the fractional anisotropic diffusion system, the stability of the system is improved and the convergence speed of the system is accelerated.

Acknowledgements

This work was supported by NSFC Grant Nos. 61701060 and 61801067, Guangxi Colleges and Universities Key Laboratory of Intelligent Processing of Computer Images and Graphics Project No. GIIP1806, and the Science and Technology Research Project of Higher Education of Hebei Province (Grant No. QN2019069), and Chongqing Key Lab of Computer Network and Communication Technology (CY-CNCL-2017-02).

Conflict of interest

All authors declare no conflict of interest in this paper.