Hybrid High-Order and Fractional-Order Total Variation with Nonlocal Regularization for Compressive Sensing Image Reconstruction

: Total variation often yields staircase artifacts in the smooth region of the image reconstruction. This paper proposes a hybrid high-order and fractional-order total variation with nonlocal regularization algorithm. The nonlocal means regularization is introduced to describe image structural prior information. By selecting appropriate weights in the fractional-order and high-order total variation coefﬁcients, the proposed algorithm makes the fractional-order and the high-order total variation complement each other on image reconstruction. It can solve the problem of non-smooth in smooth areas when fractional-order total variation can enhance image edges and textures. In addition, it also addresses high-order total variation alleviates the staircase artifact produced by traditional total variation, still smooth the details of the image and the effect is not ideal. Meanwhile, the proposed algorithm suppresses painting-like effects caused by nonlocal means regularization. The Lagrange multiplier method and the alternating direction multipliers method are used to solve the regularization problem. By comparing with several state-of-the-art reconstruction algorithms, the proposed algorithm is more efﬁcient. It does not only yield higher peak-signal-to-noise ratio (PSNR) and structural similarity (SSIM) but also retain abundant details and textures efﬁciently. When the measurement rate is 0.1, the gains of PSNR and SSIM are up to 1.896 dB and 0.048 dB respectively compared with total variation with nonlocal regularization (TV-NLR).


Introduction
Compressive sensing (CS) [1,2] has been successfully applied in signal acquiring, processing, and compression. By exploiting the redundancy that existed in a signal, CS conducts sampling, and compression at the same time. CS theory demonstrates that a signal can be reconstructed with high probability when it exhibits sparse representation in some transformation domain, which has been widely used in many fields, such as single-pixel imaging [3], remote sensing imaging [4], and medical imaging [5].
Based on the CS theory, if an original signal x ∈ R N is sparse or after sparse transformation Ψ. The measured value y = Φx, y ∈ R M (M N) can be obtained through the measurement matrix Φ. CS theory shows that when the sparse transformation matrix Ψ and the measurement matrix Φ satisfy the restricted isometry property (RIP) [6], the original signal x can be reconstructed by solving the following optimization problem where signal α = Ψ T x, α is the sparse coefficient after sparse transformation. · 0 represents the l 0 norm. Since the solution of Equation (1) is an NP-hard problem, it can be solved by the approximate form l 1 norm min x Ψ T x 1 s.t. y = Φx.
The above-constrained problem can be constructed by the Lagrangian multiplier method to an unconstrained problem where y − Φx 2 2 is a cost function, λ is the Lagrangian parameter. Since prior characteristics exist in natural images generated by natural light imaging, the optimization problem of image compression sensing can be expressed as where u is a 2D image. R(u) represents a regularization term representing prior information of an image. The current CS recovery algorithms explore the prior knowledge that a natural image is sparse in some domains, such as discrete cosine transform (DCT) [7], wavelets [8], and gradient-domain utilized by total variation (TV) model [9,10]. Image sparse representation is a key factor to affect image reconstruction quality. Despite high effectiveness in image CS recovery, the TV-based [11] algorithms cannot recover the fine details and textures, reconstructed images suffer from undesirable staircase artifact. This problem has sparked numerous studies in designing regularizers that could efficiently suppress the staircase artifacts while retaining sharp edges. To overcome this drawback, a weighted total variation [12] is proposed to improve the sparsity of the TV norm. Furthermore, an improved total variation by intraprediction is proposed [13]. Besides, total generalized variation (TGV) [14,15] is widely discussed, which is more precise in describing pixels variations in smooth regions, and thus reduces oil painting artifacts while still being able to preserve sharp edges. Except for these TV and its variants, there also exists some combination of TV (or TGV) and different waves, such as TV + wavelet [16] and TGV + shearlet [17]. However, the methods described above are improved on TV, there is still a staircase artifact effect in the results. Further, nonlocal means (NLM) [18] uses the similarity between the surrounding pixels of the image to perform weighted filtering, which effectively utilizes the non-local self-similarity. NLM is successfully applied in the CS image reconstruction, but it often results in the painting-like effect.
Recently, a class of fractional-order TV regularization models has received considerable interest and it is widely used in image denoising [19,20]. Adaptive weighted high-frequency iterative fractional-order TV is proposed, the high frequency gradient of an image is reweighted in iterations adaptively when using fractional-order TV [21]. Another conventional way to suppress the staircase artifact is to use a high-order TV regularization [22,23], there exists a high-order total variation minimization model that removes undesired artifacts for restoring blurry and noisy images [24]. High-order TV regularization can reconstruct piecewise linear regions, but high-order TV may also smooth out the image details, and it may reduce the ability of edges-preserving [23]. In order to make full use of image prior information as much as possible to construct a regularized model. A generalized hybrid non-convex variational regularization model [25] and unidirectional hybrid total variation with nonconvex low-rank regularization [26] have been proposed.
Motivated by the aforementioned studies, both the fractional-order TV and highorder TV can solve the problem of staircase artifacts, but fractional-order TV will cause non-smooth effects just like noise in smooth areas of the image. The ability of high-order TV to reduce the staircase effect is not very ideal since it may smooth the image details. In this paper, a two-dimensional compressive sensing image reconstruction model based on hybrid high-order and fractional-order TV (HoFrTV) with nonlocal regularization is proposed. The proposed algorithm makes the fractional-order and the high-order TV complement each other on image reconstruction, which can effectively solve the problems caused by high-order and fractional-order TV. The proposed algorithm effectively reduces the problem of staircase artifacts while preserving the edges, and suppresses painting-like effects produced by nonlocal means regularization.
To effectively solve the proposed algorithm model, we introduce auxiliary variables to a construct constrained optimization problem in Section 2.2. We divide the proposed model into four subproblems to solve. These subproblems are fractional-order TV model, high-order TV model, nonlocal means regularization model, and the reconstructed image by iterative update respectively. The augmented Lagrangian method (ALM) and the alternating direction method of multipliers (ADMM) are incorporated to solve these problems. Each parameter of the Lagrangian function is optimized for the experiment results empirically. The experimental results show that the proposed algorithm outperforms the current state-of-the-art algorithms, the edges and details of the reconstructed image are more abundant, and the visual effect of the reconstructed image is better.
The remainder of this paper is organized as follows. Section 2 introduces the related regularization model and the proposed HoFrTV model. In Section 3, parameter selection and experimental results are presented. In Section 4, the conclusions are drawn. Fractional total variation can be regarded as the generalization of total variation composed of fractional orders. There are three widely used definitions of the model, including Riemann-Liouville (R-L), Grünwald-Letnikov (G-L), and Caputo model [27]. Here, we choose the G-L model because it is easier to implement for image reconstruction. The G-L model can be defined as

