Fast Second Degree Total Variation Method for Image Compressive Sensing

This paper presents a computationally efficient algorithm for image compressive sensing reconstruction using a second degree total variation (HDTV2) regularization. Firstly, a preferably equivalent formulation of the HDTV2 functional is derived, which can be formulated as a weighted L 1-L 2 mixed norm of second degree image derivatives under the spectral decomposition framework. Secondly, using the equivalent formulation of HDTV2, we introduce an efficient forward-backward splitting (FBS) scheme to solve the HDTV2-based image reconstruction model. Furthermore, from the averaged non-expansive operator point of view, we make a detailed analysis on the convergence of the proposed FBS algorithm. Experiments on medical images demonstrate that the proposed method outperforms several fast algorithms of the TV and HDTV2 reconstruction models in terms of peak signal to noise ratio (PSNR), structural similarity index (SSIM) and convergence speed.


Introduction
Recently, the compressive sensing (CS) technique which condenses the information in sparse or compressible images into a small amount of data and yet reconstructs them accurately has been developed for data acquisition in many practical applications, including radar receivers, magnetic resonance imaging (MRI) and microscopy [1,2].
In this paper, we focus on the issue of reconstructing a desired image f:O!< from its noisy and undersampled Fourier measurement g, where O&< 2 is the spatial support of the image. In general, the observed measurement g is often linearly modeled as where A is a linear operator representing the randomized partial Fourier sampling pattern, and e is often assumed to be Gaussian white noise with standard deviation σ. As is well known, the reconstruction of f from its observed measurement g is an ill-posed problem. Variational regularization approach is usually employed to cope with the ill-posedness by formulating image CS reconstruction as an optimization problem: where the first term is called the data fidelity term, while the second one J(f) is called the regularization term which enforces some prior knowledge on the desired image, and the regularization parameter λ!0 balances the contributions of the two terms.
As an important and valid sparsity regularizer, the total variation (TV) regularizer defined by the L 1 norm of image gradient has been applied to image CS reconstruction problems and turned out to be very efficient for preserving image edges [3][4][5]. However, TV favors piecewise constant solutions so that TV-based CS reconstruction methods easily tend to produce severe staircase effect in the smooth regions of reconstructed image. To reduce the staircase effect and preserve image edges sharpness as far as possible, higher order regularization schemes have been gradually investigated by more and more researchers for linear inverse problems [6][7][8][9][10][11]. Recently, a novel family of higher order regularizers termed as higher degree total variation (HDTV) has been proposed for image CS reconstruction, and their experimental results demonstrate that HDTV can effectively minimize the staircase effect while better preserving the edges and singularities [12]. Specifically, HDTV preserves strong directional derivatives at a specified direction but attenuates the small ones at other directions thus encouraging smoothing along line-like features and enhancing them. Moreover, the use of higher degree derivatives yields piecewise linear solutions so that they can minimize the staircase effect and provide state-of-the-art results.
In spite of this, the main challenge of the current HDTV method consists in its high computational complexity compared to the popular TV methods. Due to the non-differentiability of HDTV, the corresponding HDTV-based reconstruction model proposed in [12] is solved by using a majorization-minimization (MM) algorithm, which is similar to the iteratively reweighted algorithms used in the TV minimization problems [13]. Under the MM framework, the MM algorithm proceeds by successively minimizing a sequence of quadratic surrogate functionals that majorize the HDTV regularizer and the data fidelity term, then the minimization of quadratic surrogate functionals is solved by using a conjugate gradient (CG) algorithm. Specifically, the iteratively reweighted MM algorithm alternates between the CG optimization and the update of the weights from current iterate. However, the spatially varying weights often tend to be very large values, then the minimization subproblems of quadratic surrogate functionals are with large condition numbers. As a result, the iteratively reweighted MM algorithm is computational expensive so as to result in slow convergence speed. To reduce the computational complexity and accelerate the convergence speed, we are motivated to investigate more efficient HDTV-based algorithms.
In this paper, we only consider the special case of second degree total variation (HDTV2)based reconstruction model and focus on introducing a computationally efficient HDTV2 reconstruction algorithm. Compared to the previous work [12], our main contributions are listed as follows: 1. Under the spectral decomposition framework, we reinterpret the original HDTV2 regularizer which is the L 1 -L 2 mixed norm of the second degree directional derivatives as a novel weighted L 1 -L 2 mixed norm of the second degree image derivatives.
3. A complete convergence proof of the proposed FBS algorithm is shown relying on the properties of averaged and firmly non-expansive operators. The experimental results of our proposed scheme compared with the recently iteratively reweighted MM scheme [12] demonstrate the significant performance improvements in terms of peak signal to noise ratio (PSNR), structural similarity index (SSIM) and convergence speed.
The rest of this paper is outlined as follows. In next section, we first review the HDTV2-based image reconstruction model and analyze in detail the shortcomings of the iteratively reweighted MM scheme for HDTV2-based reconstruction model. To overcome these shortcomings, a novel equivalent formulation of HDTV2 is derived. Owing to the equivalent formulation, an efficient algorithm using FBS scheme is introduced for HDTV2-based minimization problem and the convergence of the proposed FBS algorithm is also proved. Section Experiments shows the experimental results for medical image CS reconstruction problems. Finally, we conclude the paper in Section Conclusions.

