A Novel Iterative Thresholding Algorithm Based on Plug-and-Play Priors for Compressive Sampling

: We propose a novel fast iterative thresholding algorithm for image compressive sampling (CS) recovery using three existing denoisers—i.e., TV (total variation), wavelet, and BM3D (block-matching and 3D ﬁltering) denoisers. Through the use of the recently introduced plug-and-play prior approach, we turn these denoisers into CS solvers. Thus, our method can jointly utilize the global and nonlocal sparsity of images. The former is captured by TV and wavelet denoisers for maintaining the entire consistency; while the latter is characterized by the BM3D denoiser to preserve details by exploiting image self-similarity. This composite constraint problem is then solved with the fast composite splitting technique. Experimental results show that our algorithm outperforms several excellent CS techniques.

To establish these complicated models for image CS recovery, penalty functions are frequently used to encourage solutions of a certain form [5,6,11].In FCSA (Fast Composite Splitting Algorithm) [5], sparsity penalty terms in wavelet and gradient domain are jointly used to constrain the solution space; while WaTMRI (wavelet tree sparsity magnetic resonance imaging) [6] superadds a tree sparse regularization to the objective function of FCSA, further forcing the parent-child wavelet coefficients to be zeros or non-zeros together.On the other hand, in Bayesian CS [8,12], graphical models are employed to describe the probabilistic relationship between measurements and the original signal.In [8], the structure of the wavelet coefficients is modeled by a hidden Markov tree (HMT).
Among these structures, nonlocal sparsity [11,12,15], which exploits the self-repetitive structure exhibited often in natural images has shown great potential.In NLR-CS (CS via nonlocal low-rank regularization) [11], a patch-based low-rank regularization model is built to enforce the low-rank property over the sets of nonlocal similar patches; while RCoS (CS recovery via collaborative sparsity) [14] uses a 3D sparsity term to maintain image nonlocal consistency.In contrast, the D-AMP (denoising-based approximate message passing) algorithm [12] exploits the image self-similarity through the use of nonlocal based denoisers, such as the BM3D (block-matching and 3D filtering) denoiser [16], which performs hard or soft thresholding with a 3D orthogonal dictionary (3D filtering) on 3D image blocks built by stacking similar patches together (block-matching).Closely relating to D-AMP, BM3D-CS [15] iteratively adds noise to the missing part of the spectra and then applies Future Internet 2017, 9, 24 2 of 10 BM3D to the result.However, the above signal models cannot simultaneously utilize the global and nonlocal sparsity.Although the nonlocal sparse regularization can preserve image fine details by clustering similar components, it is based on block matching, and thus introduces visible artifacts in homogeneous regions, manifesting as low-frequency noise [17].On the other hand, TV and wavelet sparse regularization can maintain the consistency of the entire image.
In this paper, we propose an iterative thresholding method that is competitive in quality with D-AMP, but is much simpler to implement.We combine three popular denoisers through the use of the plug-and-play prior approach [18,19], which is able to turn a denoiser into reconstruction solver.For charactering the nonlocal sparsity, we use the nonlocal wavelet denoiser.For representing the global sparsity, we use the wavelet and gradient denoiser.This multiple denoisers problem is then solved based on fast composite splitting technique [5] and the fast proximal method FISTA (fast iterative shrinkage-thresholding algorithm) [20,21].In the frame of FISTA, the original composite regularization problem is firstly decomposed into three simpler regularization sub-problems via fast composite splitting technique; then each of them is separately solved by thresholding methods.Experimental results show that our method impressively outperforms previous methods in terms of reconstruction accuracy, and also achieves competitive advantage in computation complexity.