The Proposed Algorithm
where α fractional order, t and a are the upper and lower boundaries of the independent variable respectively, and h is the differential step-size, α k is the binomial coefficient is the Gamma function. Without loss of generality, let u denote the image. The fractional-order total variation in the horizontal and vertical directions can be expressed as and where i and j denote the pixel in the i-th row and the j-th column of an N × N size image. K ≥ 3 is the number of items involved in the computation of the fractional-order derivative, which is usually set as K = N. Based on this definition, the fractional-order semi-norm is defined as On the basis of the first-order total variation operators of the image, we can obtain the high-order total variation, defined as where (D vv u), (D vh u), (D hv u), and (D hh u) are high-order total variation operators in the horizontal and vertical directions, respectively. The high-order total variation can be seemed as performing the differential operation again on the first-order total variation.

Nonlocal Mean Regularization Model
The pixels in the image can be estimated by the weighted average of the surrounding pixels of similar neighborhoods in the similar neighborhood of the search window, which can be expressed as followsû where Ω is a search window area D s × D s . let u i and u j denote the central pixel of similar neighborhood window d s × d s respectively. The weight of u j to u i , denoted by w i,j , which is determined by the similarity of the pixels in a similar neighborhood window. It can be calculated by the Gaussian l 2 distance between pixels in the neighborhood window. It can be written in the following form Specifically, assume that u j lies in the search window of u i . The neighborhood window centered on the central pixel u j can slide in the search window to calculate the similarity between two neighboring windows. h is the factor controlling the Gaussian kernel and c is the normalization constant. Applying nonlocal regularization in the CS reconstruction process, furthermore, it can be rewritten in a matrix form as here, W represents a matrix composed of w i,j in Equation (10), let R(u) denote nonlocal regularization term. When calculating the second norm, vectorizing the elements in the norm is employed and then the second norm of the vectorized matrix is calculated.