Second Degree Total Variation (HDTV2)-Based Image Reconstruction Model
In this paper, we just consider the special case of HDTV corresponding to the second degree total variation (HDTV2) for image reconstruction.

Review and analysis of HDTV2-based reconstruction model
As detailed in [12], the HDTV2-based reconstruction model is specified aŝ where HDTV 2 (f) is the HDTV2 regularizer defined as the L 1 -L 2 mixed norm of second degree directional derivatives of image f, which is formulated as with f y;2 ðx; yÞ ¼ ½cos 2 y; 2cosysiny; sin 2 y |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} denoted the second degree directional derivative of image f at coordinate (x, y) along the unit vector u θ = [cosθ,sinθ] T , where G 2 (x, y) is a vector consisting of the second degree partial derivatives of image f at coordinate (x, y), and (Á) T denotes the transpose operator. Based on the definitions (4) and (5), the HDTV2 regularizer can be rewritten as where

5.
Since HDTV2 is nondifferentiable, the authors in [12] use an iteratively reweighted MM algorithm to solve model (3) at k th iterate as follows: which then can be solved by using a conjugate gradient (CG) algorithm. Specifically, the MM algorithm alternates between CG optimization (computing f (k+1) subproblem) and the recomputation of the weights from the current iterate (computing E (k) ). However, the spatially varying weights in E (k) often tend to be very large values so that the f (k+1) subproblem is with large condition number. Hence, the resulting MM algorithm is also computationally expensive so as to show slow convergence speed.
To reduce the computational cost and accelerate the convergence speed, we introduce a more efficient HDTV2-based reconstruction algorithm. To this end, we first derive an equivalent formulation of HDTV2, and then use the equivalent formulation to design a computationally efficient HDTV2 reconstruction algorithm in the following sections.

Reinterpretation of HDTV2
From now on, we will reinterpret the HDTV2 regularizer defined in Eq (4) as a weighted L 1 -L 2 mixed norm of second degree image derivatives under the spectral decomposition framework.
As shown in Eq (6), since C 2 is symmetric and positive definite, we can use the spectral decomposition theorem to rewrite C 2 as and we further obtain that 5 which is also symmetric and positive definite. Based on that, we will derive an equivalent formulation of the HDTV2 regularizer in the following proposition. 5 and the second degree vectorial differential operator @ 2 = (@ xx , @ xy , @ yy ) T , then the HDTV2 can be formulated as where || Á || 2 stands for the Euclidean norm. Proof: Combining Eqs (6) and (8), then HDTV2 can be simplified as As described in Eq (10), the HDTV2 regularizer can be reinterpreted as weighted L 1 -L 2 mixed norm of second degree image derivatives, while W 2 can be understood as the weighting matrix. According to the equivalent formulation of HDTV2, it is more preferable for the algorithm we will propose in next section.

Forward-Backward (FB) Splitting Algorithm for HDTV2-Based Image Reconstruction
In this section, we focus on using the equivalent formulation of HDTV2 to design a more efficient HDTV2 reconstruction algorithm.