Compressive Sampling and Fast Iterative Shrinkage-Thresholding Algorithm
The compressive sampling recovery is an underdetermined problem.More than one solution can yield the same measurements, but if x is sparse, we can recover the original signal by minimizing the following L1 problem.
where λ is a regularization parameter.
The above problem can be efficiently solved by various methods [22][23][24][25], one of which is ISTA (iterative shrinkage-thresholding algorithm) [26].Specifically, the general step of ISTA is summarized in Algorithm 1, where f = y − Ax 22 /2, and ∇ f (x k ) = A * (Ax k − y) denotes the gradient of the function f at the point x k ; A * denotes the conjugate transpose of A; c is a step size.The undersampled Fourier matrix A = SF is used as the measurement matrix, where F is the 2D Fourier transform and S is a selection matrix containing m rows of the identity matrix.Step (b) is a proximity operator.It has a close form solution, which can be expressed as: FISTA [21] offers an accelerated version of ISTA by adding an accelerated step (step (c)), which is summarized in Algorithm 2.
Algorithm 1. ISTA [26] Input: the CS measurements y and the measurement matrix A Initialization: FISTA [21] Input: the CS measurements y and the measurement matrix A Initialization:

The New Composite Model
Images are generally sparse in wavelet and gradient domain, although not strictly.Moreover, there exist abundant similar image patches in natural images.These sparse priors can be utilized by using three regularizations (TV, wavelet sparse, and nonlocal wavelet regularization).The nonlocal sparse regularization can preserve image fine details, while TV and wavelet sparse regularization can maintain the consistency of the entire image.In order to gain the advantages of these regularizations, we combine the nonlocal self-similarity of patches and the global sparsity in wavelet and gradient domain to jointly enforce the solution in image CS recovery, and formulate a composite regularization problem as: where α, β, and γ are three regularization parameters.The first item is the fidelity term; while the last three items are priori constraint terms, denoting TV norm, wavelet sparse norm, and nonlocal wavelet norm respectively.Usually, the TV norm of a 2D image X with size p × q is computed by: and the wavelet sparse norm can be represented as Ψx 1 where Ψ is the wavelet transform.However, these regularizations are equal to denoising operations in the solving process by adopting the plug-and-play prior approach.We characterize the image self-similarity by means of the last item in Equation (3) derived from BM3D [15], which is used to seek patch correlation in image denoising.
To construct x NL , we first stack similar image patches to 3D groups, then 3D wavelet transform is conduct on each 3D group, which is achieved by conducting 2D wavelet transform on each patch and then conducting 1D wavelet transform along the third axis.

Solving The Composite Model with Fast Composite Splitting Technique
We solve the above composite regularization problem Equation (3) based on FISTA algorithm and fast composite splitting technique [5] which decomposes the original composite problem into simpler regularization sub-problems, separately solves each of them, and then obtains the final solution by a linear combination.Algorithm 3 outlines the solving procedure of our algorithm named as PPP-CS (CS with plug-and-play priors).In step (a), the fidelity term is computed with gradient descent method to get an initial solution x g .In steps (b)-(d), variable x is split into three variables x 1 , x 2 and x 3 ; the corresponding proximity operator is performed over each of them.In step (e) the solution x k is obtained by linear combination of x 1 , x 2 , and x 3 .Step (f) is the acceleration step borrowed from FISTA.The convergence of PPP-CS is guaranteed by fast composite splitting technique and FISTA algorithm.

Algorithm 3. CS with Plug-and-play Priors
Input: the CS measurements y and the measurement matrix A Initialization: Steps (b)-(d) are the so-called proximal map, which is formally defined in [27], i.e., given a continuous convex function g(x) and any scalar ρ > 0, the proximal map associated to function g is defined as: In the plug-and-play prior approach [18,19], this proximal map is equal to a denoising operation, where x is a noisy image, g(u) is a regularization corresponding to a denoiser, and ρ is the denoising threshold.In the steps (b)-(d), we regard x g as a noisy image; then the corresponding denoising methods can be adopted to obtain the solution of each sub-problem.For step (b), we use a TV-based denoising method which is based on the FISTA algorithm [27].For step (c), the wavelet denoisers are available; but we can solve this sub-problem easily, since it has a close form solution, which can be computed with soft thresholding (see Equation ( 2)).For step (d), we use the BM3D denoiser [15].
The key to success in solving Equation ( 3) is to determine the denoising threshold of the BM3D denoiser.The original threshold of BM3D is proportional to the noise variance of image; however, it is hard to estimate the noise variance of the current image in the CS reconstruction.Apparently, the denoising threshold δ should be reduced along with the iteration as the recovered image becomes more and more clear.Therefore, a natural idea is that it should be proportional to the observation residual in CS reconstruction.
where s is a scale factor and can be set empirically.