The Proposed Algorithm Model
Incorporating Equations (8), (9), and (12) into the CS optimization problem jointly, the proposed hybrid high-order and fractional-order total variation with nonlocal regularization model image CS recovery algorithm are formulated as arg min Using variable splitting, we introduce auxiliary variables and change Equation (13) to a constrained optimization problem arg min where τ and ε control the weights of fractional-order and high-order total variation. Note that the problem of Equation (14) is quite difficult to solve directly due to the non- differentiability and non-linearity of the combined regularization terms. An augmented Lagrangian based approach is developed to solve the problem where µ 1 , µ 2 , µ 3 , µ 4 , and β are regularization parameters associated with the quadratic penalty terms. γ 1 , γ 2 , γ 3 , and γ 4 are the Lagrange multipliers associated with the constraints of Equation (15). The idea of using ADMM is to find a saddle point of L A (w, v, z, u) which is the solution to the original problem in Equation (13). It can minimize the augmented Lagrangian function L A (w, v, z, u) alternately, the problem in Equation (15) is decomposed into the four subproblems. We investigated the subproblems one by one.

w and v Subproblems
Given v, z, and u, the optimization problem associated with w can be expressed as According to Lemma 2 [28], the closed solution form of Equation (16) is When solving v subproblem, given w, v and z, similarly, the v subproblem becomes Same as w subproblems, the closed solution form of (18) can write as

u Subproblems
Fixed w, v, and z, u subproblems is equivalently expressed as We can see that the problem in Equation (20) is a quadratic function optimization problem. To reduce the calculation, the gradient descent method was used η is the step size. u can be reconstructed by every iterative update, d indicates its gradient, where η = abs(d T d/d T Gd) is the optimal step and G = µ 1 (D α ) T D α + µ 4 I + µ 3 Φ T Φ.
According to [18], the Equation (23) can be further transformed into where r = ( u − γ 4 µ 4 ), r can be regarded as an approximation of z. Since the weight matrix W represents the nonlocal means operator, the Equation (24) can be rewritten as Setting the gradient of the Equation (25) to zero, we acquire the closed-form solution as follows Finally, the Lagrange multipliers are updated by the following Once we obtained an efficient solution for each separated subproblem, the overall algorithm will be more efficient to get better reconstruction. Given all the derivations above, the specific implementation process of the proposed algorithm is described in Algorithm 1.

Experimental Results and Discussion
In this section, the experimental results are presented to demonstrate the performance of the proposed algorithm HoFrTv model. In our implementation, we chose the USC-SIPI image database, which is a collection of digitized images. The database was maintained primarily to support research in image processing, image analysis, and machine vision. Ten images were selected for verification, the original images are shown in Figure 1.

Experimental Results and Discussion
In this section, the experimental results are presented to demonstrate the performance of the proposed algorithm HoFrTv model. In our implementation, we chose the USC-SIPI image database, which is a collection of digitized images. The database was maintained primarily to support research in image processing, image analysis, and machine vision. Ten images were selected for verification, the original images are shown in Figure 1.  Reconstructed image quality was measured using the peak signal-to-noise ratio (PSNR) and the structure similarity (SSIM) [29]. The calculation expression is as follows The root mean square error is the arithmetic square root of the mean square error, it can measure the deviation between the reconstructed image and the original image. Where u ij and u ij denote the pixel values of the reconstructed image and the original image, R is the maximum value of the image gray level range. SSIM is defined as where µ x and µ y are the mean value of the reconstructed image and the original image. σ x and σ y are the standard deviations of the reconstructed image and the original image. C 1 and C 2 are small positive numbers to avoid µ 2 x + µ 2 y and σ 2 x + σ 2 y being zero, and The stopping criterion for all of the algorithms tested was set to where u t+1 and u t are the restored images at the current iterate and previous iterate respectively.

