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

Abstract

A hybrid variational model combined first- and second-order total variation for image reconstruction from its finite number of noisy compressive samples is proposed in this paper. Inspired by majorization-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).

1. 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., analog-digital 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 𝐮∈𝑅^𝑁1×𝑁2 is encoded into noisy compressive samples by 𝐲=Φ𝐮+𝐧. (1) Here, Φ∈𝐶^𝑀×𝑁 is the sensing matrix (𝑀 is the number of measurements and 𝑁=𝑁₁×𝑁₂ is the number of total pixels, 𝑀≪𝑁), 𝐲∈𝐶^𝑀 is the compressive samples of original image 𝐮, and 𝐧 is Gaussian white noise with standard deviation 𝜎. While in the decoding phase, a nonlinear recovery algorithm is used to reconstruct the original image 𝐮 from compressive samples 𝐲 [1–3]. One common approach is to formulate the CS decoding as the following optimization problem [4]: 𝐮=argmin_𝐮‖𝐲−Φ𝐮‖²₂+𝜆ℛ(𝐮), (2) where ℛ(𝐮) is the regularization term and 𝜆 is a tunable regularization parameter that trades off the fidelity term ‖𝐲−Φ𝐮‖²₂ and the regularization term ℛ(𝐮). One of popular choices of the regularization term is the total variation (TV) regularizer defined by the 𝑙₁ norm of the image gradient [5, 6]: ℛ₁(𝐮)=_Ω|∇𝐮|𝑑𝑥𝑑𝑦. (3)

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–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 𝑙₁ 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–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.

2. Second-Order TV Penalty

In this section, we reinterpret classic TV regularizer and extend it to second-order TV regularizer. We denote first-order directional derivative of 𝐮 along the unit vector 𝐞_𝜃=(cos𝜃,sin𝜃) as 𝐮_𝜃,1, which can be expressed as the linear combination of horizontal and vertical directional derivatives as follows: 𝐮_𝜃,1𝑥,𝑦=𝐮_𝑥cos𝜃+𝐮_𝑦sin𝜃=[cos𝜃,sin𝜃]𝜕𝐮𝜕𝑥𝜕𝐮𝜕𝑦. (4) Similarly, for the second-order directional derivative 𝐮_𝜃,2, we have 𝐮_𝜃,2𝑥,𝑦=cos²𝜃,2cos𝜃sin𝜃,sin²𝜃𝜕²𝐮𝜕𝑥²𝜕²𝐮𝜕𝜕𝑥𝜕𝑦₂𝐮𝜕𝑦². (5)

In order to extend the standard TV scheme to higher-order TV regularizer, we reinterpret the TV regularizer as a group separable 𝑙₁−𝑙₂ mixed norm of directional derivatives of the specified image. To do this, we rewrite gradient magnitude as follows: |∇𝐮|=𝐮_𝑥²+𝐮_𝑦²=1𝜋₀^2𝜋𝐮_𝑥cos𝜃²1𝑑𝜃+𝜋₀^2𝜋𝐮_𝑦sin𝜃²=1𝑑𝜃𝜋₀^2𝜋𝐮_𝑥cos𝜃²𝑑𝜃+2₀^2𝜋𝐮_𝑥𝐮_𝑦+cos𝜃sin𝜃𝑑𝜃_02𝜋𝐮_𝑦sin𝜃²𝑑𝜃^1/2=1𝜋₀^2𝜋𝐮_𝑥cos𝜃+𝐮_𝑦sin𝜃²=𝑑𝜃1𝜋₀^2𝜋𝐮_𝜃,1²𝑑𝜃. (6)

Hence, TV regularizer defined in (3) can be expressed as a functional involving the directional derivatives of 𝐮 as follows: ℛ₁(𝐮)=_Ω|∇𝐮|𝑑𝑥𝑑𝑦=_Ω1𝜋₀^2𝜋𝐮_𝜃,1²𝑑𝜃𝑑𝑥𝑑𝑦. (7)

Based on the above reinterpretation, we introduce second-order TV regularizer by replacing first-order directional derivative 𝐮_𝜃,1 in (7) with second-order directional derivative 𝐮_𝜃,2, thus we have ℛ₂(𝐮)=_Ω1𝜋₀^2𝜋𝐮_𝜃,2²𝑑𝜃𝑑𝑥𝑑𝑦. (8)

