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

Authors: Can Feng; Liang Xiao; Zhihui Wei
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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).


---

## Body

## 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 imageu∈RN1×N2 is encoded into noisy compressive samples by
(1)y=Φu+n.
Here, Φ∈CM×N is the sensing matrix (M is the number of measurements and N=N1×N2 is the number of total pixels, M≪N), y∈CM 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–3]. One common approach is to formulate the CS decoding as the following optimization problem [4]:
(2)u^=argminu∥y-Φu∥22+λℛ(u),
where ℛ(u) is the regularization term and λ is a tunable regularization parameter that trades off the fidelity term ∥y-Φu∥22 and the regularization term ℛ(u). One of popular choices of the regularization term is the total variation (TV) regularizer defined by the l1 norm of the image gradient [5, 6]:
(3)ℛ1(u)=∫Ω|∇u|dxdy.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 l1 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. Section2 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 ofu 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:
(4)uθ,1(x,y)=uxcosθ+uysinθ=[cosθ,sinθ][∂u∂x∂u∂y].
Similarly, for the second-order directional derivative uθ,2, we have
(5)uθ,2(x,y)=[cos2θ,2cosθsinθ,sin2θ][∂2u∂x2∂2u∂x∂y∂2u∂y2].In order to extend the standard TV scheme to higher-order TV regularizer, we reinterpret the TV regularizer as a group separablel1-l2 mixed norm of directional derivatives of the specified image. To do this, we rewrite gradient magnitude as follows:
(6)|∇u|=|ux|2+|uy|2=1π∫02π|uxcosθ|2dθ+1π∫02π|uysinθ|2dθ=(1π(∫02π|uxcosθ|2dθ+2∫02πuxuycosθsinθdθpppppp+∫02π|uysinθ|2dθ))1/2=1π∫02π|uxcosθ+uysinθ|2dθ=1π∫02π|uθ,1|2dθ.Hence, TV regularizer defined in (3) can be expressed as a functional involving the directional derivatives of u as follows:
(7)ℛ1(u)=∫Ω|∇u|dxdy=∫Ω1π∫02π|uθ,1|2dθdxdy.Based on the above reinterpretation, we introduce second-order TV regularizer by replacing first-order directional derivativeuθ,1 in (7) with second-order directional derivative uθ,2, thus we have
(8)ℛ2(u)=∫Ω1π∫02π|uθ,2|2dθdxdy.Substituting (5) into (8), we can rewrite the second-order TV regularizer ℛ2(u) as follows:
(9)ℛ2(u)=∫Ω(1π∫02π|uxxcos2θ+2cosθsinθuxyppppppppppi+uyysin2θ|2dθ1π)1/2dxdy=∫Ω123uxx2+4uxy2+3uyy2+2uxxuyydxdy.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 samplesy and sensing matrix Φ, the proposed reconstruction model can be formulated as
(10)u^=argminu∥y-Φu∥22+λℛ1(u)+μℛ2(u),
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 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(11)ℛ1(u)=∑i=1N(Dihu)2+(Divu)2,
where Dih and Div are operators corresponding to horizontal and vertical first-order differences at pixel i, respectively, and the summation ∑i=1N‍ takes over all pixels. Similarly, we can write discrete version of second-order TV regularizer as
(12)ℛ2(u)=∑i=1N34(Dihhu)2+(Dihvu)2+34(Divvu)2+12DihhuDivvu.
Here Dihh, Divv, and Dihv are operators corresponding to, respectively, horizontal-horizontal, vertical-vertical, horizontal-vertical second-order differences at pixel i.For the sake of simplicity, we denote(13)η1i(u)≜(Dihu)2+(Divu)2,η2i(u)≜34(Dihhu)2+(Dihvu)2+34(Divvu)2+12DihhuDivvu.
Thus, (10) can be rewritten as
(14)u^=argminu∥y-Φu∥22+λ∑i=1Nη1i(u)+μ∑i=1Nη2i(u).