Parameter Selection
Since τ and ε were related to the weights of fractional-order and high-order total variation, the respective ranges of τ and ε were greater than 0 and less than 1. We constrained the summation of these two weights to be 1 in our proposed algorithm model. In the following content, we analyzed the influence of high-order and fractional-order weights parameters on the reconstructed results.

The Influence of High-Order
This subsection discusses the influence of high-order TV when Equation (14) does not exist in fractional-order, which means that the fractional-order is the traditional total variation. In order to investigate its effect, experiments were implemented with different weight parameters ε to high-order TV, the weights parameters value range was greater than 0 and less than 1. We randomly selected Lena image, the results are shown in Table 1. From Table 1, with the parameter ε increasing, the higher the parameter ε was, the better the result of reconstruction. This is due to high-order TV that can effectively reduce the stair casing artifacts in the reconstructed image. According to the experimental results, when parameter ε was between 0.6 and 0.9, the reconstruction results could achieve stable results. So, the parameter ε value range could be selected from 0.6 to 0.9. According to the experiment in the database, we could make appropriate selection values according to different images. In our experiments, we set ε = 0.7 and τ = 0.3 to balance the results.

The Influence of Fractional-Order
In this experiment, α is the fractional-order. The effect of this parameter on the image reconstruction performance was tested without high-order TV in Equation (14), ranging from 0.5 to 2. Image Lena was tested in the experiment. The results are shown in Table 2. From Table 2, we could see that when α < 1, the reconstructed results were poor. The reason is fractional-order TV loses more details and textures. When α = 1, fractionalorder TV converts to traditional TV. When α > 1, the larger the parameter α is, the better the textures and image details. We can see that the PSNR at α = 1 is lower than the PSNR at α > 1. We can select the fractional-order α between α = 1.3 and α = 2. However, when a value of α is too close to 2, the frequency of textures would be enhanced excessively that becomes a kind of noise. Finally, to achieve a good trade-off, in our experiments, we set α = 1.7 to balance the results. Non-local means regularization indicates the estimated value of the current pixel obtained by the weighted average of the pixels in the image that have a similar neighborhood structure. There is no doubt that the radius of the neighborhood kernel window and the radius of the neighborhood search window were vital to the experiment. If their values are too small, the self-similarity of the image cannot be fully utilized, and the characteristics of the non-local means cannot be used fully. If their values are too large, the image search area will become larger and the time will be longer, which leads to the fact that the algorithm efficiency will be lower. There are different kernel windows to search window ratio (k:s) result values for different measurement rates in Table 3. Figure 2 shows the PSNR and time curves with measurement rates = 0.1 and 0.15 for different k:s. We can see from Table 3 and Figure 2 that the performance was the best when k:s was 3:7 considering from time and PSNR. Therefore, the kernel window and search window were set to be 3 and 7 in the following experiments.

The Influence of Non-Local Mean Regularization Kernel Window and Search Window
Non-local means regularization indicates the estimated value of the current pixel obtained by the weighted average of the pixels in the image that have a similar neighborhood structure. There is no doubt that the radius of the neighborhood kernel window and the radius of the neighborhood search window were vital to the experiment. If their values are too small, the self-similarity of the image cannot be fully utilized, and the characteristics of the non-local means cannot be used fully. If their values are too large, the image search area will become larger and the time will be longer, which leads to the fact that the algorithm efficiency will be lower. There are different kernel windows to search window ratio (k:s) result values for different measurement rates in Table 3. Figure  2 shows the PSNR and time curves with measurement rates = 0.1 and 0.15 for different k:s. We can see from Table 3 and Figure 2 that the performance was the best when k:s was 3:7 considering from time and PSNR. Therefore, the kernel window and search window were set to be 3 and 7 in the following experiments.

Verify Fractional Order Existence Performance
According to Section 3.1, we set fractional-order α = 1.7 in this experiment, in order to verify the validity of fractional-order existence. On the basis of Section 3.1.1, the fractional-order α = 1.7 was added. The test image selection was the same as the previous tested image Lena. The results are shown in Table 4.