Substituting (5) into (8), we can rewrite the second-order TV regularizer ℛ₂(𝐮) as follows: ℛ₂(𝐮)=_Ω1𝜋₀^2𝜋𝐮_𝑥𝑥cos²𝜃+2cos𝜃sin𝜃𝐮_𝑥𝑦+𝐮_𝑦𝑦sin²𝜃²𝑑𝜃^1/2=𝑑𝑥𝑑𝑦_Ω123𝐮²_𝑥𝑥+4𝐮²_𝑥𝑦+3𝐮²_𝑦𝑦+2𝐮_𝑥𝑥𝐮_𝑦𝑦𝑑𝑥𝑑𝑦. (9)

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 second-order derivatives are steerable which enables us to obtain analytical expressions for the new regularizer as the standard TV regularizer. Furthermore, the minimization of second-order 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.

3. The Proposed Variational Approach for Image Compressive Sensing

3.1. 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 𝐲 and sensing matrix Φ, the proposed reconstruction model can be formulated as 𝐮=argmin_𝐮‖𝐲−Φ𝐮‖²₂+𝜆ℛ₁(𝐮)+𝜇ℛ₂(𝐮), (10) 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.

3.2. 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 𝐮 has been uniformly sampled. At this point, we switch the notation and start assuming that 𝐮 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 ℛ₁(𝐮)=^𝑁_𝑖=1𝐷^ℎ_𝑖𝐮²+𝐷^𝑣_𝑖𝐮², (11) where 𝐷^ℎ_𝑖 and 𝐷^𝑣_𝑖 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 ℛ₂(𝐮)=^𝑁_𝑖=134𝐷_𝑖^ℎℎ𝐮²+𝐷_𝑖^ℎ𝑣𝐮²+34𝐷_𝑖^𝑣𝑣𝐮²+12𝐷_𝑖^ℎℎ𝐮𝐷_𝑖^𝑣𝑣𝐮. (12) Here 𝐷_𝑖^ℎℎ, 𝐷_𝑖^𝑣𝑣, and 𝐷_𝑖^ℎ𝑣 are operators corresponding to, respectively, horizontal-horizontal, vertical-vertical, horizontal-vertical second-order differences at pixel 𝑖.

For the sake of simplicity, we denote 𝜂^𝑖₁(𝐮)≜𝐷^ℎ_𝑖𝐮²+𝐷^𝑣_𝑖𝐮²,𝜂^𝑖₂(𝐮)≜34𝐷_𝑖^ℎℎ𝐮²+𝐷_𝑖^ℎ𝑣𝐮²+34𝐷_𝑖^𝑣𝑣𝐮²+12𝐷_𝑖^ℎℎ𝐮𝐷_𝑖^𝑣𝑣𝐮. (13) Thus, (10) can be rewritten as 𝐮=argmin_𝐮‖𝐲−Φ𝐮‖²₂+𝜆^𝑁_𝑖=1𝜂^𝑖₁(𝐮)+𝜇^𝑁_𝑖=1𝜂^𝑖₂(𝐮). (14)

3.3. 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 𝐹(𝐱) be the cost function to be minimized. Instead of minimizing 𝐹(𝐱) directly, the MM approaches solve a sequence of surrogate problems {𝑄_𝑘(𝐱),𝑘=0,1,2,…} which are solved easier than 𝐹(𝐱). As the result, a sequence of {𝐱_𝑘,𝑘=1,2,…} obtained by 𝐱_𝑘+1=argmin_𝐱𝑄_𝑘(𝐱) (15) are then applied to estimate the solution of 𝐹(𝐱). In (15), 𝑄_𝑘(𝐱) is the majorizer of function 𝐹(𝐱) at a fixed point 𝐱_𝑘 which satisfies the following two properties: (a)𝑄_𝑘(𝐱)(𝐱)(b)𝑄≥𝐹,∀𝐱,_𝑘𝐱_𝑘=𝐹𝐱_𝑘. (16)(17)

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: 𝑄_𝑘(𝐮)=^𝑁_𝑖=1𝜆𝜂^𝑖₁𝐮_𝑘+𝜇𝜂^𝑖₂𝐮_𝑘+^𝑁_𝑖=1𝜆2𝜂^𝑖₁(𝐮)−𝜂^𝑖₁𝐮_𝑘𝜂^𝑖₁𝐮_𝑘+𝜇2𝜂^𝑖₂(𝐮)−𝜂^𝑖₂𝐮_𝑘𝜂^𝑖₂𝐮_𝑘+‖𝐲−Φ𝐮‖²₂. (18)

Using the inequality √𝑎≤𝑎^′+(𝑎−𝑎^′)/(2𝑎^′), we have 𝑄_𝑘(𝐮)≥‖𝐲−Φ𝐮‖²₂+𝜆^𝑁_𝑖=1𝜂^𝑖₁(𝐮)+𝜇^𝑁_𝑖=1𝜂^𝑖₂(𝐮), (19) for any 𝐮, with equality for 𝐮_𝑘.