### 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.LetF(x) be the cost function to be minimized. Instead of minimizing F(x) directly, the MM approaches solve a sequence of surrogate problems {Qk(x),k=0,1,2,…} which are solved easier than F(x). As the result, a sequence of {xk,k=1,2,…} obtained by
(15)xk+1=argminxQk(x)
are then applied to estimate the solution of F(x). In (15), Qk(x) is the majorizer of function F(x) at a fixed point xk which satisfies the following two properties:
(16)(a)Qk(x)≥F(x),∀x,(17)(b)Qk(xk)=F(xk).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:(18)Qk(u)=∑i=1N(λη1i(uk)+μη2i(uk))+∑i=1N(λ2η1i(u)-η1i(uk)η1i(uk)+μ2η2i(u)-η2i(uk)η2i(uk))+∥y-Φu∥22.Using the inequalitya≤a′+(a-a′)/(2a′), we have
(19)Qk(u)≥∥y-Φu∥22+λ∑i=1Nη1i(u)+μ∑i=1Nη2i(u),
for any u, with equality for uk.Equation (19) means that majorizer Qk(u) satisfies the two properties (16) and (17). So, the solution of optimization problem (14) can be obtained by finding the minimizer of Qk(u), k=0,1,2,…. Consider
(20)uk+1=argminuQk(u),Notice that, for fixedk, the terms η1i(uk) and η2i(uk) 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
(21)uk+1=argminuλ2∑i=1Nη1i(u)η1i(uk)+μ2∑i=1Nη2i(u)η2i(uk)+∥y-Φu∥22.LetDh and Dv denote the horizontal and vertical global finite difference matrices with size N×N such that Dhu and Dvu are the vectors of all horizontal and vertical first-order differences of image u. We then have
(22)λ2∑i=1Nη1i(u)η1i(uk)=λ2∑i=1N(Dihu)2+(Divu)2η1i(uk)=(D1hu,…,DNhu,D1vu,…,DNvu)×(w1(1)⋱wN(1)w1(1)⋱wN(1))×(D1hu⋮D1huD1vu⋮D1vu)=(D1u)TΛ1kD1u=uTD1TΛ1kD1u,
where
(23)D1=[DhDv],wi(1)=λ2η1i(uk),i=1,2,…,N,Λ1k=diag{wk(1);wk(1)},
here, wk(1) is a column vector whose ith component is wi(1).Similarly, for the second term in (21), we have
(24)μ2∑i=1Nη2i(u)η1i(uk)=μ2∑i=1N((η1i(uk))-1(34(Dihhu)2+(Dihvu)2+34(Divvu)2ppppppipp+12DihhuDivvu(Dihvu)2)(η1i(uk))-1)=(Dhhu,Dhvu,Dvvu)×14(3diag{wk(2)}0diag{wk(2)}04diag{wk(2)}0diag{wk(2)}03diag{wk(2)})×(DhhuDhvuDvvu)=(D2u)TΛ2kD2u=uTD2TΛ2kD2u,
where D2=[Dhh;Dhv;Dvv], Dhh, Dhv and Dvv are three second-order difference matrices with size N×N such that
(25)Dhhu=(D1hhu⋮D1hhu),Dhvu=(D1hvu⋮D1hvu),Dvvu=(D1vvu⋮D1vvu),
and the weighting matrix Λ2k is a block matrix defined by
(26)Λ2k=14(3diag{wk(2)}0diag{wk(2)}04diag{wk(2)}0diag{wk(2)}03diag{wk(2)}),
here wk(2) is a column vector whose ith component is
(27)wi(2)=μ2η2i(uk),i=1,2,…,N.Substituting (22) and (24) into (21), we can rewrite (21) in the more compact form as follows:
(28)uk+1=argminuuTD1TΛ1kD2u+uTD2TΛ2kD2u+∥y-Φu∥22.Since the objective function in optimization (28) is a quadratic function, according to the theory of variational methods, uk+1 is the solution to the following Euler-Lagrange equation:
(29)(D1TΛ1kD1+D2TΛ2kD2+ΦTΦ)u=ΦTy.
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.
Inputy, Φ, λ, μ, ε  (ε>0 is a tolerance). Set k=0. Initialize u0=ΦTy.Step 2.
(a)
ComputeΛ1k according to (23) for fixed uk.(b)
ComputeΛ2k according to (26) and (27) for fixed uk.Step 3.
Solve Euler-Lagrange equation (29) using CG algorithm, obtain uk+1.Step 4.
Terminate if relative change∥uk+1-uk∥2/∥uk∥2<ε. Otherwise, set k=k+1 and go to Step 2.Here, we make one comment about implementation of this algorithm. In Step3, since the linear system of equations cannot be solved analytically due to its huge dimension (MN×MN for an M×N image), we take D1TΛ1kD1+D2TΛ2kD2+ΦTΦ as an operator and solve the system by CG algorithm in our implementation.


## 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 samplesy and sensing matrix Φ, the proposed reconstruction model can be formulated as
(10)u^=argminu∥y-Φu∥22+λℛ1(u)+μℛ2(u),
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 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(11)ℛ1(u)=∑i=1N(Dihu)2+(Divu)2,
where Dih and Div are operators corresponding to horizontal and vertical first-order differences at pixel i, respectively, and the summation ∑i=1N‍ takes over all pixels. Similarly, we can write discrete version of second-order TV regularizer as
(12)ℛ2(u)=∑i=1N34(Dihhu)2+(Dihvu)2+34(Divvu)2+12DihhuDivvu.
Here Dihh, Divv, and Dihv are operators corresponding to, respectively, horizontal-horizontal, vertical-vertical, horizontal-vertical second-order differences at pixel i.For the sake of simplicity, we denote(13)η1i(u)≜(Dihu)2+(Divu)2,η2i(u)≜34(Dihhu)2+(Dihvu)2+34(Divvu)2+12DihhuDivvu.
Thus, (10) can be rewritten as
(14)u^=argminu∥y-Φu∥22+λ∑i=1Nη1i(u)+μ∑i=1Nη2i(u).


