Compressed Sensing MRI Reconstruction Algorithm Based on Contourlet Transform and Alternating Direction Method

. Compressed sensing (CS) based methods have recently been used to reconstruct magnetic resonance (MR) images from undersampled measurements, which is known as CS-MRI. In traditional CS-MRI, wavelet transform can hardly capture the information of image curves and edges. In this paper, we present a new CS-MRI reconstruction algorithm based on contourlet transform and alternating direction method (ADM). The MR images are firstly represented by contourlet transform, which can describe the images’ curves and edges fully and accurately. Then the MR images are reconstructed by ADM, which is an effective CS reconstruction method. Numerical results validate the superior performance of the proposed algorithm in terms of reconstruction accuracy and computation time.


Introduction
CS is a new sampling and compression theory.It utilizes the sparseness of a signal in a particular domain and can reconstruct the signal from significantly fewer samples than Nyquist sampling, which has been the fundamental principle in signal processing for many years [1][2][3].Due to the above advantages, CS has received considerable attentions in many areas, one of which is MRI reconstruction [4,5].MRI is safer, more frequent, and accurate for clinical diagnosis.However, conventional MRI needs to spend much time scanning body regions, causing the expensive cost and the nonidealized space resolution.In the meantime, the physiological property in the tested body will make the image blurry and distortional.Therefore, under the premise of guaranteeing the image quality, speeding up the MRI compression and reconstruction has been the powerful impetus to promote the development of MRI techniques.
For CS-MRI, there are two key points that need further investigation.The first one is sparse transform.In MRI reconstruction, the MR images themselves are not sparse but have sparse representations in some transform domains.In traditional CS-MRI, wavelet transform is commonly used as a sparse transform [6,7].However, as the limitations of direction, wavelet transform can hardly capture the information of image curves and edges fully and accurately.In contrast, curves and edges are mainly features of MR images.Therefore, more effective sparse transform should be considered for CS-MRI.Contourlet transform, also known as Pyramid Directional Filter Bank (PDFB), is put forward to make up for the inadequacy of the wavelet transform [8].Contourlet transform can describe the image's contour and directional texture information fully and accurately since it realizes any directional decomposition at each scale.Furthermore, contourlet is constructed directly in a discrete domain and has low computing complexity.Thus, contourlet transform can be easily implemented for MR images [9,10].
The second one is the reconstruction algorithm.In recent years, a number of algorithms have been put forward for the signal reconstruction in CS, for example, interior-point algorithm [11], iterative shrinkage/thresholding algorithm (ISTA) [12], fast iterative shrinkage/threshold algorithm (FISTA) [13], Sparse Reconstruction by Separable Approximation (SpaRSA) [14].But not all of these algorithms are suitable for CS-MRI since the dimensions of the MR images are huge.Alternating direction method (ADM) is an efficient CS reconstruction algorithm that has a faster convergence speed than some traditional methods [15].Meanwhile, ADM is able to solve large-scale CS problem since all the iterations of ADM only contain the first-order information of the objective function, which have low computing complexity.
In this paper, we present a new CS-MRI reconstruction algorithm based on contourlet transform and ADM.The proposed algorithm can recover the curves and edges of a MR image more precisely and suit large-scale MRI reconstruction.
The organization of the rest of this paper is as follows: in Section 2, we first introduce CS-MRI model briefly and then present our new algorithm.Numerical results will demonstrate the effectiveness of the proposed algorithm in Section 3. Finally, we conclude the paper in Section 4.

Contourlet-Based CS-MRI Model.
The basic problem of CS is to recover a signal x from underdetermined linear measurement y = Φx where y ∈ R  , x ∈ R  , and Φ ∈ R × ,  < .This underdetermined linear system has infinite solutions when seen from the aspect of algebra.However, according to the CS theory, under the assumption that x is sparse, x can be reconstructed by the following optimization problem: where ‖x‖ 0 is the ℓ 0 norm which means the nonzero numbers of x.
Problem (1) is difficult to solve since it is NP-hard.It can be relaxed as the following convex problem: where ‖x‖ 1 = ∑  |  | is the sum of absolute values of x.Problem (2) is a convex optimization problem and can be solved by many algorithms.
Problems (1) and ( 2) are under the assumption that x is sparse.However, in many applications, the signal itself is not sparse but has a sparse representation in some transform domains.For example x = Ψ, where x is the original signal which is not sparse and  is the sparse coefficient with respect to the sparse transform matrix Ψ.In this case, CS model should be y = Φ = ΦΨ * x = Ax, where Ψ * denotes the inverse of Ψ.The optimization problem should be changed as follows: (3)