Discrete Model Formulation
First of all, we will give the discrete formulation of model (3). To simplify our analysis, an image with size r×c is stacked in a vector of size N = r×c, we assume the reflexive boundary condition for images and use the forward finite differences to approximate the second degree derivatives [14].
Next, we define operator V 2 :< N !< 3×N to refer to the discrete version of the vectorial differential operator @ 2 defined in proposition 1, and its adjoint operator is defined as Hereafter, we denote X 2 = < 3×N . Then, let p = (p 1 ,p 2 ,Á Á Á,p N )2X 2 and q = (q 1 ,q 2 ,Á Á Á,q N )2X 2 , we define the inner product h Á ; Á i X 2 and norm jj Á jj X 2 in X 2 as Moreover, let us define operator Λ 2 :X 2 !X 2 and its adjoint operator L Τ 2 : X 2 ! X 2 , that act on vector field p = (p 1 ,p 2 ,Á Á Á,p N )2X 2 as (Λ 2 p) i = W 2 p i and ðL Τ 2 pÞ i ¼ W Τ 2 p i , where (Á) i denotes the i th column element of the corresponding vector filed.
With the l 1 −l 2 mixed norm of a vector field p = (p 1 ,p 2 ,Á Á Á,p N )2X 2 defined as jjp i jj 2 , then the HDTV2 described in Eq (10) can be discretized as In this paper, we concentrate on the issue of image reconstruction from undersampled Fourier data. Therefore, we can write the discrete formulation of the model (3) aŝ where A = PF2C M×N is the sensing matrix with F2C N×N denoted the two-dimensional discrete Fourier transform (DFT) matrix and P2C M×N denoted the selection matrix containing M rows of the identity matrix of order N, g2C M is the observed undersampled Fourier measurement and f2< N is the image to be reconstructed. In the following, we will describe in detail the minimization algorithm for solving above model (12).

The Proposed Algorithm
Let us denote R 1 ðf Þ ¼ 1 2 jjg À Af jj 2 2 and R 2 (f) = λ||Λ 2 V 2 f|| 1,2 , it is clear that R 1 and R 2 are both proper and lower semicontinuous convex functions from < N to [0,+1), then the model (12) becomes the minimization of functions of the form R = R 1 +R 1 . Specifically, let G = {f2< N | 02@R(f)} be the set of solutions of model (12), where @R is the subdifferential of function R, then f2G if and only if where prox βR :< N !< N is proximity operator of function βR, which is defined as where β is known as the proximal step size. Thus, the proximal-type algorithm for model (12) can be recursively constructed as However, the main difficulty with above proximal-type algorithm is that prox βR is hard to compute. Under the operator splitting framework, splitting methods for model (12) do not attempt to compute the proximal operator prox βR of the combined function βR, but instead perform a sequence of calculations involving separately the individual proximal operators prox bR 1 and prox bR 2 which are relatively easier to compute. As detailed in model (12), we note that the gradient of R 1 (f) is Lipschitz continuous with a constant upper bound by ||P|| 2 ||F|| 2 . Therefore, we employ the popular forward-backward (FB) splitting algorithm [15][16][17][18] to solve model (12), then the recursive form of FB splitting becomes where rR 1 (f (k) ) = A T (Af (k) −g) is the gradient of function R 1 at f (k) and the proximal step size τ2(0,2/(||P|| 2 ||F|| 2 )). From the numerical analysis point of view, the gradient descent step within the parentheses is the forward step, and the application of proximity operator prox tR 2 is the backward step.
By denoting z (k) = f (k) −τA T (Af (k) −g) and α = τλ, then Eq (16) becomes Specially, the computation of proximity operator prox tR 2 described in Eq (17) can be considered as an HDTV2-based denoising problem with z (k) being the noisy observation. Considering that f is a digital image whose intensities are finite and bounded, we can equivalently convert the unconstrained denoising problem Eq (17) to a constrained one, which is formulated as N} is a convex closed set which enforces bounded constraint on the solution, a and b are two constants, and ι s (f) is the indicator function of S which takes value 0 for f2S and 1 otherwise. As referred in [19], || Á || 1,2 is the dual norm of || Á || 1,2 , then the HDTV2 regularizer in Eq (18) can be equivalently written as where B 2 = {ω = (ω 1 ,ω 2,Á Á Á, ω N )2X 2 | ||ω i || 2 1,8i = 1,Á Á Á, N} denotes the l 1 −l 2 unit norm ball, and h Á,Á i 2 denotes the inner product in Euclidean space < N . Hereafter, we denote Η 2 = Λ 2 V 2 for simplicity. Thus, using Eq (19), we can solve the problem Eq (18) by the following minimax problem: Noting that the function E(f,ω) can be also written as Due to the fact that E(f,ω) is strictly convex in f and concave in ω, we can exchange the order of the minimum and maximum. Then, the saddle-point ðf ðkþ1Þ ;ωÞ of Eq (20), which makes can be obtained by solving where P S denotes the orthogonal projection onto the convex closed set S, andω is the maximizer of the dual problem, which is formulated aŝ Since the objective function h(ω) in Eq (24) is differentiable and has Lipschitz continuous gradient, as detailed in [20], we can compute the gradient of h(ω) as rhðωÞ ¼ a H 2 P S ðz ðkÞ À a H Τ 2 ωÞ: In the following proposition, we will derive an upper bound of the Lipschitz constant of rh (ω).
To accelerate the convergence rate, we use the Nesterov's iterative method [22] to solve problem Eq (24), which exhibits convergence rates of one order higher than the standard gradient ascent method. In detail, we employ the operator P B 2 that returns the orthogonal projection onto the l 1 −l 2 unit norm ball B 2 to compute the solution of the dual problem Eq (24), and then obtain the solution of Eq (18) by using Eq (23).
In summary, a detailed description of our complete forward-backward splitting (FBS) algorithm for HDTV2-based image reconstruction is provided in Table 1.