Equation (19) means that majorizer 𝑄_𝑘(𝐮) satisfies the two properties (16) and (17). So, the solution of optimization problem (14) can be obtained by finding the minimizer of 𝑄_𝑘(𝐮), 𝑘=0,1,2,…. Consider 𝐮_𝑘+1=argmin_𝐮𝑄_𝑘(𝐮), (20)

Notice that, for fixed 𝑘, the terms 𝜂^𝑖₁(𝐮_𝑘) and 𝜂^𝑖₂(𝐮_𝑘) 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 𝐮_𝑘+1=argmin_𝐮𝜆2^𝑁_𝑖=1𝜂^𝑖₁(𝐮)𝜂^𝑖₁𝐮_𝑘+𝜇2^𝑁_𝑖=1𝜂^𝑖₂(𝐮)𝜂^𝑖₂𝐮_𝑘+‖𝐲−Φ𝐮‖²₂. (21)

Let 𝐷^ℎ and 𝐷^𝑣 denote the horizontal and vertical global finite difference matrices with size 𝑁×𝑁 such that 𝐷^ℎ𝐮 and 𝐷^𝑣𝐮 are the vectors of all horizontal and vertical first-order differences of image 𝐮. We then have 𝜆2^𝑁_𝑖=1𝜂^𝑖₁(𝐮)𝜂^𝑖₁𝐮_𝑘=𝜆2^𝑁_𝑖=1𝐷^ℎ_𝑖𝐮²+𝐷^𝑣_𝑖𝐮²𝜂^𝑖₁𝐮_𝑘=𝐷^ℎ₁𝐮,…,𝐷^ℎ_𝑁𝐮,𝐷^𝑣₁𝐮,…,𝐷^𝑣_𝑁𝐮×𝑤₁⁽¹⁾⋱𝑤_𝑁⁽¹⁾𝑤₁⁽¹⁾⋱𝑤_𝑁⁽¹⁾×𝐷^ℎ₁𝐮⋮𝐷^ℎ₁𝐮𝐷^𝑣₁𝐮⋮𝐷^𝑣₁𝐮=𝐷₁𝐮^𝑇Λ^𝑘₁𝐷₁𝐮=𝐮^𝑇𝐷^𝑇₁Λ^𝑘₁𝐷₁𝐮, (22) where 𝐷₁=𝐷^ℎ𝐷^𝑣,𝑤_𝑖⁽¹⁾=𝜆2𝜂^𝑖₁𝐮_𝑘Λ,𝑖=1,2,…,𝑁,_𝑘1=diag𝐰_𝑘⁽¹⁾;𝐰_𝑘⁽¹⁾, (23) here, 𝐰_𝑘⁽¹⁾ is a column vector whose 𝑖th component is 𝑤_𝑖⁽¹⁾.

Similarly, for the second term in (21), we have 𝜇2^𝑁_𝑖=1𝜂^𝑖₂(𝐮)𝜂^𝑖₁𝐮_𝑘=𝜇2^𝑁_𝑖=134𝐷_𝑖^ℎℎ𝐮²+𝐷_𝑖^ℎ𝑣𝐮²+34𝐷_𝑖^𝑣𝑣𝐮²+12𝐷_𝑖^ℎℎ𝐮𝐷_𝑖^𝑣𝑣𝐮𝜂^𝑖₁𝐮_𝑘⁻¹=𝐷^ℎℎ𝐮,𝐷^ℎ𝑣𝐮,𝐷^𝑣𝑣𝐮×143diag𝐰_𝑘⁽²⁾0diag𝐰_𝑘⁽²⁾04diag𝐰_𝑘⁽²⁾0diag𝐰_𝑘⁽²⁾03diag𝐰_𝑘⁽²⁾×𝐷^ℎℎ𝐮𝐷^ℎ𝑣𝐮𝐷^𝑣𝑣𝐮=𝐷₂𝐮^𝑇Λ^𝑘₂𝐷₂𝐮=𝐮^𝑇𝐷^𝑇₂Λ^𝑘₂𝐷₂𝐮, (24) where 𝐷₂=[𝐷^ℎℎ;𝐷^ℎ𝑣;𝐷^𝑣𝑣], 𝐷^ℎℎ, 𝐷^ℎ𝑣 and 𝐷^𝑣𝑣 are three second-order difference matrices with size 𝑁×𝑁 such that 𝐷^ℎℎ𝐷𝐮=^1ℎℎ𝐮⋮𝐷₁^ℎℎ𝐮,𝐷^ℎ𝑣𝐷𝐮=^1ℎ𝑣𝐮⋮𝐷₁^ℎ𝑣𝐮,𝐷^𝑣𝑣𝐷𝐮=^1𝑣𝑣𝐮⋮𝐷₁^𝑣𝑣𝐮, (25) and the weighting matrix Λ^𝑘₂ is a block matrix defined by Λ^𝑘₂=143diag𝐰_𝑘⁽²⁾0diag𝐰_𝑘⁽²⁾04diag𝐰_𝑘⁽²⁾0diag𝐰_𝑘⁽²⁾03diag𝐰_𝑘⁽²⁾, (26) here 𝐰_𝑘⁽²⁾ is a column vector whose 𝑖th component is 𝑤_𝑖⁽²⁾=𝜇2𝜂^𝑖₂𝐮_𝑘,𝑖=1,2,…,𝑁. (27)