ADM for Contourlet-Based CS-MRI.
In this subsection, we solve (3) by ADM.We first introduce an auxiliary variable r and transform (3) into an equivalent problem: The augmented Lagrangian function of problem ( 4) is given by min x,r { (x, r) =     Ψ * x    1 −  T (Ax + r − y) where  ∈ R  is a Lagrangian multiplier and  > 0 is a penalty parameter.
The objective function (x k , r) is only connected with r; then (5) is equivalent to For problem (7), since (/r)(x k , r) = −( k ) T + (Ax k + r − y) T = 0, we have r =  k / − (Ax k − y).The minimization of (7) with respect to r is shown by where Ρ   is the projection onto the set B  : {r ∈ R  : ‖r‖ ≤ }.
Then we fix r = r k+1 ,  =   ; that is, The objective function (x, r k+1 ) is only relevant to x; then problem ( 5) is equivalent to min x  (x, r k+1 ) .
Equation (10) Require: A, y, r 0 , x 0 ,  0 ,  > 0,  > 0, Γ > 0. Ensure: x while "stopping criterion is not met" do where g k ΔA T (Ax k + r k+1 − y −  k /) is the gradient of quadratic items about (12), Γ > 0 is a parameter, and the first equation of ( 12) is based on the Taylor expansion.Then ( 11) is equivalent to Problem ( 13) is equivalent approximately to which has a close form solution by shrinkage (or soft thresholding) formula Finally, we update the multiplier  through where  > 0 is a constant.Now, we present alternating direction method for problem (3) as Algorithm 1.

Numerical Experiments
In this section, we present numerical results to illustrate the performance of the proposed algorithm for MRI reconstruction.All experiments are made by using MATLAB 7.8.0 on the PC with Intel Core 3.4 GHz and 4 G memory.
We compare the proposed Algorithm 1 with the stateof-the-art MRI reconstruction algorithms: SparseMRI [5], ICOTA [9], FICOTA [10], and SpaRSA [14].SparseMRI is based on conjugate gradient (CG), ICOTA is based on contourlet transform and ISTA, FICOTA is based on contourlet transform and FISTA, and SpaRSA is an effective algorithm to solve the CS problem.We also compare the reconstruction performance of contourlet and wavelet.We quantify the reconstruction performance by peak signal to noise ratio (PSNR) and CPU time.PSNR is defined as where MSE = (1/) ∑ −1 =0 ∑ −1 =0 ( ori (, ) −  rec (, )) 2 ,  ori and  rec are the original image and reconstructed image, respectively, and ,  are the size of the images.
In the following experiments, we choose the measurement matrix Φ as a partial Fourier transform matrix with  rows and  columns and define the sampling ratio as /.The MR scanning time is less if the sampling ratio is lower.We use the variable density sampling strategy just as the previous works in papers [4,5,9,10], which randomly choose more Fourier coefficients from low frequencies and less coefficients from high frequencies.
In the first experiment, we use the MR image as Figure 1(a) shows, and the variable density sampling pattern as Figure 1(b) shows.For wavelet transform, we use Daubechies wavelet with 4 vanishing moments, and contourlet transform with decomposition [5,4,4,3], just the same as paper [9].
Figure 1 shows the reconstructed images using different algorithms.Table 1 summarizes the comparisons of different algorithms.From Table 1 we can see that Algorithm 1 with contourlet transform outperforms Algorithm 1 with wavelet transform in terms of PSNR, although its running time is slightly slower, and Algorithm 1 outperforms other algorithms in terms of both PSNR and CPU time.
In the second experiment, we illustrate the reconstruction performance of each algorithm as the sampling rate varies from 0.1 to 0.9.
Figure 2 shows the variations of PSNR and CPU time of CS-MRI reconstruction versus sampling rates for different algorithms.Figure 2(a) shows that PSNR of Algorithm 1 with contourlet is better than other algorithms with sampling growing.Furthermore, Figure 2(a) also shows that the advantage of contourlet becomes less obvious with the increase of sampling rate.Figure 2(b) indicates that Algorithm 1 with contourlet is slightly slower than Algorithm 1 with wavelet but is much faster than other algorithms.

Conclusion
In this paper, we propose a novel algorithm based on contourlet transform and the classic alternating direction method to solve CS-MRI problem.The proposed algorithm has low computational complexity and is suitable for large scale problem.Our numerical results show that the proposed algorithm compares favorably with these algorithms referred to in terms of PSNR and CPU time.

Table 1 :
Comparisons of different algorithms.