Adaptive Weighted High Frequency Iterative Algorithm for Fractional-Order Total Variation with Nonlocal Regularization for Image Reconstruction

We propose an adaptive weighted high frequency iterative algorithm for a fractional-order total variation (FrTV) approach with nonlocal regularization to alleviate image deterioration and to eliminate staircase artifacts, which result from the total variation (TV) method. The high frequency gradients are reweighted in iterations adaptively when we decompose the image into high and low frequency components using the pre-processing technique. The nonlocal regularization is introduced into our method based on nonlocal means (NLM) filtering, which contains prior image structural information to suppress staircase artifacts. An alternating direction multiplier method (ADMM) is used to solve the problem combining reweighted FrTV and nonlocal regularization. Experimental results show that both the peak signal-to-noise ratios (PSNR) and structural similarity index (SSIM) of reconstructed images are higher than those achieved by the other four methods at various sampling ratios less than 25%. At 5% sampling ratios, the gains of PSNR and SSIM are up to 1.63 dB and 0.0114 from ten images compared with reweighted total variation with nuclear norm regularization (RTV-NNR). The improved approach preserves more texture details and has better visual effects, especially at low sampling ratios, at the cost of taking more time.


Introduction
Compressed sensing (CS) [1,2] is an emerging framework for data acquisition and reconstruction, which permits us to reconstruct the original sparse or compressible signals from only a small number of linear measurements. CS has been exploited in image processing, such as 3D video [3], medical imaging [4], single-pixel imaging [5]. This is based on the principle that, through optimization, the sparsity of a signal can be recovered from far fewer samples than required by the Nyquist-Shannon sampling theorem when the measurement matrix satisfies the restricted isometry property (RIP) [6]. There is an attractive advantage that CS-based methods decrease data storage and transmission costs significantly for systems requiring large data. The most notable application, the single-pixel imaging system, reconstructed images from only a small amount of data from a single photodetector, with the result that it implements the mixing of image signals with a random mask such as the Hadamard matrix generated by a digital micromirror device (DMD) [7].
More specifically, the CS model reconstructs the image x ∈ R n from measurements y and it can be expressed as: min x Ψx s.t. y = Ax (1) where Ψ is transformation domain, A ∈ R M×N represents measurement matrices. The image x can be recovered by an inverse problem:x = argmin where Ψ(x) is a sparse prior, µ is the regularization parameter, which can control the trade-off between the regularization and data fidelity. Sparse prior knowledge [8] plays an important role in signal reconstruction. Generally, the current CS algorithms exploit the prior knowledge of original images under some suitable transformation such as DCT [9], wavelets [10], and total variation model [11,12]. The most popular approach is TV owing to the advantages of preserving image edges and reconstruction performance, which can be defined as: where D x and D y denote horizontal and vertical difference operators, respectively. TV models sets the same penalty for all gradients, which fails to preserve image details accurately and often leads to undesirable serious staircase artifacts. So many improved solvers have been proposed such as TVAL3 [13], TVNLR [14], RTV-NNR [15]. Li CB et al. proposed TVAL3, introducing an augmented Lagrangian method and an alternating technique with a nonmonotone line search to preserve edges. Zhang et al. [14] used an improved nonlocal regularization constraint operator, reducing reconstruction errors significantly by averaging the weights of TV regularization. Candes [16] proposed a reweighted TV, which penalized the gradients in each iteration adaptively with the weights w i . The model can be expressed as: where weights w can be defined as: Weights w are updated by x in each iteration, setting the different penalty for different regions. The parameter ε is a small positive constant to avoid division by zero. The regions which have large gradients (e.g., texture details) have small penalty and the others have large penalty. This method preserves the image edges effectively.
To improve the quality of reconstructed images, another strategy is nonlocal regularization based on nonlocal means (NLM) filtering. This strategy is effective for preserving image details and sharp edges by exploiting structural information. Dong [17] proposed a nonlocal low-rank regularization method using a smooth surrogate function for the rank as structured sparsity instead of the convex nuclear norm. This method combines nonlocal self-similarity and low rank to eliminate redundant information and artifacts. However, the reconstructed images are over-smoothed owing to the average tendency of different similar patches.
In this paper, we propose an adaptive reweighted fractional-order TV with nonlocal regularization for image reconstruction. This approach improves the TV model and only weights the high frequency gradients by extracting the high frequency components (e.g., texture details) from the images using pre-processing technique [18]. Fractional-order differential operators are introduced into the TV model replacing integer-order differential operators, enhancing high frequency components of images. We adopted the Grünwald-Letnikov (G-L) model using four different directions to handle the fractional-order gradients. Prior knowledges are exploited by introducing the nonlocal regularization constraint as structural information. An efficient augmented Lagrangian is developed to solve the above problem. ADMM is used to decompose this objective function into four sub-problems. We evaluate reconstruction performance by using peak a signal-to-noise (PSNR) and structural similarity index (SSIM) compared with other TV-based CS reconstruction methods.
We decompose the low and high frequency components, analyze the fractional differential model with different directions and give the definitions of nonlocal regularization in Section 2 of this paper. In Section 3, we describe equations of the proposed models. In Section 4, experimental results demonstrate the effectiveness of the proposed models. In Section 5, we give the conclusion.

