Image denoising method based on variable exponential fractional-integer-order total variation and tight frame sparse regularization

This study presents a variational image restoration algorithm based on variable exponential fractional-order total variation (TV), variable exponential integer-order TV and tight frame sparse regularization. The energy functional of this variational problem is composed of four parts: a fractional-order TV regularization term with a variable exponent, an integer-order TV regularization term with a variable exponent, a data ﬁdelity term and a tight frame regularization term. The variable exponent is a function of the gradient information of the image. The ﬁrst three parts of the energy functional is the smooth part, and the last part is the non-smooth part. To minimize the smooth part, the optimization problem can be simply transformed into a gradient descent ﬂow using the variational method. For the non-smooth part, we use the soft-thresholding method over the sparse framelet coefﬁcients. As the combination of the fractional-order derivative and integer-order derivative, the variable exponent and the sparse regularisation, the proposed method can effectively remove noise of the image, protect the image boundary, and keep image texture details well. Experiments with simulated data and real


INTRODUCTION
In practical applications, signals are always deteriorated by different types of noises in the process of acquisition, transmission and conversion, which results in image quality degradation. Therefore, image denoising is still an important problem in image processing. It is hard to develop a universal algorithm to deal with all kinds of noise. Here we consider the noise model where u 0 represents the noisy image, u represents a clean image to be recovered, and n is the noise dependent on the image u or not. The so-called image denoising is to remove the noise n as much as possible with the loss of important features of images such as edges, irregular small scale texture details as little as possible. Hence it is of great importance to develop a suitable regularization parameter) which controls the smoothness of image, and Ω is the compact support domain. The first term in the energy functional of (2) is a regularization term which ensures the regularity of the solution and the suppression of noise, and the second term is a fidelity term which reduces excessive loss of image information. However, this model has some inherent disadvantages: it is an isotropic diffusion model and the function in the W 1,2 (Ω) space has strong smooth characteristics, so high-frequency components of images like boundaries feature information cannot be well preserved. Therefore, it is necessary to consider anisotropic diffusion models.
In 1992, Rudin, Osher and Fatemi proposed the following representative total variation regularization problem, called ROF model in [2] min u∈W 1,1 Considering the image denoising problem in space W 1,1 (Ω), this model avoids the problem of excessive smoothness for high frequency information of image and effectively maintains the boundary information of the image. However, due to a priori piecewise constant assumption, it can result in a staircase effect [3][4][5].
In order to effectively reduce the staircase effect and oversmoothing problem, taking the advantages of models (2) and (3), an anisotropic diffusion model in the Sobolev space W 1,p (Ω), p ∈ (1, 2) is a natural idea min u∈W 1,p Then, Blomgren et al. made some improvement to the above model, and proposed a variable exponential functional model [6] min u∈W 1,p with p a monotonically decreasing function that satisfies lim s→0 p(s) = 2, and lim s→∞ p(s) → 1. This model can maintain the boundary features and deal with the flat region well. And due to the variable exponential p, this model can control the diffusion adaptively. However, the exponent p depends on the gradient information of the image, which increases the amount of computations and it is hard to establish the lower semi-continuity of the energy function to analyse the existence, uniqueness and the stability of the model solution.
Then, Chen, Levine and Rao proposed the following model [7] min u∈BV (Ω)∩L 2 (Ω) where , |r| > , is the Gaussian kernel, k and are two fixed positive parameters, and is a user-defined threshold. They analysed the model (6) in the bounded variation function space BV (Ω). However, there are too many parameters and the energy functional is much more complicated, which increases the calculation amount of the algorithm and reduce the practical applicability.
Inspired by the ideas above, Li et al. proposed the following simplified model (VarPTV) [8] min where . Model (7) is simpler than model (6) and has fewer parameters. When the model (7) deals with the boundary region, the gradient value is large, and then p → 1. Hence this model approximates the ROF model (3), and it can preserve the image boundary information well. For the smooth region, the gradient value is small, and then p → 2. Consequently, this model can take the advantage of the smooth denoising of the model (2). For other regions, the adaptive processing can be performed by the exponent p(x). They analysed this model in the variable exponential function space W 1,p(x) (Ω) and proved the existence and the uniqueness of the minimization problem (7). What mentioned above are all based on integer-order partial differential equations. Although integer-order partial differential equations can effectively remove noise, due to the locality of integer-order derivatives, there are still some shortcomings [9] in the preservation of small-scale texture details. In recent years, fractional-order partial differential equations have been extensively studied in different computational fields [10][11][12]. As a generalization of integer-order derivative [13][14][15], the fractional-order derivative has properties of memory and non-locality. Moreover, the function space with fractional-order derivative is larger, which makes more types of image information can be processed. Therefore, applying the fractional derivative model to the image denoising problem [16,17] can eliminate the staircase effect and preserve fine details of the image better.
In [17], Chen et al. proposed a fractional-order total variation image restoration based on primal-dual algorithm, They solved the minimization problem by primal-dual method, used Grünwald-Letnikov fractional derivative definition, and chose the maximum permissible error as the iterative stopping criterion. Different from the method in [17], we previously proposed the following model combing the fractional-order derivative with the integer-order derivative [18].
where D v denotes the fractional-order derivative operator, ,ũ 0 is the pre-processed image by a Gaussian filter. The gradient of the image is large on the boundary, so function g → 0, and then p → 1. In the smooth region, the gradient is relatively small, so the function g → 1, and then p → 2. And we solve this minimization problem by solving the associated Euler-Lagrange equation, use the fractional derivative in Fourier domain for computation simplicity, and take the maximum peak signal-to-noise ratio as the iterative stopping criterion. The model has a good suppression of noise, but there is still some noise in the smooth domain. If we want the noise to be removed further, image details will be lost much more. Hence we should make a good tradeoff between noise reduction and details preservation.
Inspired by these models and in order to get a better denoising algorithm, we proposed an improved denoised model based on Fractional-order total variation (FTV) regularization with a variable exponent, TV regularization with a variable exponent and tight frame sparse regularization. Here TV regularization with a variable exponent can effectively avoid over-smoothness problem and staircase effect and is able to preserve discontinuities and image structures. Tight frame sparse regularization can reduce the noise further. And we solve this minimization problem by solving the associated Euler-Lagrange equation, use the fractional derivative in Fourier domain for computation simplicity, and take the maximum peak signal-to-noise ratio as the iterative stopping criterion. Taking advantages of each part and the interaction of these three parts, our proposed method can remove noise well, greatly preserve the texture details and at the same time reduce computation time a lot.
The reminder of this paper is organized as follows. In Section 2, we elaborate the overall definition of the energy function and the derivation of the gradient descent flow equation. In Section 3, the numerical solution of the gradient descent flow equation is given and five test images and two real head phantom CBCT reconstructed images are tested. Experimental results validate the effectiveness of our proposed model. Conclusions are given in Section 4.