Verify Fractional Order Existence Performance
According to Section 3.1, we set fractional-order α = 1.7 in this experiment, in order to verify the validity of fractional-order existence. On the basis of Section 3.1.1, the fractionalorder α = 1.7 was added. The test image selection was the same as the previous tested image Lena. The results are shown in Table 4. Comparing Table 1 with Table 4, we could see that the image reconstructed results in Table 4 were overall higher. The performance of fractional-order existence was better than that without fractional-order.

Verify High-Order Existence Performance
Similarly, in this experiment. In Section 3.1.2, the weight of high-order TV was 0.7. In order to verify the validity of high-order existence. On the basis of Section 3.1.2, the highorder was added. The test image was Lena. The results are shown in Table 5. From comparing Table 2 with Table 5, we could see that the performance of high-order existence was better than that without high-order by comparing data one by one.
In order to verify the effect of the proposed algorithm, we gave the visual comparison of different TV modes in algorithms to test images Lena and Peppers. From Figures 3 and 4, traditional TV with nonlocal means regularization in Figures 3a and 4a usually produced undesirable staircase artifact and painting-like effects. Even though fractional-order TV in Figures 3b and 4b could enhance image edges and textures, it could cause non-smooth of smooth areas. Although, the high-order TV in Figures 3c and 4c was capable of alleviating the problem caused by traditional TV, still smoothed out the details of the image, which led to the fact that the effect was not ideal. In Figures 3d and 4d, this was obvious that our proposed algorithm made the fractional-order and the high-order TV complement each other on image reconstruction. From the figures, we could see that our proposed algorithm preserved image edges and textures more effectively, and could alleviate the staircase artifact and painting-like effects produced by nonlocal means regularization. Electronics 2021, 10, x FOR PEER REVIEW 12 of 19

Comparison to Other Reconstruction Algorithms
In this experiment, we compared the proposed method with several state-of-the-art algorithms: TVNLR [18], BCS-TV [21], TVAL3 [28], and two non-TV based algorithms: BCS-SPL [30] and MH-BCS [31]. In order to reduce the computational complexity and memory requirements. Images were divided into non-overlapping blocks of size 32 × 32 for all algorithms. Table 6 displays their PSNR and SSIM values with the growth of the measurement rates for each image. We added additional test images and reconstruction results in the Appendix A. Figures 5 and 6 display PSNR and SSIM values respectively with the growth of the measurement rates for two images. From Table 6 and Figures 5  and 6, it can be seen that the HoFrTV algorithm performed best for the images at different measurement rates. The visual quality of the reconstructed images was used to further verify the effectiveness of the proposed algorithm. To compare the results visually, Figures 7 and 8 optionally display some reconstructed images and local enlargements obtained using the different algorithms with a 0.2 measurement rate. We can clearly see that the visual quality of the images recovered by HoFrTV was better than that of others. HoFrTV could efficiently reconstruct fine details and textures while preserving sharp edges and avoid producing painting-like effects.

Comparison to Other Reconstruction Algorithms
In this experiment, we compared the proposed method with several state-of-the-art algorithms: TVNLR [18], BCS-TV [21], TVAL3 [28], and two non-TV based algorithms: BCS-SPL [30] and MH-BCS [31]. In order to reduce the computational complexity and memory requirements. Images were divided into non-overlapping blocks of size 32 × 32 for all algorithms. Table 6 displays their PSNR and SSIM values with the growth of the measurement rates for each image. We added additional test images and reconstruction results in the Appendix A. Figures 5 and 6 display PSNR and SSIM values respectively with the growth of the measurement rates for two images. From Table 6 and Figures 5  and 6, it can be seen that the HoFrTV algorithm performed best for the images at different measurement rates. The visual quality of the reconstructed images was used to further verify the effectiveness of the proposed algorithm. To compare the results visually, Figures 7 and 8 optionally display some reconstructed images and local enlargements obtained using the different algorithms with a 0.2 measurement rate. We can clearly see that the visual quality of the images recovered by HoFrTV was better than that of others. HoFrTV could efficiently reconstruct fine details and textures while preserving sharp edges and avoid producing painting-like effects.