Convergence Analysis
In this section, we focus on analyzing the convergence of the proposed FBS algorithm for HDTV2-based image reconstruction problem relying on the properties of averaged and firmly non-expansive operators, which is very different from the convergence proof of the FBS algorithm [17]. To simplify our analysis, according to the iteration form of FBS in Eqs (16) and (17), we first introduce the following notations: where TðÁÞ ¼ D S HDTV2 ðS G ðÁÞÞ. Specifically, Eq (31) indicates that f (k) = T k (f (0) ), for any initial vector f (0) . It is clear from Eq (30) that S G (Á) = (I−τrR 1 )(Á) and S HDTV2 (Á) = prox αR2 (Á). According to the well-known Krasnoselskii-Mann (KM) theorem [21], if T is an averaged non-expansive operator and has fixed points, then the sequence {T k (x (0) )} k converges to a fixed point of T, for any initial vector x (0) .
In order to analyze the convergence of the proposed FBS algorithm, we need to show that the operator T defined in Eq (31) is averaged non-expansive and the set of the fixed points of T is nonempty.

backward step
and I À 1 L rR 1 is also non-expansive, thus S G is τL-averaged non-expansive.
Based on Lemma 1, we can obtain that the operator T defined in Eq (31) is 1 þ tL 2 -averaged non-expansive. Next, we will show that the objective function R(f) described in Eq (12) is coercive under certain conditions and the set of the fixed points of T is nonempty in the following Lemma. Proof. It is clear that D is not a full-rank matrix. Thus, the null space of D is the set {c1} with c being a scalar, where 12< N . With the HDTV2 regularizer described in Eq (11), we further have that where |Á| denotes the absolute value of corresponding scalar and ||Á|| 1 denotes the l 1 norm in Euclidean space. Using Eq (33), we obtain that Let y 1 then we can have that Null(U) = {0}, i.e. Null(A)\Null(D) = {0}. As if for arbitrary f 0 2Null(U), then Af 0 = 0 and Df 0 = 0. Based on the null space of D is the set {c1}, A = PF2C M×N is the sensing matrix with F2C N×N denoted the two-dimensional discrete Fourier transform (DFT) matrix and P2C M×N denoted the selection matrix containing M rows of the identity matrix of order N(Here we note that P must include the row corresponding to the (0,0) Fourier coefficient), then we have f 0 = c1 but Af 0 6 ¼ 0 when c 6 ¼ 0, thus f 0 = 0, i.e. Null(U) = {0}. Then, matrix U is full rank. Thus, it follows that when ||f|| 2 !1 in Eq (34), either ||y 1 || 2 !1 or ||y 2 || 2 !1, then R(f)!1. Based on that, Lemma 2 follows. Based on Lemma 2, we obtain that the objective function R(f) is coercive and the set of the minimizers of R(f) is nonempty. Denotingf as a minimizer of R(f) and using Eq (30), we havê f ¼ S HDTV2 ðS G ðf ÞÞ ¼ Tðf Þ so thatf is a fixed point of T. Thus, the set of the fixed points of T is nonempty.
Based on above analysis, we obtain that the operator T is 1 þ tL 2 -averaged non-expansive and the set of the fixed points of T is nonempty. According to the KM theorem [21], we can have that the sequence {f (k) } k generated by Eq (30) converges to a fixed point of T as k!1, for any initial vector f (0) 2< N .

