Combined First- and Second-Order Variational Model for Image Compressive Sensing

Ahybrid variationalmodel combinedfirstand second-order total variation for image reconstruction from its finite number of noisy compressive samples is proposed in this paper. Inspired bymajorization-minimization scheme, we develop an efficient algorithm to seek the optimal solution of the proposed model by successively minimizing a sequence of quadratic surrogate penalties. Both the nature and magnetic resonance (MR) images are used to compare its numerical performance with four state-of-the-art algorithms. Experimental results demonstrate that the proposed algorithm obtained a significant improvement over related state-of-the-art algorithms in terms of the reconstruction relative error (RE) and peak signal to noise ratio (PSNR).


Introduction
Traditional approaches to sampling signals or images follow the basic principle of the Nyquist-Shannon sampling theorem; the sampling rate must be at least twice of the frequency bandwidth.In many practical applications, including image and video cameras, MRI scanners, and radar receivers, requirements of high resolution imaging lead to very high Nyquist sampling rate which brings a series of realistic difficult problems in the field of data conversion (e.g., analogdigital conversion), storage, and transmission.The technique called compressive sensing (CS), which goes against conventional wisdoms in data acquisition, has recently been developed to overcome those problems.
CS theory designs nonadaptive sampling techniques that condense the information in sparse or compressible images into a small amount of data and yet reconstruct them accurately.The general framework of compressive sensing consists of two phases: encoding and decoding.In the encoding phase, a sparse or compressible image u ∈   1 × 2 is encoded into noisy compressive samples by y = Φu + n. ( Here, Φ ∈  × is the sensing matrix ( is the number of measurements and  =  1 ×  2 is the number of total pixels,  ≪ ), y ∈   is the compressive samples of original image u, and n is Gaussian white noise with standard deviation .While in the decoding phase, a nonlinear recovery algorithm is used to reconstruct the original image u from compressive samples y [1][2][3].One common approach is to formulate the CS decoding as the following optimization problem [4]: where R(u) is the regularization term and  is a tunable regularization parameter that trades off the fidelity term ‖y − Φu‖ 2 2 and the regularization term R(u).One of popular choices of the regularization term is the total variation (TV) regularizer defined by the  1 norm of the image gradient [5,6]: Total variation regularization was firstly introduced by Rudin et al. [7] to cope with image denoising problem in 1992.Their approach had a significant impact in the area of inverse problems in image processing.Since then, various improvements have been proposed by the community, for example, adaptive TV regularization [8,9], nonlocal TV regularization [10], and anisotropic and directional TV regularization [11][12][13].TV regularization is widely used in image processing problems mainly due to its good properties such as convexity, invariance to image shifts, high-accuracy of recovery piecewise constant images, and desirable ability to preserve edges.However, it also has some limitations.The penalization of  1 norm of gradient usually encourages the recovery of images in a piecewise-constant form, thus, resulting in reconstructed images with patchy of painting like staircase artifacts [14,15].To reduce the staircase effect and give better denoising performance, several higher-order TV regularization schemes were reported in the context of image denoising over the last few years [14][15][16][17][18].However, the higher-order regularizers are prone to causing blurring across sharp edges, since they prefer to suppress large partial gradients.To complement each other for most nature images, it is thus needed to further improve the image restoration capability by combining first-and higher-order methods into one framework.
In this paper, we propose a hybrid compressive sensing method using a new hybrid TV measure by a combined first-and second-order TV penalty for recovering a piecewise smooth image containing all possible sharper edges from limited compressive samples.To seek the optimal solution of the proposed model, we develop an efficient image reconstruction algorithm by using the majorization-minimization scheme.The novelty in this work is that our hybrid TV regularization method is able to utilize the best properties of first-and second-order TV regularization and manage to overcome the weaknesses of both.
The rest of the paper is organized as follows.Section 2 provides a brief review of second-order TV regularization method for CS image reconstruction problem.In Section 3, we give a detailed description of our hybrid TV-based CS construction model and corresponding optimization algorithm.Experimental results and performance analysis are given in Section 4. Finally, Section 5 concludes this paper.

Second-Order TV Penalty
In this section, we reinterpret classic TV regularizer and extend it to second-order TV regularizer.We denote firstorder directional derivative of u along the unit vector e  = (cos , sin ) as u ,1 , which can be expressed as the linear combination of horizontal and vertical directional derivatives as follows: Similarly, for the second-order directional derivative u ,2 , we have In order to extend the standard TV scheme to higherorder TV regularizer, we reinterpret the TV regularizer as a group separable  1 −  2 mixed norm of directional derivatives of the specified image.To do this, we rewrite gradient magnitude as follows: Hence, TV regularizer defined in (3) can be expressed as a functional involving the directional derivatives of u as follows: Based on the above reinterpretation, we introduce second-order TV regularizer by replacing first-order directional derivative u ,1 in (7) with second-order directional derivative u ,2 , thus we have Substituting ( 5) into (8), we can rewrite the second-order TV regularizer R 2 (u) as follows: Since the second-order TV regularizer sums the square magnitude of the directional derivatives of the image along all directions and orientations, this penalty is invariant to rotations and translations and is also convex.The secondorder derivatives are steerable which enables us to obtain analytical expressions for the new regularizer as the standard TV regularizer.Furthermore, the minimization of secondorder TV regularizer will preserve the strong directional derivatives and attenuate the small ones at other directions; thus, it can provide better preservation along line-like features.

The Proposed Variational Approach for
Image Compressive Sensing

Hybrid TV-Based Reconstruction Model.
Combining first-order and second-order TV regularizers, we present a hybrid TV regularization model to reconstruct original image from its noisy compressive samples.Given noisy compressive samples y and sensing matrix Φ, the proposed reconstruction model can be formulated as where  and  are two positive parameters that control contribution of the fidelity term and two regularization terms.
In the first-order TV regularization method, images are often assumed to be piecewise-constant, while the second-order TV regularization method implicitly adopts a piecewise-linear assumption.The piecewise-linear approximation of images usually yields a lesser approximation error than piecewise-constant representation.Therefore, assuming that the images under consideration can be better decomposed into piecewise-constant and piecewise-linear functions, we expect the combination of first-and second-order regularization to be more accurate in terms of reconstruction quality.

Discretization of Hybrid TV Regularizer.
Let us now consider the discrete formulation of two regularization terms in (10).A discrete version can be obtained by considering that the image u has been uniformly sampled.At this point, we switch the notation and start assuming that u denote vector containing all these pixels arranged in column lexicographic ordering.
Using local first-order differences to approximate the two orthogonal components of the gradient, we obtain a discrete version of first-order TV regularizer.Consider where  ℎ  and  V  are operators corresponding to horizontal and vertical first-order differences at pixel , respectively, and the summation ∑  =1 takes over all pixels.Similarly, we can write discrete version of second-order TV regularizer as Here  ℎℎ  ,  VV  , and  ℎV  are operators corresponding to, respectively, horizontal-horizontal, vertical-vertical, horizontal-vertical second-order differences at pixel .
For the sake of simplicity, we denote Thus, (10) can be rewritten as

An Efficient Algorithm.
In this subsection, we derive an efficient optimization algorithm which belongs to the class of majorization-minimization (MM) approaches [19,20] to solve proposed model (10).The MM approaches have been introduced to solve optimization problems that are too difficult to solve directly.These schemes reformulate the original problem as the solution to a sequence of surrogate problems.Optimizing the surrogate functions will drive the objective function downward until a local optimum is reached.
Let (x) be the cost function to be minimized.Instead of minimizing (x) directly, the MM approaches solve a sequence of surrogate problems {  (x),  = 0, 1, 2, . ..} which are solved easier than (x).As the result, a sequence of {x  ,  = 1, 2, . ..} obtained by Lena Pepper Barbara Brain are then applied to estimate the solution of (x).In (15),   (x) is the majorizer of function (x) at a fixed point x  which satisfies the following two properties: To develop a MM-based algorithm, we derive a quadratic majorizer for the proposed hybrid TV regularizer.The choice of a quadratic majorizer is motivated by the fact that the minimization of a quadratic function amounts to solving a system of linear equations, a task for which there are excellent methods available in the literature.We define the following majorizer: Using the inequality √ ≤ √   + ( −   )/(2 √   ), we have for any u, with equality for u  .Equation ( 19) means that majorizer   (u) satisfies the two properties ( 16) and (17).So, the solution of optimization problem ( 14) can be obtained by finding the minimizer of   (u),  = 0, 1, 2, . ... Consider Notice that, for fixed , the terms   1 (u  ) and   2 (u  ) in ( 18) are simply additive constants which can be ignored since they do not affect the solution of the optimization problem.Thus, problem (20) is equivalent to Let  ℎ and  V denote the horizontal and vertical global finite difference matrices with size  ×  such that  ℎ u and  V u are the vectors of all horizontal and vertical first-order differences of image u.We then have where here, w (1)   is a column vector whose th component is  (1)   .Similarly, for the second term in (21), we have where  2 = [ ℎℎ ;  ℎV ;  VV ],  ℎℎ ,  ℎV and  VV are three second-order difference matrices with size  ×  such that and the weighting matrix Λ  2 is a block matrix defined by here w (2)   is a column vector whose th component is Substituting ( 22) and ( 24) into ( 21), we can rewrite (21) in the more compact form as follows: Since the objective function in optimization (28) is a quadratic function, according to the theory of variational methods, u +1 is the solution to the following Euler-Lagrange equation: System (29) can be solved efficiently using the conjugate gradient (CG) algorithm.In summary, the proposed algorithm for our hybrid TV model ( 10) is composed of four steps.
Otherwise, set  =  + 1 and go to Step 2.
Here, we make one comment about implementation of this algorithm.In Step 3, since the linear system of equations cannot be solved analytically due to its huge dimension ( ×  for an  ×  image), we take as an operator and solve the system by CG algorithm in our implementation.

Experimental Results
In this section, we present some numerical results to illustrate the performance of our method.We compare our results with those obtained by four state-of-the-art methods: (a) TVAL3 [21]; (b) RecPF [22]; (c) Hybrid TVL1 [23]; and (d) second-order TV [17].MATLAB codes of those algorithms are available at the following web sites: http://www.caam.rice.edu/∼optimization/L1/TVAL3/http://www.caam.rice.edu/∼optimization/L1/RecPF/http://www.dssz.com/495330.htmlhttp://www.engineering.uiowa.edu/∼jcb/Software/HDTV/Software.htmlTo compare the performance fairly, all parameters of those algorithms are set as the suggestion values by the authors in [17,[21][22][23], respectively.The reconstruction performance of each method is evaluated in terms of two indicators: the reconstruction relative error (RE) and peak signal to noise ratio (PSNR).The relative error of reconstructed image is defined as follows: where x and x are the original and the reconstructed signal/image, respectively.Apparently, the lower the value of relative error is, the better reconstructed performance will be.The PSNR of reconstructed images is defined by where x and x are the original and the reconstructed image, respectively, and ,  is the size of image.
In our experiments, we consider a common task of reconstructing an image from their undersampled Fourier samples which is an important problem in magnetic resonance imaging (MRI).We generate our test sets using four images (see Figure 1): three natural images and one MR image.In each test, we acquire compressive samples by first applying randomized partial discrete Fourier transform (DFT) encoding and then adding Gaussian white noise with different signal noise ratio (SNR) level (30 dB and 40 dB).The encoding procedure can be formulated as where F represents Fourier transform matrix,  ∈  × is a selection matrix containing  rows of the identity matrix of order , n ∈   is Gaussian white noise, and F serves as the sensing matrix.The reconstructed images of "Lena" image with sampling ratio 25% and noise level 30 dB are shown in Figure 2. We note that the reconstructed image of TVAL3 is contaminated by noise.The results of RecPF and Hybrid TVL1 provide better performance than that of TVAL3 but with some loss of fine details, for example, the texture in hair and hat regions.Comparing the results of those five methods, we observe that second-order TV and our method preserve fine details more effectively than other three methods.From Figures 2(e) and 2(f), we can see that the result of our method have better image quality than that of second-order TV.
In order to highlight the differences, we zoom in the region (marked by the red box in Figure 2(a)) of all the images shown in Figure 2 and present the details in Figure 3.It further demonstrates that second-order TV and our method are superior to other three methods.Furthermore, it can be seen from the regions marked by red boxes and green arrows in Figures 3(e) and 3(f) that second-order TV provide better preservation of smooth image regions and our method provide better preservation of fine details.
Figure 4 compares the reconstructed results of "Brain" image by different methods with sampling ratio 20% and SNR level 40 dB.We can see that RecPF and Hybrid TVL1based reconstruction result in a little patchy artifact and some loss of fine details.But on the whole, all five methods provide very good reconstruction of image features, and it is hard to tell which one is better.To make the difference between those methods more evident, we zoom in the region  marked by the red box in Figure 4(a) of each image and show them in Figure 5.The regions marked by green arrows clearly demonstrate that our method provides consistently improved reconstructions.
To characterize the performance quantitatively, we show the PSNRs of the reconstructed images at various sampling ratios and SNR levels in Table 1.We observe that our method provides the best overall SNR for most of the cases.Our method obtains an improvement of around 0.35 dB on average over second-order TV method for "Lena" image reconstruction, and 0.33 dB, 0.30 dB, and 0.42 dB for "Pepper", "Barbara", and "Brain" images reconstructions, respectively.Moreover, to show the influence of sampling ratio on the performance of reconstruction, we plot relative errors of reconstructed images with various sampling ratio (from 15% to 50%, with interval 5%) in Figure 6.We can see that relative errors of our method keep the lowest level among these methods and our method obtain notable improvements when sampling ratio is less than 35%.
We compare the average central processing unit (CPU) times of the different methods in Table 2.All experiments are performed under Windows 7 and MATLAB V8.0 (R2012b) running on a Lenovo workstation with an Intel (R) Xeon (R) CPU W3520 at 2.67 GHz and 4 GB of memory.
From Table 2, we see that RecPF is faster than other four algorithms.This is primarily because RecPF solves the TV problem using shrinkage and fast Fourier transforms (FFT).However, this technique can also be applied to accelerate our algorithm, which we plan to pursue in the future.

Conclusion
This paper presents a hybrid variational model for image compressive sensing.The model minimizes the sum of three terms corresponding to least squares data fitting, firstorder TV and second-order TV.We propose an efficient majorization-minimization algorithm to determine the solution of our model.We test our method with nature images and MR image.Comparisons of the proposed regularization method with four state-of-the-art algorithms demonstrate the significant improvement in the quality of the reconstructed images.Although achieving better performance in terms of RE and PSNR, our algorithm is slower than FFT-based method such as RecPF.How to accelerate the algorithm is an important problem which we will investigate in forthcoming research.We hope that our method is useful in relevant areas of image compressive sensing such as SAR/ISAR imaging and MR images reconstruction.

Figure 3 :
Figure 3: Zoom into a region of images shown in Figure 2 (the region is marked in Figure 2(a) by the red box).

Figure 5 :
Figure 5: Zoom into a region of images shown in Figure 4 (the region is marked in Figure 4(a) by the red box).

Figure 6 :
Figure 6: Relative errors of various methods versus sampling ratio.

Table 1 :
PSNRs of the reconstructed images.

Table 2 :
Average CPU times of different methods (s).