PROPOSED MODEL
In this section, we give an energy functional in the continuous function space with a tight support domain Ω, where E f is a fractional-order TV regularization term with a variable exponent p(x), E t is TV regularization term with a variable exponent p(x), E t f is tight frame (TF) regularization term, and E d is a data fidelity term. , and are three positive parameters, which determines the contribution of each term to the entire functional.

Fractional-order TV regularization term with a variable exponent p(x)
Boundary and texture details are important information of the image. The more boundaries and textures, the richer the structural information. And then the image will be much clearer. From a mathematical point of view, in order to better preserve these important details in image denoising, a related operator must be found to suppress noise without excessively smoothing boundary and texture details.
The integer-order derivative is defined in a tiny neighbourhood of a certain point, and it has locality. However, the texture of the image is global, so the integer-order derivative cannot satisfy the above requirements. A low-order integer derivative will cause a staircase effect on the denoised image, while high-order integer derivative will oversmooth the image, destroy the image boundary, and blur the image texture. Unlike the integer-order derivative, the fractional derivative is defined over the entire domain. And it has properties of memory and non-locality, which is more appropriate to deal with medium-low-frequency details such as textures.
Inspired by the ideas above, it is especially important to choose an appropriate order v between the low order and the high order, so that the advantages of the high and low orders of derivatives can be taken, and their disadvantages can be reduced. Moreover, the exponent of the norm of the fractional-order TV is also very important. Therefore, we consider a fractional-order TV model with a variable exponent p(x) (VarPFTV).
, and D v denotes the fractional-order 2 . In order to obtain the gradient descent flow equation of (11), we calculate the Euler-Lagrange equation by the variational method. Take an arbitrary test function (x, y) ∈ C ∞ (Ω), and define Taking the derivative of Φ with respect to , we have Let = 0, and we get Let D v * x and D v * y be the adjoint operators of D v x and D v y respectively. Since Φ( ) achieves its minimum at = 0, we get Thus, by the arbitrariness of , we obtain the Euler-Lagrange equation of (11) And the gradient of E f is given by . (17) Consequently, the gradient descent flow equation of (11) is formulated as .