Experimental Results
In this section, we present many numerical experiments to illustrate the performance of our proposed FBS algorithm for HDTV2-based image CS reconstruction problem. Specifically, we compare our proposed FBS method for HDTV2 (denoted by HDTV2-FBS) with the fast iterative shrinkage-thresholding algorithm (FISTA) for TV [18] (denoted by TV-FISTA) and stateof-the-art MM method for HDTV2 [12] (denoted by HDTV2-MM). Finally, the reconstruction performance of each method is evaluated in terms of peak signal to noise ratio (PSNR) and structural similarity index (SSIM), which are defined as PSNR ¼ 10log 10 max i;j jf ði; jÞj wheref and f2< r×c are the reconstructed image and the original image, respectively, ||Á|| F is the Frobenius norm, mf and μ f are the average values off and f, respectively, s 2 f and s 2 f are the variance off and f, respectively, sf f is the covariance off and f, c 1 and c 2 are two small positive constants.

Experimental Setting
In our experiments, the test images including magnetic resonance (MR) image and cell images are Brain, Fluorescent cell, CIL 240 and CIL 248, respectively, where the Brain image is available at http://www.strokeeducation.co.uk/?page_id=300, the Fluorescent cell image is available at http://rsb.info.nih.gov/ij, and the CIL 240 and CIL 248 images are part of the biomedical image database [25]. Moreover, the intensities of all test images are normalized in the range of [0,1]. Specifically, the compressive samples for each test are acquired by first applying variable density random Fourier encoding and then adding the Gaussian white noise corresponding to a SNR level of 30 dB.
To compare the performance fairly, we constrain the reconstructed images of all methods to lie in the convex closed set S = {f2< N |f i 2[0,1],8i = 1,. . ., N}. Moreover, all parameters of the TV-FISTA and HDTV2-MM methods are set to be the default values in [18], [12], but the regularization parameter for each method is chosen to give the best PSNR results. All the methods are implemented in MATLAB 7.11.0 (R2010b) and run on an Intel (R) Xeon(R) CPU 2.67-GHz PC with 4-GB memory.