Comparison to Other Reconstruction Algorithms
In this experiment, we compared the proposed method with several state-of-the-art algorithms: TVNLR [18], BCS-TV [21], TVAL3 [28], and two non-TV based algorithms: BCS-SPL [30] and MH-BCS [31]. In order to reduce the computational complexity and memory requirements. Images were divided into non-overlapping blocks of size 32 × 32 for all algorithms. Table 6 displays their PSNR and SSIM values with the growth of the measurement rates for each image. We added additional test images and reconstruction results in the Appendix A. Figures 5 and 6 display PSNR and SSIM values respectively with the growth of the measurement rates for two images. From Table 6 and Figures 5 and 6, it can be seen that the HoFrTV algorithm performed best for the images at different measurement rates. The visual quality of the reconstructed images was used to further verify the effectiveness of the proposed algorithm. To compare the results visually, Figures 7 and 8 optionally display some reconstructed images and local enlargements obtained using the different algorithms with a 0.2 measurement rate. We can clearly see that the visual quality of the images recovered by HoFrTV was better than that of others. HoFrTV could efficiently reconstruct fine details and textures while preserving sharp edges and avoid producing painting-like effects.

Computational Complexity
The experiments were performed under the MATLABR2018a environment with Intel Corei5-4200 CPU of 3.4 GHz and 4.0 GB RAM. Table 7 and Figure 9 give the time of reconstructing image at different measurement rates. From this result, we could see that BCS-SPL, TVAL3, and MH-BCS were faster than the others. Comparing with TVNLR, HoFrTV, and BCS-TV, HoFrTV was faster than BCS-TV and slower than the TVNLR algorithm. It is because at each iteration needs to calculate the high-order and fractionalorder total variation in the reconstruction process. However, the HoFrTV algorithm had higher reconstruction quality.

Computational Complexity
The experiments were performed under the MATLABR2018a environment with Intel Corei5-4200 CPU of 3.4 GHz and 4.0 GB RAM. Table 7 and Figure 9 give the time of reconstructing image at different measurement rates. From this result, we could see that BCS-SPL, TVAL3, and MH-BCS were faster than the others. Comparing with TVNLR, HoFrTV, and BCS-TV, HoFrTV was faster than BCS-TV and slower than the TVNLR algorithm. It is because at each iteration needs to calculate the high-order and fractional-order total variation in the reconstruction process. However, the HoFrTV algorithm had higher reconstruction quality.

Conclusions
In this paper, a hybrid high-order and fractional-order total variation with nonlocal means regularization model was proposed for compressive sensing image reconstruction. ALM and ADMM methods were used to solve this model. Our proposed algorithm makes the fractional-order and the high-order total variation complement each other on image reconstruction. It can solve the problem of non-smooth in smooth areas when fractional-order total variation can enhance image edges and textures. In addition, it also addresses high-order total variation alleviates the staircase artifact produced by traditional total variation, still smoothes the details of the image and the effect is not ideal. Meanwhile, the proposed algorithm suppresses painting-like effects produced by nonlocal means regularization. Experimental results show that, comparing with several state-of-the-art algorithms, the reconstructed images obtained by the proposed approach not only outperformed many existing methods in terms of PSNR and SSIM but also had better visual quality.

Conclusions
In this paper, a hybrid high-order and fractional-order total variation with nonlocal means regularization model was proposed for compressive sensing image reconstruction. ALM and ADMM methods were used to solve this model. Our proposed algorithm makes the fractional-order and the high-order total variation complement each other on image reconstruction. It can solve the problem of non-smooth in smooth areas when fractional-order total variation can enhance image edges and textures. In addition, it also addresses highorder total variation alleviates the staircase artifact produced by traditional total variation, still smoothes the details of the image and the effect is not ideal. Meanwhile, the proposed algorithm suppresses painting-like effects produced by nonlocal means regularization. Experimental results show that, comparing with several state-of-the-art algorithms, the reconstructed images obtained by the proposed approach not only outperformed many existing methods in terms of PSNR and SSIM but also had better visual quality.  Funding: This work was supported by the National Natural Science Foundation of China (NSFC) (61675184 and 61405178).
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Figure A1. Five additional test images.