Fractional-Order Differential Model
We introduce a fractional-order differential operator into TV-based algorithm. A fractional-order gradient, regarded as a generalization of the integer-order gradient, is composed of the fractional-order derivative at different directions. Here, we use the Grünwald-Letnikov (G-L) model [19,20] for image reconstruction and this model can be defined as: where v is fractional orders of function, t and a are the upper and lower boundaries of independent variables respectively, and h is the differential step-size. The total variation model of the image X can be viewed as the sum of two-dimensional discrete signal gradients in Equation (3). Fractional-order differential operators are introduced into the TV model, and the new model can be expressed as: where D v x , D v y represent fractional-order gradients in horizontal and vertical directions respectively: More specifically, we use four different directions to handle fractional-order gradients simply, corresponding to negative x and y axes, positive x and y axes, which the approximate extensive backward difference and the forward difference can be deduced easily. If the higher accuracy is required, more different directions (i.e., eight directions or sixteen directions) can be used with time costs permitting. In Equation (8), the corresponding coefficients C v k are expressed as: According to Equation (9), four different expressions can be written as: The fractional differential model enhances high frequency components more effectively, which preserves image details meanwhile losing some low frequency components.

Adpative Reweighted Total Variation Model
Natural images can be decomposed into smooth regions (low frequency components) and texture details (high frequency components), such as Lena in Figure 1. There is a fuzzy contour in Figure 1a, and it contains the most of image energy. Sharp edge textures can be seen in Figure 1b, which is crucial for visual effects. Equation (4) sets the different penalty for different gradients to preserve image edges without considering the structural information that results in the false textures and artifacts.
The fractional differential model enhances high frequency components more effectively, which preserves image details meanwhile losing some low frequency components.

Adpative Reweighted Total Variation Model
Natural images can be decomposed into smooth regions (low frequency components) and texture details (high frequency components), such as Lena in Figure 1. There is a fuzzy contour in Figure 1a, and it contains the most of image energy. Sharp edge textures can be seen in Figure 1b, which is crucial for visual effects. Equation (4) sets the different penalty for different gradients to preserve image edges without considering the structural information that results in the false textures and artifacts. We propose a new adaptive reweighting strategy for total variation model to solve the problem which loses the high frequency parts of the images. Decomposing low and high frequency components is a critical technique for the new strategy. Image x is decomposed into the smooth regions (low frequency) L x and image details (high frequency) H x . We only weight the high frequency components H x in iterations. The new TV model can be defined as: (11) To extract the low frequency components of images, we can solve the following deconvolution problem: where L f is a 3 3 × low pass filter and all coefficients are 1/9, ⊗ denotes convolution operation.  We propose a new adaptive reweighting strategy for total variation model to solve the problem which loses the high frequency parts of the images. Decomposing low and high frequency components is a critical technique for the new strategy. Image x is decomposed into the smooth regions (low frequency) x L and image details (high frequency) x H . We only weight the high frequency components x H in iterations. The new TV model can be defined as: To extract the low frequency components of images, we can solve the following deconvolution problem: where f L is a 3 × 3 low pass filter and all coefficients are 1/9, ⊗ denotes convolution operation. Z L is a low frequency feature map, and g d = [1, − 1] is the horizontal and vertical gradient operator. κ is a user-defined parameter. The solution of Equation (12) can solved by fast Fourier transform (FFT): where F and F −1 is FFT and inverse FFT, ' * ' denotes complex conjugate, and '•' denotes component-wise multiplication. Therefore, the low and high frequency components can be expressed as: We only employ the high frequency components to update the weights:

Nonlocal Regularization Model
For nature images, there are lots of same image patches in different regions and the same patches may be far apart in space domain, that is the nonlocal similarity. The similarity between two patches depends on the similarity of intensity gray levels by measuring Gaussian weighted Euclidean distances. We consider a nonlocal regularization model [21,22] based on nonlocal self-similarity to obtain more appropriate weight coefficients.
Given an image x = x(i)|i ∈ Ω , the estimate valuex(i) for a pixel i can be defined as: where S(i) is a i-centered search window L s × L s , x( j) is the intensity at pixel j. N i and N j denote the center pixel of similarity window d s × d s . The weights w ij depending on the similarity between gray value vectors P(N i ) and P(N j ), can be defined as: where h is an attenuation factor, altering the decay of the exponential function, and α is a standard deviation for Gaussian kernels. Z(i) is a normalizing constant, defined as: The nonlocal regularization (NR) model is expressed as:

Reweighted Fractional-Order TV Method with Nonlocal Regularization
In this section, we consider two strategies to improve the TV method for solving texture deficiency. First, low frequency and high frequency components of images are decomposed by Equations (13) and (14). Fractional differential operators are introduced into the TV model by using four different directions to handle fractional-order gradients approximately. We defined the FrTV weights by using Equation (15) and only adaptively reweight the high frequency gradients. Second, the nonlocal regularization constraint is employed as structural information to eliminate staircase artifacts.
The proposed method can be expressed as the optimization problem: where λ and β are regularization parameters, w is the high frequency parts weights and W is a regularization matrix. It is difficult to solve directly due to the non-differentiability. To solve this problem efficiently, we add constrained auxiliary variables: The corresponding augmented Lagrangian converting constrained Equation (21) into unconstrained objective functions: where a, b, c are the Lagrangian multipliers, α, γ 1 , γ 2 are corresponding meta-parameters. We solve the Equation (22) by solving Equations (23) and (24) iteratively: Alternating direction method of multipliers (ADMM) is drawn into decomposing four sub-problems using a variable splitting technique. We update each parameter iteratively until convergence. Assuming other parameters are fixed, u problem can be expressed as: Equation (25) has the closed form solution. Derivative Equation (25) and make derivative zero, we can obtain the solution: Fixing u, x, Z H , the sub-problem Z L is equivalent to: According to the Shrinkage-like lemma, solving Equation (22) with respect to Z L gives a closed-form solution at the t-th iteration: The sub-problem Z H has the similar expression with the same lemma: According to Equation (14) and fixing the other parameters, the last sub-problem x can be expressed as: For Equation (30), taking the derivative and finding zero of the derivatives, we can get the solution: Calculating Equation (31) directly will cost much time, so we find the solution by conjugate gradient method. The complete procedure of the improved method is shown in Algorithm 1.

Experimental Results and Analysis
The improved algorithm is compared with the other CS algorithms by evaluating the reconstructed performance on ten standard natural images with size 256 × 256, which are from the university of southern California image library. The original images are shown in Figure 2. In general, we quantify the reconstruction quality in terms of PSNR and SSIM. PSNR and SSIM are defined as: where MSE is the mean square errors between original images I(i, j) and reconstructed images K(i, j). µ x and µ y are gray mean of original images and reconstructed images, σ x and σ y are variance of original images and reconstructed images, σ xy is covariance. C 1 and C 2 are constant.
Electronics 2020, 9, x FOR PEER REVIEW 8 of 15 In our experiments, a Gaussian random matrix is used as a measurement matrix. We select the most appropriate parameters as follows: Lagrangian multipliers a, b, c are initialized to zero matrices. The proposed algorithm is compared with four algorithms, i.e., BCS-TV [23], TVNLR, NLR-CS, RTV-NNR. We obtain the experimental results at various sampling ratios (5%, 10%, 15%, 20%, 25%, 30%). All the experiments are performed on the Lenovo computer with Inter(R) Core (TM) i5-10210U CPU(1.6 GHz) and 16 G memory, running Windows 10 and Matlab 2012a.