Experiments
To evaluate reconstruction performances of the proposed algorithm in items of average PSNR (peak signal to noise ratio), runtime, and visual quality, we compare our algorithm with four excellent algorithms.There are two algorithms based on iterative shrinkage-thresholding method FCSA [1] (using wavelet sparsity and gradient sparsity), WaTMRI [2] (using wavelet sparsity, gradient sparsity and tree sparsity) and two Bayesian CS methods Turbo-AMP [3] (using tree sparsity), D-AMP [6] (using nonlocal sparsity via BM3D denoiser).Three experiments are carried on both natural images and MR (Magnetic Resonance) images with size 256 × 256 at four sampling ratios (18, 20, 22, and 25%).The eight natural images and four MR images used in our experiments are shown in Figure 1.The undersampled Fourier matrix is generated by following the previous works [2,3], which randomly choose more Fourier coefficients from low frequency and less on high frequency.All experiments are on a desktop with 3.80 GHz AMD A10-5800K CPU.Matlab version is 7.11 (R2010b).sparsity), D-AMP [6] (using nonlocal sparsity via BM3D denoiser).Three experiments are carried on both natural images and MR (Magnetic Resonance) images with size 256 × 256 at four sampling ratios (18, 20, 22, and 25%).The eight natural images and four MR images used in our experiments are shown in Figure 1.The undersampled Fourier matrix is generated by following the previous works [2,3], which randomly choose more Fourier coefficients from low frequency and less on high frequency.All experiments are on a desktop with 3.80 GHz AMD A10-5800K CPU.Matlab version is 7.11 (R2010b).We set the maximum iterations K = 10 for Turbo-AMP and set K = 50 for the rest algorithms.For our algorithm PPP-CS, we use the setting  = 0.001,  = 0.023, c = 1.The scale factor in Equation (6) s is set to 30 in the first 20 iterations, and is fixed at 5 after 20 iterations; since the estimated threshold  is basically steady at that moment.The 3D wavelet transform Г3D is composed of 2-D bior1.5 and 1-D Haar.To construct 3D groups by stacking similar patches, we need to set the following parameters: the size of each patch is 8 × 8 pixels the size of window for searching matched patches is 25 × 25 pixels; the number of best matched patches is 16; and the sliding step to process every next reference patch is set to 6.After the first five iterations, step (d) is only updated every three iterations in order to speed up the algorithm.For the rest algorithms, the default settings in their codes are used.