## 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.LetF(x) be the cost function to be minimized. Instead of minimizing F(x) directly, the MM approaches solve a sequence of surrogate problems {Qk(x),k=0,1,2,…} which are solved easier than F(x). As the result, a sequence of {xk,k=1,2,…} obtained by
(15)xk+1=argminxQk(x)
are then applied to estimate the solution of F(x). In (15), Qk(x) is the majorizer of function F(x) at a fixed point xk which satisfies the following two properties:
(16)(a)Qk(x)≥F(x),∀x,(17)(b)Qk(xk)=F(xk).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:(18)Qk(u)=∑i=1N(λη1i(uk)+μη2i(uk))+∑i=1N(λ2η1i(u)-η1i(uk)η1i(uk)+μ2η2i(u)-η2i(uk)η2i(uk))+∥y-Φu∥22.Using the inequalitya≤a′+(a-a′)/(2a′), we have
(19)Qk(u)≥∥y-Φu∥22+λ∑i=1Nη1i(u)+μ∑i=1Nη2i(u),
for any u, with equality for uk.Equation (19) means that majorizer Qk(u) satisfies the two properties (16) and (17). So, the solution of optimization problem (14) can be obtained by finding the minimizer of Qk(u), k=0,1,2,…. Consider
(20)uk+1=argminuQk(u),Notice that, for fixedk, the terms η1i(uk) and η2i(uk) 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
(21)uk+1=argminuλ2∑i=1Nη1i(u)η1i(uk)+μ2∑i=1Nη2i(u)η2i(uk)+∥y-Φu∥22.LetDh and Dv denote the horizontal and vertical global finite difference matrices with size N×N such that Dhu and Dvu are the vectors of all horizontal and vertical first-order differences of image u. We then have
(22)λ2∑i=1Nη1i(u)η1i(uk)=λ2∑i=1N(Dihu)2+(Divu)2η1i(uk)=(D1hu,…,DNhu,D1vu,…,DNvu)×(w1(1)⋱wN(1)w1(1)⋱wN(1))×(D1hu⋮D1huD1vu⋮D1vu)=(D1u)TΛ1kD1u=uTD1TΛ1kD1u,
where
(23)D1=[DhDv],wi(1)=λ2η1i(uk),i=1,2,…,N,Λ1k=diag{wk(1);wk(1)},
here, wk(1) is a column vector whose ith component is wi(1).Similarly, for the second term in (21), we have
(24)μ2∑i=1Nη2i(u)η1i(uk)=μ2∑i=1N((η1i(uk))-1(34(Dihhu)2+(Dihvu)2+34(Divvu)2ppppppipp+12DihhuDivvu(Dihvu)2)(η1i(uk))-1)=(Dhhu,Dhvu,Dvvu)×14(3diag{wk(2)}0diag{wk(2)}04diag{wk(2)}0diag{wk(2)}03diag{wk(2)})×(DhhuDhvuDvvu)=(D2u)TΛ2kD2u=uTD2TΛ2kD2u,
where D2=[Dhh;Dhv;Dvv], Dhh, Dhv and Dvv are three second-order difference matrices with size N×N such that
(25)Dhhu=(D1hhu⋮D1hhu),Dhvu=(D1hvu⋮D1hvu),Dvvu=(D1vvu⋮D1vvu),
and the weighting matrix Λ2k is a block matrix defined by
(26)Λ2k=14(3diag{wk(2)}0diag{wk(2)}04diag{wk(2)}0diag{wk(2)}03diag{wk(2)}),
here wk(2) is a column vector whose ith component is
(27)wi(2)=μ2η2i(uk),i=1,2,…,N.Substituting (22) and (24) into (21), we can rewrite (21) in the more compact form as follows:
(28)uk+1=argminuuTD1TΛ1kD2u+uTD2TΛ2kD2u+∥y-Φu∥22.Since the objective function in optimization (28) is a quadratic function, according to the theory of variational methods, uk+1 is the solution to the following Euler-Lagrange equation:
(29)(D1TΛ1kD1+D2TΛ2kD2+ΦTΦ)u=ΦTy.
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.
Inputy, Φ, λ, μ, ε  (ε>0 is a tolerance). Set k=0. Initialize u0=ΦTy.Step 2.
(a)
ComputeΛ1k according to (23) for fixed uk.(b)
ComputeΛ2k according to (26) and (27) for fixed uk.Step 3.
Solve Euler-Lagrange equation (29) using CG algorithm, obtain uk+1.Step 4.
Terminate if relative change∥uk+1-uk∥2/∥uk∥2<ε. Otherwise, set k=k+1 and go to Step 2.Here, we make one comment about implementation of this algorithm. In Step3, since the linear system of equations cannot be solved analytically due to its huge dimension (MN×MN for an M×N image), we take D1TΛ1kD1+D2TΛ2kD2+ΦTΦ 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.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–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:
(30)RE(x)=∥x-x-∥2∥x-∥2,
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
(31)PSNR(x)=10log10(m×n×2552∥x-x-∥2),
where x- and x are the original and the reconstructed image, respectively, and m, n 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 Figure1): 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)y=Pℱu+n,
where ℱ represents Fourier transform matrix, P∈RM×N is a selection matrix containing p rows of the identity matrix of order N, n∈CN is Gaussian white noise, and Pℱ 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 Figure2. 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.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.
(a)
(b)
(c)
(d)
(e)
(f)In order to highlight the differences, we zoom in the region (marked by the red box in Figure2(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.Zoom into a region of images shown in Figure2 (the region is marked in Figure 2(a) by the red box).
(a)
(b)
(c)
(d)
(e)
(f)Figure4 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.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)Zoom into a region of images shown in Figure4 (the region is marked in Figure 4(a) by the red box).
(a)
(b)
(c)
(d)
(e)
(f)To characterize the performance quantitatively, we show the PSNRs of the reconstructed images at various sampling ratios and SNR levels in Table1. 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.
Sampling ratio (%)
20
25
33
50
Noise level (dB)
30
40
30
40
30
40
30
40
Lena (256 × 256)
TVAL3
26.99
27.14
28.12
27.65
29.03
29.19
31.51
32.34
RecPF
29.93
30.92
30.76
32.06
32.01
34.06
34.03
37.94
Hybrid TVL1
30.12
31.13
30.97
32.26
32.24
34.43
34.26
38.41
Second-order TV
30.57
31.52
31.89
33.87
32.83
36.09
34.74
40.08
Our method
31.05
32.74
32.93
34.12
32.99
36.41
34.43
40.24
Pepper (256 × 256)
TVAL3
30.08
30.41
30.75
32.23
32.78
33.93
35.08
38.08
RecPF
30.92
31.95
31.72
33.79
33.01
35.68
34.91
39.25
Hybrid TVL1
31.10
32.29
31.93
34.10
33.24
36.05
35.06
39.65
Second-order TV
31.72
33.76
31.72
35.79
33.67
37.67
35.07
41.07
Our method
32.25
34.24
32.66
36.18
33.77
37.96
34.82
41.21
Barbara (256 × 256)
TVAL3
27.25
27.57
27.78
28.27
29.14
29.45
31.42
32.45
RecPF
28.70
29.29
29.45
30.28
30.67
31.80
33.07
35.62
Hybrid TVL1
28.86
29.50
29.63
30.55
30.93
32.20
33.36
36.20
Second-order TV
29.25
30.21
30.02
31.54
31.41
33.44
33.86
37.98
Our method
29.62
30.69
30.46
31.86
31.69
33.76
33.89
38.13
Brain (256 × 256)
TVAL3
32.63
33.22
33.76
34.49
35.23
36.32
37.52
40.18
RecPF
31.49
32.11
32.51
33.50
33.89
35.27
36.31
39.12
Hybrid TVL1
31.69
32.38
32.75
33.83
34.23
35.75
36.63
39.71
Second-order TV
32.65
33.64
33.66
35.13
35.20
37.29
37.52
41.11
Our method
33.19
34.72
34.02
35.63
35.46
37.70
37.35
41.55Relative errors of various methods versus sampling ratio.
(a)
Lena, Noise level: 40 dB
(b)
Pepper, Noise level: 40 dB
(c)
Barbara, Noise level: 40 dB
(d)
Brain, Noise level: 40 dBWe compare the average central processing unit (CPU) times of the different methods in Table2. 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).
Sampling ratio (%)
20
25
0.33
0.5
Noise level (dB)
30
40
30
40
30
40
30
40
TVAL3
2.21
2.15
2.12
1.90
1.82
1.81
1.81
1.76
RecPF
0.80
0.68
0.69
0.67
0.69
0.67
0.67
0.58
Hybrid TVL1
1.76
1.72
1.68
1.67
1.56
1.56
1.51
1.50
Second-order TV
81.25
78.36
77.18
75.63
75.02
75.01
74.45
74.56
Our method
29.50
26.79
26.37
24.13
24.01
23.76
23.24
23.12From Table2, 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.



---

*Source: 470165-2013-10-23.xml*