TV regularization term with a variable exponent p(x)
Due to the high sensitivity of the fractional-order derivative, a little of noise in the smooth region is considered as details of the image in processing, so the processed image has some noise particles [18]. In order to denoise further and improve the sharpness of the image, we add an integer-order TV regularization term with a variable exponent for the model. It can make full use of the advantages of integer-order TV with a variable exponent to deal with local noise elimination and preserve boundary information. And this model will make the denoised image look more natural.
The integer-order TV regularization term with a variable exponent p(x) [2] is given by where ∇ denotes the gradient operator, ∇u = (∇ x u, ∇ y u) = (u x , u y ) the gradient of u and |∇u| = ). By the variational method, similar to the fractional-order TV term with a variable exponent, we can obtain the gradient descent flow equation of (19) as ) .

TF regularization
TF is a tight frame system that can sparsely represent piecewise smooth functions. Following our previous experiences in tight frame based image processing [19], we use framelets introduced in [20][21][22] to reduce the noise. The TF regularization term is formulated as where W is a tight frame operator that satisfies W T W = I . And then the minimization of E t f reduces to where T (⋅) is the soft-thresholding operator defined by x if |x| > , and zero otherwise. Equation (22) can be considered as a denoising procedure using framelet sparsity, and it was revealed in [20] that the solution minimizes an objective function involving a balanced sparsity of framelet coefficients.

Data fidelity
In the process of image denoising, there is always some inconsistency between the denoised image u and the noisy image u 0 .
In order to penalize this inconsistency, the data fidelity term is necessary, which can be formulated as The gradient of the energy function E d can be easily calculated as ∇E d = u − u 0 . And the gradient descent flow of (23) is given by

Our proposed formulation
From overall analysis, our proposed model mainly includes two modules: the smooth part , and the non-smooth part E 2 (u) = E t f (u). For the smooth part E 1 (u), by the linear property of the gradient operator, we can calculate its gradient by Then the gradient descent flow equation of E 1 (u) can be formulated as As for the minimization of the non-smooth part E 2 (u), we use the soft-thresholding method