Quantative Evaluation
We carry all algorithms on the 12 test images with various sampling ratios, the average PSNR results by running each image five times are shown in Figure 2, from which we can see that the nonlocal sparsity-based methods PPP-CS and D-AMP are better than others.Particularly, the proposed algorithm PPP-CS achieves the highest PSNR results.By jointly using of the global and nonlocal sparsity, the proposed PPP-CS method performs much better than the D-AMP method on all sensing rates.The PSNR results of reconstructed images from the CS measurements with 20% sampling are included in Figure 3.One can see that the proposed PPP-CS method produces the highest PSNRs on almost all test images (except for three cases where PPP-CS slightly falls behind D-AMP).The average PSNR improvement over D-AMP is about 10.54 dB, which validates the superiority of the proposed algorithm in objective quality.
Figure 4 gives the average CPU time of different algorithms carried on the 12 test images by running each 5 times.FCSA spends the least CPU time; however, its PSNR results in Figure 1 are relatively poor.Note that PPP-CS spends 62.05% less CPU time to achieve higher PSNR than D-AMP, for the result that the Bayesian algorithm AMP is slower than our fast iterative shrinkage-thresholding algorithm.We set the maximum iterations K = 10 for Turbo-AMP and set K = 50 for the rest algorithms.For our algorithm PPP-CS, we use the setting α = 0.001,β = 0.023, c = 1.The scale factor in Equation ( 6) s is set to 30 in the first 20 iterations, and is fixed at 5 after 20 iterations; since the estimated threshold δ is basically steady at that moment.The 3D wavelet transform Г 3D is composed of 2-D bior1.5 and 1-D Haar.To construct 3D groups by stacking similar patches, we need to set the following parameters: the size of each patch is 8 × 8 pixels the size of window for searching matched patches is 25 × 25 pixels; the number of best matched patches is 16; and the sliding step to process every next reference patch is set to 6.After the first five iterations, step (d) is only updated every three iterations in order to speed up the algorithm.For the rest algorithms, the default settings in their codes are used.

Quantative Evaluation
We carry all algorithms on the 12 test images with various sampling ratios, the average PSNR results by running each image five times are shown in Figure 2, from which we can see that the nonlocal sparsity-based methods PPP-CS and D-AMP are better than others.Particularly, the proposed algorithm PPP-CS achieves the highest PSNR results.By jointly using of the global and nonlocal sparsity, the proposed PPP-CS method performs much better than the D-AMP method on all sensing rates.The PSNR results of reconstructed images from the CS measurements with 20% sampling are included in Figure 3.One can see that the proposed PPP-CS method produces the highest PSNRs on almost all test images (except for three cases where PPP-CS slightly falls behind D-AMP).The average PSNR improvement over D-AMP is about 10.54 dB, which validates the superiority of the proposed algorithm in objective quality.
Figure 4 gives the average CPU time of different algorithms carried on the 12 test images by running each 5 times.FCSA spends the least CPU time; however, its PSNR results in Figure 1 are relatively poor.Note that PPP-CS spends 62.05% less CPU time to achieve higher PSNR than D-AMP, for the result that the Bayesian algorithm AMP is slower than our fast iterative shrinkage-thresholding algorithm.Figure 5 gives the performance comparisons between different methods in terms of the CPU time over the PSNR with 20% sampling.The PPP-CS obtains the best reconstruction results by achieving the highest PSNR in less CPU time after five seconds.Since PPP-CS is based on FISTA, it can converge fast.The D-AMP is inferior to the PPP-CS for most of the time, because of its higher cost in each iteration.These results show the effectiveness of acceleration strategies in the PPP-CS inherited from the FISTA for the image reconstruction.
Future Internet 2017, 9, 24 7 of 10 Figure 5 gives the performance comparisons between different methods in terms of the CPU time over the PSNR with 20% sampling.The PPP-CS obtains the best reconstruction results by achieving the highest PSNR in less CPU time after five seconds.Since PPP-CS is based on FISTA, it can converge fast.The D-AMP is inferior to the PPP-CS for most of the time, because of its higher cost in each iteration.These results show the effectiveness of acceleration strategies in the PPP-CS inherited from the FISTA for the image reconstruction.

Visual Quality Evaluation
Figure 6 provides the CS recovery results for image Cameraman in the case of sample rate = 20% using the test algorithms.For each algorithm, we compute the PSNR as well as the structural similarity index (SSIM) which better reflects the visual quality of the images.From Figure 6, it is clear to see that the nonlocal sparsity-based methods PPP-CS and D-AMP are still better than others.The proposed algorithm enjoys great advantages over other competing algorithms in producing more clean screens, e.g., on the area of sky.Due to the joint use of the global and nonlocal sparsity and adaptive thresholding strategy, the elimination of artifacts is more effective by applying PPP-CS.