The Influence of Fractional-Order v
The optimization of the parameters plays a vital role in image processing to obtain satisfied reconstructed results within acceptable time. Fractional-order differential operators can enhance high frequency components comparing with integer-order differential model. We analyze the effect of the parameter v on reconstruction performance to get the best quality. The image Barbara is used to evaluate the performance at 30% sampling ratios and several different fractional orders are chosen. The results are shown in Figure 3 and Figure 4. For each v, we obtain the reconstructed images and the residual images. It is obviously that when 1 v < , images reconstructed by the proposed method lost many texture details and the contour can be seen in residual images. When 1 v = , this method computes the integer order gradients and cannot preserve details effectively. When 1 v > , high frequency components are enhanced and reconstructed images have better visual effects. We cannot see the contour in residual images. However, when v is close to 2, PSNR reduce significantly and reconstructed images are fuzzier because texture details are enhanced excessively and become the noise. To obtain the best results, we set 1

The Influence of Fractional-Order v
The optimization of the parameters plays a vital role in image processing to obtain satisfied reconstructed results within acceptable time. Fractional-order differential operators can enhance high frequency components comparing with integer-order differential model. We analyze the effect of the parameter v on reconstruction performance to get the best quality. The image Barbara is used to evaluate the performance at 30% sampling ratios and several different fractional orders are chosen. The results are shown in Figures 3 and 4. For each v, we obtain the reconstructed images and the residual images. It is obviously that when v < 1, images reconstructed by the proposed method lost many texture details and the contour can be seen in residual images. When v = 1, this method computes the integer order gradients and cannot preserve details effectively. When v > 1, high frequency components are enhanced and reconstructed images have better visual effects. We cannot see the contour in residual images. However, when v is close to 2, PSNR reduce significantly and reconstructed images are fuzzier because texture details are enhanced excessively and become the noise. To obtain the best results, we set v = 1.4 in our experiments.
Electronics 2020, 9, x FOR PEER REVIEW 8 of 15 In our experiments, a Gaussian random matrix is used as a measurement matrix. We select the most appropriate parameters as follows: κ , ε = 0.1 empirically. The Lagrangian multipliers a, b, c are initialized to zero matrices. The proposed algorithm is compared with four algorithms, i.e., BCS-TV [23], TVNLR, NLR-CS, RTV-NNR. We obtain the experimental results at various sampling ratios (5%, 10%, 15%, 20%, 25%, 30%). All the experiments are performed on the Lenovo computer with Inter(R) Core (TM) i5-10210U CPU(1.6 GHz) and 16 G memory, running Windows 10 and Matlab 2012a.

The Influence of Fractional-Order v
The optimization of the parameters plays a vital role in image processing to obtain satisfied reconstructed results within acceptable time. Fractional-order differential operators can enhance high frequency components comparing with integer-order differential model. We analyze the effect of the parameter v on reconstruction performance to get the best quality. The image Barbara is used to evaluate the performance at 30% sampling ratios and several different fractional orders are chosen. The results are shown in Figure 3 and Figure 4. For each v, we obtain the reconstructed images and the residual images. It is obviously that when 1 v < , images reconstructed by the proposed method lost many texture details and the contour can be seen in residual images. When 1 v = , this method computes the integer order gradients and cannot preserve details effectively. When 1 v > , high frequency components are enhanced and reconstructed images have better visual effects. We cannot see the contour in residual images. However, when v is close to 2, PSNR reduce significantly and reconstructed images are fuzzier because texture details are enhanced excessively and become the noise. To obtain the best results, we set 1