Substituting (22) and (24) into (21), we can rewrite (21) in the more compact form as follows: 𝐮_𝑘+1=argmin_𝐮𝐮^𝑇𝐷^𝑇₁Λ^𝑘₁𝐷₂𝐮+𝐮^𝑇𝐷^𝑇₂Λ^𝑘₂𝐷₂𝐮+‖𝐲−Φ𝐮‖²₂. (28)

Since the objective function in optimization (28) is a quadratic function, according to the theory of variational methods, 𝐮_𝑘+1 is the solution to the following Euler-Lagrange equation: 𝐷^𝑇₁Λ^𝑘₁𝐷₁+𝐷^𝑇₂Λ^𝑘₂𝐷₂+Φ^𝑇Φ𝐮=Φ^𝑇𝑦. (29) 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.

Step 1. Input 𝑦, Φ, 𝜆, 𝜇, 𝜀  (𝜀>0 is a tolerance). Set 𝑘=0. Initialize 𝑢₀=Φ^𝑇𝑦.

Step 2.  (a)Compute Λ^𝑘₁ according to (23) for fixed 𝐮_𝑘.(b)Compute Λ^𝑘₂ according to (26) and (27) for fixed 𝐮_𝑘.

Step 3. Solve Euler-Lagrange equation (29) using CG algorithm, obtain 𝐮_𝑘+1.

Step 4. Terminate if relative change ‖𝐮_𝑘+1−𝐮_𝑘‖₂/‖𝐮_𝑘‖₂<𝜀. 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.

4. 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.html http://www.engineering.uiowa.edu/~jcb/Software/HDTV/Software.html

To compare the performance fairly, all parameters of those algorithms are set as the suggestion values by the authors in [17, 21–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: (𝐱)=‖RE𝐱−𝐱‖₂‖𝐱‖₂, (30) where 𝐱 and 𝐱 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 PSNR(𝐱)=10log₁₀𝑚×𝑛×255²‖𝐱−𝐱‖₂, (31) where 𝐱 and 𝐱 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 𝐲=𝑃ℱ𝐮+𝐧, (32) where ℱ represents Fourier transform matrix, 𝑃∈𝑅^𝑀×𝑁 is a selection matrix containing 𝑝 rows of the identity matrix of order 𝑁, 𝐧∈𝐶^𝑁 is Gaussian white noise, and 𝑃ℱ serves as the sensing matrix.

Figure 1 

Test images.

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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Reconstructed results of “Lena” image from noisy and undersampled Fourier data (noise level: 30 dB, sampling ratio: 25%). (a) Actual image; (b) TVAL3; (c) RecPF; (d) hybrid TVL1; (e) second-order TV; (f) our method.

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.

(a)

(b)

(c)

(d)

(e)

(f)

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 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 TVL1-based 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Reconstructed results of “Brain” image from noisy and undersampled Fourier data (noise level: 40 dB, sampling ratio: 20%). (a) Actual image; (b) TVAL3; (c) RecPF; (d) hybrid TVL1; (e) second-order TV; (f) our method.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

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

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%.

Table 1 

PSNRs of the reconstructed images.

(a) Lena, Noise level: 40 dB

(b) Pepper, Noise level: 40 dB

(c) Barbara, Noise level: 40 dB

(d) Brain, Noise level: 40 dB

Figure 6 

Relative errors of various methods versus sampling ratio.

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.

Table 2 

Average CPU times of different methods (s).

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.

5. 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, first-order 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the anonymous reviewers for their careful reading and valuable comments on their paper, which helped them to improve that paper. This research is supported by the National Nature Science Foundation of China (61171165, 60802039), the Nature Science Foundation of Jiangsu Province (BK2010488), and the Qinglan Outstanding Scholar Project of Jiangsu Province.