Visual Quality Evaluation
Figure 6 provides the CS recovery results for image Cameraman in the case of sample rate = 20% using the test algorithms.For each algorithm, we compute the PSNR as well as the structural similarity index (SSIM) which better reflects the visual quality of the images.From Figure 6, it is clear to see that the nonlocal sparsity-based methods PPP-CS and D-AMP are still better than others.The proposed algorithm enjoys great advantages over other competing algorithms in producing more clean screens, e.g., on the area of sky.Due to the joint use of the global and nonlocal sparsity and adaptive thresholding strategy, the elimination of artifacts is more effective by applying PPP-CS.
Future Internet 2017, 9, 24 7 of 10 Figure 5 gives the performance comparisons between different methods in terms of the CPU time over the PSNR with 20% sampling.The PPP-CS obtains the best reconstruction results by achieving the highest PSNR in less CPU time after five seconds.Since PPP-CS is based on FISTA, it can converge fast.The D-AMP is inferior to the PPP-CS for most of the time, because of its higher cost in each iteration.These results show the effectiveness of acceleration strategies in the PPP-CS inherited from the FISTA for the image reconstruction.The corresponding the iterations vs. PSNR curves and iterations vs. SSIM curves are given in Figures 7 and 8.Note that our algorithm achieves significant improvement in comparison with other algorithms both visually and quantitatively (PSNR and SSIM).The corresponding the iterations vs. PSNR curves and iterations vs. SSIM curves are given in Figures 7 and 8.Note that our algorithm achieves significant improvement in comparison with other algorithms both visually and quantitatively (PSNR and SSIM).The corresponding the iterations vs. PSNR curves and iterations vs. SSIM curves are given in Figures 7 and 8.Note that our algorithm achieves significant improvement in comparison with other algorithms both visually and quantitatively (PSNR and SSIM).

Conclusions
A new image compressive sampling recovery algorithm based on the FISTA algorithm and plug-and-play priors is proposed in this paper.Our contributions are as follows: (1) jointly using of the global and nonlocal sparsity to maintain entire consistency and preserve details in image CS recovery, which can be achieved through the use of the plug-and-play prior approach; (2) applying the FISTA algorithm and fast composite splitting technique to turn the complicated composite constraint problem into three independent shrinkage-thresholding problems, and obtain a simple implementation and a fast convergence rate.Finally, we combine the FISTA algorithm with plug-and-play priors for the first time, making it easy to use the powerful denoisers in CS reconstruction.In the experiments, our algorithm is shown to outperform four excellent algorithms in reconstructed quality, and also achieve competitive advantage in runtime.

Conclusions
A new image compressive sampling recovery algorithm based on the FISTA algorithm and plug-and-play priors is proposed in this paper.Our contributions are as follows: (1) jointly using of the global and nonlocal sparsity to maintain entire consistency and preserve details in image CS recovery, which can be achieved through the use of the plug-and-play prior approach; (2) applying the FISTA algorithm and fast composite splitting technique to turn the complicated composite constraint problem into three independent shrinkage-thresholding problems, and obtain a simple implementation and a fast convergence rate.Finally, we combine the FISTA algorithm with plug-and-play priors for the first time, making it easy to use the powerful denoisers in CS reconstruction.In the experiments, our algorithm is shown to outperform four excellent algorithms in reconstructed quality, and also achieve competitive advantage in runtime.

Figure 4 .
Figure 4. Average CPU time (s) comparisons on the 12 test images with various sampling ratios.

Figure 6
Figure6provides the CS recovery results for image Cameraman in the case of sample rate = 20% using the test algorithms.For each algorithm, we compute the PSNR as well as the structural similarity index (SSIM) which better reflects the visual quality of the images.From Figure6, it is clear to see that the nonlocal sparsity-based methods PPP-CS and D-AMP are still better than others.The proposed algorithm enjoys great advantages over other competing algorithms in producing more clean screens, e.g., on the area of sky.Due to the joint use of the global and nonlocal sparsity and adaptive thresholding strategy, the elimination of artifacts is more effective by applying PPP-CS.