Numerical algorithm
In Section 2, we have discussed how the model is obtained. Now in this section, we introduce how to solve this model numerically. For the easiness to implement, we use the frequency domain definition of the fractional-order derivative. Assume that the original image u is of size m × n. D v x u and D v y u are the fractional-order derivatives of order v for the given image u. 2 is the norm of fractional gradient operator. In practical calculations, using the characteristic of the 2D discrete Fourier transform (DFT), fractionalorder central partial difference formulation (refer to [16]) can be given by where w 1 , w 2 = 0, 1, … n − 1 are DFT frequencies that correspond to x and y, respectively. F −1 is the inverse operator of 2D DFT. The adjoint operators D v * x u and D v * y u can be calculated by the following formulation respectively.
where conj(⋅) is the complex conjugation.
The purely diagonal operator of K in the frequency domain is defined by So, we can obtain Since the adjoint of is F is F −1 , we can get where D v * y u can be calculated by the same method.
There are a family of framelet systems constructed in [22] from B-splines with different orders. In the univariate case, the piecewise constant B-spline framelet system is the well-known Haar wavelet system, whose filters are The filters of the piecewise linear B-spline framelet system are  Filters for a bivariate framelet system can be constructed by tensor products of the corresponding univariant filters. Because the image is 2-dimensional, we will use this bivariate framelet system in our experiments. In order to approximate (10), we first review the notations applied in the finite differences schemes. Let Δt be the time step, h = 1 be the space step, and k be the maximum iteration number. Then noise removal method for solving (10) can be written as Step 1: Let the input noisy image be u 0 and set n = 1, u n = u 0 , k, Δt, t = nΔt , , , , k.
Step 2: Calculate the 2D DFT of u n asû n . Calculate the v-order central partial differences D v x u n and D v y u n using (28) and (29).
Step 3: Calculate the Step 4: Calculate the d n = F −1 ⋅ K * 1 ⋅ F (l n D v y u n ) + F −1 ⋅ K * 1 ⋅ F (l n D v y u n ) according to the results of Step 2, Step 3.
Step 6: Calculate u n+1 = W T T W (u n − Δt ( ⋅ d n − ⋅ b n + (u n − u 0 )) and set n = n + 1; if n = k or satisfy the iterative stopping criterion (the maximum peak signalto-noise ratio), output u n+1 , stop; else go to Step 2.

Experimental results
In this subsection, we provide several experiments with the simulation data of five test images and two real head phantom CBCT reconstructed images. The test images "Shepp-Logan", "Barbara", "Boat", "Baboon", "Elaine" without and with additive white Gaussian noise of standard deviation = 20 are shown in Figures 1-3. Two slices real head phantom CBCT reconstructed images at normal dose 20 mAs and low dose 10 mAs are shown in Figure 4. Firstly, we did experiments with the five test images with different noise levels. Then, we compared our proposed method with VarPTV model [2], VarPFTV model (our model with = 0 and = 0) and IPM model [23]. Secondly, in order to further validate the effectiveness and increase the possibility of application of our proposed method, we also did experiments with two real images of low-dose (10 mAs) head phantom CBCT reconstructed data provided by the laboratory of Shandong Xinhua Medical Instrument Co., Ltd., and gave comparison results. In addition to subjective evaluation, the peak signal-to-noise ratio (PSNR) and the structure similarity (SSIM) are used as objective indexes to quantify the denoising effect of experimental results. PSNR and SSIM are defined as follows PSNR = 10 × log 10 where u is the original image without noise,ū is the denoised image, u and̄u are mean values of u andū respectively, u and̄u are standard deviation of u andū respectively, uū is the covariance of u andū, c 1 and c 2 are two small positive number to avoid zeros of the denominator, and m, n are the size of u. Generally the larger the value of PSNR and SSIM are, the better the performance is. For all of these algorithms, the numerical iteration will stop at the point where the PSNR value is the highest. In order to elaboratively show good performance of our method, we also give the enlarged parts of denoised images. As there are three regularization term in our method, to analyse the computation complexity, we present the running time.

Experiments with simulated data
We present experimental results at noise level -the noise standard deviation = 20. Figures 5-14 show the experimental results on the test images. In each group of experiments, the noise image of the test image, denoised image and enlarged image by different methods are presented.
In Figure 5, from the denoising image, there remains some noise in Figure 5(a) and 5(b) processed by the VarPTV method and the VarPFTV method. The edges processed by IPM method in Figure 5(c) are unclear and contain some artefacts. In Figure 5(d) of our proposed method, we can hardly see any residual noise and the edges are preserved very well, which can be further validated from the enlarged images in Figure 6(a)-6(d). It can be seen that our method can remove noise well and at the same time preserve edge details.
In Figure 7, from the denoising image, there remains some noise in Figure 7   proposed method, we can hardly see any residual noise. In addition, the flat area of the ground, face and arms are very smooth. Moreover, from the enlarged images, compared with the images in Figure 8(a)-8(c), it can be seen that our method in Figure 8(d) can remove noise well and at the same time preserve the texture details.
In Figure 9, there are more detail elements in the test image, and the gray clouds in the sky are relatively low in gray, which is extremely easy to be affected by the additional noise. Because of the high sensitivity of the fractional-order derivative, the clouds area are easy to be treated as noise in the denoising process. Consequently, the resulting image could show some gain effect. From the denoising image in Figure 9(a)-9(c), one can see that the VarPTV, VarPFTV and IPM methods are not good at dealing with the sky part, and the gray scale contrast of the characters on the hull and the figure is inferior to our proposed method. The image in Figure 9(d) processed by our method contains the least noise, the clouds in the sky are smooth, and     Figure 10(a)-10(d). These results validate the superiority of our proposed method.
In Figure 11, there are fewer smooth regions in this test image, and there are more types of details around the hair, so it is beneficial for us to verify the ability of our proposed method to preserve the details of the image. From Figure 11 From the enlarged images in Figure 12(a)-12(d), the superiority of our method is much more clear.   In Figure 13, it can be seen that the image in Figure 13(a) and 13(b) processed by the VarPTV and VarPFTV methods still have clear noise. Comparing images in Figure 13(c) and 13(d), the image contrast in Figure 13(d) is better than the image in Figure 13(c) on the background, the face and the scarf. From the enlarged images, one can see that the VarPTV model and the VarPFTV model in Figure 14(a) and 14(b) has much more noise, and the IPM model in Figure 14(c) causes additional black and white point-like artefacts and the contrast is worse than that of our method in Figure 14(d). This further demonstrates that our proposed method has better denoising ability than the VarPTV, VarPFTV and IPM models.
Tables 1-6 presents the objective evaluation -PSNR and SSIM value after image processing by our proposed method and    the other three methods. It can be seen from these tables that PSNR and SSIM values of our proposed method are the largest compared with the VarPTV model, VarPFTV model and IPM model. Moreover, with the increase of the noise level, the PSNR of the proposed method is also the largest, which shows the robustness of our algorithm to noise.
To analyse the computation complexity, we give the running time of test images with noise level = 20 denoised by different methods in Table 7. It can be seen that by the interaction of the three regularization terms our method takes the least time compared with other three methods with one single regularization term.
As we all know, the choice of regularization parameter is very important for the denoising performance of algorithms with regularization terms.   show the denoised images and enlarged images processed by the four different methods for real head phantom CBCT reconstructed data. From the denoised images, one can hardly see noise in Figure 15(c) and 15(d) and 17(c) and 17(d) compared with images in Figure 15(a) and 15(b) and 17(a) and 17(b). Compared with our method in Figure 15(d) and 17(d), the images in Figure 15(c) and 17(c) cause line-like artefact around edges and point-like artefacts in other area, which can be further validated from the enlarged images in Figure 16 and 18. That is to say, our proposed model have better smoothing effect and the sharpness of the boundary is stronger than the other three methods, which shows the effectiveness of the proposed method. Tables 8 and 9 present the objective evaluation -PSNR and SSIM values of real head phantom CBCT scan data after image processing by our proposed method and the other three methods. Here we calculate the PSNR and SSIM values with normal dose data as a clean image. It can be seen from these tables that the PSNR and SSIM of our proposed method is the largest compared with the VarPTV model, VarPFTV model and IPM model, which shows the superiority of our method. Table 10 shows the running time of different algorithms. It can be seen that our method takes the least time and gives the best denoising performance.

CONCLUSION
In this paper, combining the advantages of integer-order and fractional-order derivative, a variational denoising model with a variable exponential fractional-order TV is proposed. The energy functional of our proposed model consists of four parts: one is a fractional-order TV regularization term with a variable exponential function p(x), and p(x) is a function of the gradient information of the image, which can adaptively control the diffusion process; one is the integer-order TV regularization term with a variable exponential function p(x), which takes advantage of the integer-order derivative, and adaptively controls the diffusion; one is the data fidelity term; and the last part is TF regularization term. In the experiments, the calculation of fractionalorder derivative is performed in the frequency domain in order   to improve the speed of calculation by the central partial difference. Experimental results show that the proposed method can effectively remove noise, retain important local features, and make the denoised image look more natural. In addition, due to the combination of three regularization terms, our proposed method can reduce computation time a lot.
Here the regularization parameters are determined by experience. We think that determining the regularization parameters by deep learning method will be a good choice, which will be our future work.