Experimental Results
We evaluate the performance of the proposed model by comparing with other four algorithms. The PSNR values by using several algorithms are listed in Table 1. It is obviously that, for each image, BCS-TV has the worst reconstruction performance. NLR-CS, RTV-NNR and our proposed method are much better than other two methods. The proposed algorithm achieves significant performance improvements outperforming the four other algorithms and the gains are up to 7.54 dB, 7.15 dB, 3.41 dB and 1.63 dB respectively in Table 1. In the low sampling ratios, the proposed algorithm recovers more image information than others and gets the best visibility.

Experimental Results
We evaluate the performance of the proposed model by comparing with other four algorithms. The PSNR values by using several algorithms are listed in Table 1. It is obviously that, for each image, BCS-TV has the worst reconstruction performance. NLR-CS, RTV-NNR and our proposed method are much better than other two methods. The proposed algorithm achieves significant performance improvements outperforming the four other algorithms and the gains are up to 7.54 dB, 7.15 dB, 3.41 dB and 1.63 dB respectively in Table 1. In the low sampling ratios, the proposed algorithm recovers more image information than others and gets the best visibility.

Experimental Results
We evaluate the performance of the proposed model by comparing with other four algorithms. The PSNR values by using several algorithms are listed in Table 1. It is obviously that, for each image, BCS-TV has the worst reconstruction performance. NLR-CS, RTV-NNR and our proposed method are much better than other two methods. The proposed algorithm achieves significant performance improvements outperforming the four other algorithms and the gains are up to 7.54 dB, 7.15 dB, 3.41 dB and 1.63 dB respectively in Table 1. In the low sampling ratios, the proposed algorithm recovers more image information than others and gets the best visibility. We arbitrary choose four images (Lena, House, Boats and Chart2) and plot the PSNR curves as shown in Figure 5. At low sampling ratios, our algorithm has achieved the higher PSNR and improvements. When sampling ratios increase, the performance of RTV-NRR is closed to the proposed algorithm. The reason for this trend is that the magnitude of high frequency components is enhanced effectively by our proposed algorithm, while the magnitude of low frequency components is reduced, which fortunately has little effect on the visibility of the images.
We arbitrary choose four images (Lena, House, Boats and Chart2) and plot the PSNR curves as shown in Figure 5. At low sampling ratios, our algorithm has achieved the higher PSNR and improvements. When sampling ratios increase, the performance of RTV-NRR is closed to the proposed algorithm. The reason for this trend is that the magnitude of high frequency components is enhanced effectively by our proposed algorithm, while the magnitude of low frequency components is reduced, which fortunately has little effect on the visibility of the images. Intuitively, the visual results of the two images are shown in Figure 6 and the residual images (difference between original images and estimated images) are shown in Figure 7. Selecting portions of reconstructed images in red boxes and enlarging them, we compare them at 5% sampling ratios to appear significant differences. The quality of images reconstructed by BCS-TV and TVNLR are much worse than from the other three algorithms. Many block artifacts can be seen in Figure 6a,b. From Figure 7a,b, the smooth regions are preserved well but the texture details are lost a lot, therefore BCS-TV and TVNLR have bad visual effects. Images reconstructed by NLR-CS are better and have better visual effects. But this algorithm suffers from over-smoothed effects owing to the average tendency of different similar patches. Our proposed algorithm has the best visibility compared with other algorithms in Figure 6. This method preserves image texture details and eliminates staircase artifacts by using the fractional-order differential method and nonlocal regularization. We can hardly see the intact contour in residual images in Figure 7e. Experimental results demonstrate that our proposed method enhances high frequency components (image details) significantly.
Moreover, we acquire SSIM for the sampling ratios of 5% in Table 2. Table 2 shows that BCS-TV and TVNLR are much worse than others and the SSIM of the proposed algorithm are better than the Intuitively, the visual results of the two images are shown in Figure 6 and the residual images (difference between original images and estimated images) are shown in Figure 7. Selecting portions of reconstructed images in red boxes and enlarging them, we compare them at 5% sampling ratios to appear significant differences. The quality of images reconstructed by BCS-TV and TVNLR are much worse than from the other three algorithms. Many block artifacts can be seen in Figure 6a,b. From Figure 7a,b, the smooth regions are preserved well but the texture details are lost a lot, therefore BCS-TV and TVNLR have bad visual effects. Images reconstructed by NLR-CS are better and have better visual effects. But this algorithm suffers from over-smoothed effects owing to the average tendency of different similar patches. Our proposed algorithm has the best visibility compared with other algorithms in Figure 6. This method preserves image texture details and eliminates staircase artifacts by using the fractional-order differential method and nonlocal regularization. We can hardly see the intact contour in residual images in Figure 7e. Experimental results demonstrate that our proposed method enhances high frequency components (image details) significantly.
Moreover, we acquire SSIM for the sampling ratios of 5% in Table 2. Table 2 shows that BCS-TV and TVNLR are much worse than others and the SSIM of the proposed algorithm are better than the other algorithms for every image. Compared with the second-best method RTV-NNR, our proposed algorithm outperforms by up to 0.0114.
To demonstrate that our proposed algorithm is better for enhancing high frequency components, we decompose the original images and reconstructed images into smooth regions and texture details by using Equations (13) and (14). The texture details of reconstructed images are shown in Figure 8.
other algorithms for every image. Compared with the second-best method RTV-NNR, our proposed algorithm outperforms by up to 0.0114.
To demonstrate that our proposed algorithm is better for enhancing high frequency components, we decompose the original images and reconstructed images into smooth regions and texture details by using Equations (13) and (14). The texture details of reconstructed images are shown in Figure 8.    other algorithms for every image. Compared with the second-best method RTV-NNR, our proposed algorithm outperforms by up to 0.0114. To demonstrate that our proposed algorithm is better for enhancing high frequency components, we decompose the original images and reconstructed images into smooth regions and texture details by using Equations (13) and (14). The texture details of reconstructed images are shown in Figure 8.      It is clear that our proposed algorithm preserves the most textures in Figure 8 and we can see the complete edge contour. The other four algorithms lost many details, resulting in the terrible visibility in Figure 8. Furthermore, we compare the average relative errors and SSIM of texture details in Figure 9. For every image and sampling ratio, our algorithm has the smallest relative errors and the best SSIM of high frequency components. These results demonstrate that our proposed algorithm enhances the texture details significantly. Compared with other methods, our algorithm lost some low frequency components, but fortunately it has little effect on the visibility of images. We also compare the average reconstruction time at different sampling ratios in Figure 10 and it is shown that our algorithm costs a little more time than TVNLR, NLR-CS and RTV-NNR, but less than BCS-TV, which we will improve in the future.  It is clear that our proposed algorithm preserves the most textures in Figure 8 and we can see the complete edge contour. The other four algorithms lost many details, resulting in the terrible visibility in Figure 8. Furthermore, we compare the average relative errors and SSIM of texture details in Figure 9. For every image and sampling ratio, our algorithm has the smallest relative errors and the best SSIM of high frequency components. These results demonstrate that our proposed algorithm enhances the texture details significantly. Compared with other methods, our algorithm lost some low frequency components, but fortunately it has little effect on the visibility of images. We also compare the average reconstruction time at different sampling ratios in Figure 10 and it is shown that our algorithm costs a little more time than TVNLR, NLR-CS and RTV-NNR, but less than BCS-TV, which we will improve in the future. It is clear that our proposed algorithm preserves the most textures in Figure 8 and we can see the complete edge contour. The other four algorithms lost many details, resulting in the terrible visibility in Figure 8. Furthermore, we compare the average relative errors and SSIM of texture details in Figure 9. For every image and sampling ratio, our algorithm has the smallest relative errors and the best SSIM of high frequency components. These results demonstrate that our proposed algorithm enhances the texture details significantly. Compared with other methods, our algorithm lost some low frequency components, but fortunately it has little effect on the visibility of images. We also compare the average reconstruction time at different sampling ratios in Figure 10 and it is shown that our algorithm costs a little more time than TVNLR, NLR-CS and RTV-NNR, but less than BCS-TV, which we will improve in the future.