Result Analysis
To characterize the performance quantitatively, we provide the PSNR, SSIM and running time results of different methods under different sampling ratios in Table 2. It is clear to observe that the proposed HDTV2-FBS method always provides the best PSNR and SSIM results, where the PSNR result of HDTV2-FBS method is even 1 dB higher than the HDTV2-MM method for the case of the CIL 240 image. The TV-FISTA method provides very competitive PSNR and SSIM results with the HDTV2-MM method, while the TV-FISTA method takes much less running time than the HDTV2-MM method. Considering the comparisons of running time, the proposed HDTV2-FBS method always takes much less running time than the HDTV2-MM method. More specifically, the TV-FISTA method consistently takes the shortest running time. In fact, the computation of the HDTV2-FBS method is more complex than that of the TV-FISTA method, thus it is more valuable for our proposed HDTV2-FBS method to spend a little more time than the TV-FISTA method for achieving better reconstruction results. Based on these numerical results, they clearly demonstrate that our proposed HDTV2-FBS method is more effective and efficient for image CS reconstruction.
Next, from both reconstruction quality and convergence speed point of view, we representatively show some examples with different sampling ratios to illustrate the performance of our proposed HDTV2-FBS method. Fig 1 shows the reconstructed images of different methods on Brain image with sampling ratio 20% and noise level 30 dB. Obviously, we observe that the HDTV2-MM method and our proposed HDTV2-FBS method can both preserve some fine details in the reconstructed images. Furthermore, in order to highlight the differences among these methods, we zoom in the corresponding region marked by the white box in Fig 1A and present the zoomed details in Fig 1E-1H. By inspecting these images carefully, the HDTV2-MM method provides better preservation of smooth image regions but loses some fine details, and our proposed HDTV2-FBS method can better preserve some fine details, while the TV-FISTA method suffers from severe blocky artifacts in smooth regions. More specifically, we plot the cost of the objective function versus the number of iterations for the HDTV2-MM and HDTV2-FBS methods in Fig 2A to illustrate which one is better for cost minimization. As shown in Fig 2A, we can obtain that the proposed HDTV2-FBS method is more accurate than the HDTV2-MM method in terms of cost minimization, where the proposed HDTV2-FBS method exhibits a much faster convergence speed than the HDTV2-MM method. Meanwhile, we also plot the convergence curve of PSNR versus the number of iterations in Fig 2B. It is clear that our proposed HDTV2-FBS method not only provides better PSNR results than the TV-FISTA and HDTV2-MM methods, but also exhibits a much faster convergence speed than the HDTV2-MM method. The reconstructed results of different methods on image CIL 240 with sampling 30% and noise level 30 dB are shown in Fig 3. Since the CIL 240 image consists of many filament-like features, it is more desirable to well preserve these elongated features as far as possible. As we can clearly see from the residual images of TV-FISTA, HDTV2-MM and HDTV2-FBS results shown in Fig 3E-3G, the HDTV2-MM method could not well preserve the filament-like features and still leaves some noise in the smooth regions of the reconstructed image. The TV-FISTA method can remove the residual noise but fails to preserve filament-like features and produces some staircase effect in the smooth regions. By contrast, our proposed HDTV2-FBS method can better remove the residual noise and at the same time preserve the elongated filament-like features more sharpness, which makes the reconstructed image more natural. Meanwhile, considering the convergence curves of PSNR and relative error plotted in Fig 4, we can observe that our proposed HDTV2-FBS method provides the best PSNR and relative error results. Specifically, the proposed HDTV2-FBS method shows the fastest convergence speed, where it is about two folds faster than that of the HDTV2-MM method.
Finally, we show the reconstructed results of different methods on image CIL 248 with sampling 50% and noise level 30 dB in Fig 5. In this sense, the sampling ratio is enough large to   ensure a good reconstruction quality. We can see from the reconstructed images that the TV-FISTA method, the HDTV2-MM method and our proposed HDTV2-FBS method provide very good reconstruction of image features. Moreover, from the convergence speed and accuracy point of view, we clearly see from Fig 6 that our proposed HDTV2-FBS method provides the best PSNR and relative error results, and shows the fastest convergence speed, where the convergence speed of our proposed HDTV2-FBS method is much faster than that of the HDTV2-MM method. Thus, we can conclude that our proposed HDTV2-FBS method is more efficient for image CS reconstruction tasks in terms of accuracy and convergence speed.

Conclusions
In this paper, we have introduced an efficient HDTV2 algorithm for image CS reconstruction using forward-backward splitting scheme. Under the spectral decomposition framework, HDTV2 is reinterpreted as a novel weighted L 1 -L 2 mixed norm of the second degree image derivatives. Based on the resulting equivalent formulation of HDTV2, a computationally efficient forward-backward splitting scheme is thus designed to solve the HDTV2-based reconstruction problem. More specially, we prove in detail the convergence of the proposed FBS algorithm. Furthermore, many experimental results on MR image and cell images show that the proposed scheme provides better reconstruction results and especially exhibits much faster convergence speed than the iteratively reweighted MM method. Specifically, the proposed scheme can also preserve the elongated filament-like features more sharpness for cell images.

Ethics Statement
These Brain, Fluorescent cell, CIL 240 and CIL 248 images were used in this paper to verify the performance of our method. These Brain, Fluorescent cell, CIL 240 and CIL 248 images are publicly available for image compressive sensing research, and the consent was not needed. These images and experimental results are reported in this paper without any commercial purpose.