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

Authors: Zhenyu Hu; Qiuye Wang; Congcong Ming; Lai Wang; Yuanqing Hu; Jian Zou
Journal: Journal of Electrical and Computer Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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–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 Section2, 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.


## 2. ADM for Contourlet-Based CS-MRI

### 2.1. Contourlet-Based CS-MRI Model

The basic problem of CS is to recover a signalx from underdetermined linear measurement y
=
Φ
x where y
∈
R
m, x
∈
R
n, and Φ
∈
R
m
×
n,m
<
n. 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:(1)
m
i
n
x
x
0
s.t.
y
=
Φ
x
,where x
0 is the l
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:(2)
m
i
n
x
x
1
s.t.
y
=
Φ
x
,where x
1
=
∑
i
|
x
i
| 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
=
A
x, where Ψ
∗ denotes the inverse of Ψ. The optimization problem should be changed as follows:(3)
m
i
n
x
Ψ
∗
x
1
s.t.
y
=
A
x
.


### 2.2. 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: (4)
min
x
,
r
⁡
Ψ
∗
x
1
:
A
x
+
r
=
y
,
r
≤
δ
.The augmented Lagrangian function of problem (4) is given by(5)
min
x
,
r
⁡
f
x
,
r
=
Ψ
∗
x
1
-
λ
T
A
x
+
r
-
y
+
β
2
A
x
+
r
-
y
2
:
r
≤
δ
,where λ
∈
R
m is a Lagrangian multiplier and β
>
0 is a penalty parameter.If we fixx
=
x
k, λ
=
λ
k; that is,(6)
f
x
k
,
r
=
Ψ
∗
x
k
1
-
λ
k
T
A
x
k
+
r
-
y
+
β
2
A
x
k
+
r
-
y
2
.The objective function f
(
x
k
,
r
) is only connected with r; then (5) is equivalent to(7)
min
r
⁡
f
x
k
,
r
:
r
≤
δ
.For problem (7), since d
/
d
r
f
(
x
k
,
r
)
=
-
(
λ
k
)
T
+
β
(
A
x
k
+
r
-
y
)
T
=
0, we have r
=
λ
k
/
β
-
(
A
x
k
-
y
). The minimization of (7) with respect to r is shown by(8)
r
k
+
1
=
Ρ
B
δ
λ
k
β
-
A
x
k
-
y
,where Ρ
B
δ is the projection onto the set B
δ
:
r
∈
R
m
:
r
≤
δ.Then we fixr
=
r
k
+
1, λ
=
λ
κ; that is,(9)
f
x
,
r
k
+
1
=
Ψ
∗
x
1
-
λ
k
T
A
x
+
r
k
+
1
-
y
+
β
2
A
x
+
r
k
+
1
-
y
2
.The objective functionf
(
x
,
r
k
+
1
) is only relevant to x; then problem (5) is equivalent to (10)
min
x
⁡
f
x
,
r
k
+
1
.Equation (10) is equivalent approximately to(11)
min
x
⁡
Ψ
∗
x
1
+
β
2
A
x
+
r
k
+
1
-
y
-
λ
k
β
2
.Leth
(
x
)
=
β
/
2
A
x
+
r
k
+
1
-
y
-
λ
k
/
β
2; then (12)
h
x
=
β
2
A
x
k
+
r
k
+
1
-
y
-
λ
k
β
2
+
2
A
T
A
x
k
+
r
k
+
1
-
y
-
λ
k
β
T
+
1
2
x
-
x
k
T
2
A
T
A
x
-
x
k
=
β
g
k
T
x
-
x
k
+
1
2
A
x
-
x
k
2
=
β
g
k
T
x
-
x
k
+
1
2
Γ
x
-
x
k
2
,where g
k
Δ
_
_
A
T
(
A
x
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(13)
min
x
⁡
Ψ
∗
x
1
+
β
g
k
T
x
-
x
k
+
1
2
Γ
x
-
x
k
2
.Problem (13) is equivalent approximately to (14)
min
x
⁡
Ψ
∗
x
1
+
β
2
Γ
Ψ
∗
x
-
x
k
-
Γ
g
k
,which has a close form solution by shrinkage (or soft thresholding) formula(15)
x
k
+
1
=
S
h
r
i
n
k
x
k
-
Γ
g
k
,
Γ
β
Ψ
∗
.Finally, we update the multiplierλ through(16)
λ
k
+
1
=
λ
k
-
γ
β
A
x
k
+
1
+
r
k
+
1
-
y
,where γ
>
0 is a constant.Now, we present alternating direction method for problem (3) as Algorithm 1.Algorithm 1
Require: A
,
y
,
r
0
,
x
0
,
λ
0
,
β
>
0
,
γ
>
0
,
Γ
>
0.
Ensure: x
while “stopping criterion is not met” do
r
k
+
1
=
Ρ
B
δ
λ
k
β
-
A
x
k
-
y;
x
k
+
1
=
Shrink
x
k
-
Γ
g
k
,
Γ
β
Ψ
∗;
λ
k
+
1
=
λ
k
-
γ
β
A
x
k
+
1
+
r
k
+
1
-
y;
end while


## 2.1. Contourlet-Based CS-MRI Model

The basic problem of CS is to recover a signalx from underdetermined linear measurement y
=
Φ
x where y
∈
R
m, x
∈
R
n, and Φ
∈
R
m
×
n,m
<
n. 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:(1)
m
i
n
x
x
0
s.t.
y
=
Φ
x
,where x
0 is the l
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:(2)
m
i
n
x
x
1
s.t.
y
=
Φ
x
,where x
1
=
∑
i
|
x
i
| 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
=
A
x, where Ψ
∗ denotes the inverse of Ψ. The optimization problem should be changed as follows:(3)
m
i
n
x
Ψ
∗
x
1
s.t.
y
=
A
x
.


## 2.2. 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: (4)
min
x
,
r
⁡
Ψ
∗
x
1
:
A
x
+
r
=
y
,
r
≤
δ
.The augmented Lagrangian function of problem (4) is given by(5)
min
x
,
r
⁡
f
x
,
r
=
Ψ
∗
x
1
-
λ
T
A
x
+
r
-
y
+
β
2
A
x
+
r
-
y
2
:
r
≤
δ
,where λ
∈
R
m is a Lagrangian multiplier and β
>
0 is a penalty parameter.If we fixx
=
x
k, λ
=
λ
k; that is,(6)
f
x
k
,
r
=
Ψ
∗
x
k
1
-
λ
k
T
A
x
k
+
r
-
y
+
β
2
A
x
k
+
r
-
y
2
.The objective function f
(
x
k
,
r
) is only connected with r; then (5) is equivalent to(7)
min
r
⁡
f
x
k
,
r
:
r
≤
δ
.For problem (7), since d
/
d
r
f
(
x
k
,
r
)
=
-
(
λ
k
)
T
+
β
(
A
x
k
+
r
-
y
)
T
=
0, we have r
=
λ
k
/
β
-
(
A
x
k
-
y
). The minimization of (7) with respect to r is shown by(8)
r
k
+
1
=
Ρ
B
δ
λ
k
β
-
A
x
k
-
y
,where Ρ
B
δ is the projection onto the set B
δ
:
r
∈
R
m
:
r
≤
δ.Then we fixr
=
r
k
+
1, λ
=
λ
κ; that is,(9)
f
x
,
r
k
+
1
=
Ψ
∗
x
1
-
λ
k
T
A
x
+
r
k
+
1
-
y
+
β
2
A
x
+
r
k
+
1
-
y
2
.The objective functionf
(
x
,
r
k
+
1
) is only relevant to x; then problem (5) is equivalent to (10)
min
x
⁡
f
x
,
r
k
+
1
.Equation (10) is equivalent approximately to(11)
min
x
⁡
Ψ
∗
x
1
+
β
2
A
x
+
r
k
+
1
-
y
-
λ
k
β
2
.Leth
(
x
)
=
β
/
2
A
x
+
r
k
+
1
-
y
-
λ
k
/
β
2; then (12)
h
x
=
β
2
A
x
k
+
r
k
+
1
-
y
-
λ
k
β
2
+
2
A
T
A
x
k
+
r
k
+
1
-
y
-
λ
k
β
T
+
1
2
x
-
x
k
T
2
A
T
A
x
-
x
k
=
β
g
k
T
x
-
x
k
+
1
2
A
x
-
x
k
2
=
β
g
k
T
x
-
x
k
+
1
2
Γ
x
-
x
k
2
,where g
k
Δ
_
_
A
T
(
A
x
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(13)
min
x
⁡
Ψ
∗
x
1
+
β
g
k
T
x
-
x
k
+
1
2
Γ
x
-
x
k
2
.Problem (13) is equivalent approximately to (14)
min
x
⁡
Ψ
∗
x
1
+
β
2
Γ
Ψ
∗
x
-
x
k
-
Γ
g
k
,which has a close form solution by shrinkage (or soft thresholding) formula(15)
x
k
+
1
=
S
h
r
i
n
k
x
k
-
Γ
g
k
,
Γ
β
Ψ
∗
.Finally, we update the multiplierλ through(16)
λ
k
+
1
=
λ
k
-
γ
β
A
x
k
+
1
+
r
k
+
1
-
y
,where γ
>
0 is a constant.Now, we present alternating direction method for problem (3) as Algorithm 1.Algorithm 1
Require: A
,
y
,
r
0
,
x
0
,
λ
0
,
β
>
0
,
γ
>
0
,
Γ
>
0.
Ensure: x
while “stopping criterion is not met” do
r
k
+
1
=
Ρ
B
δ
λ
k
β
-
A
x
k
-
y;
x
k
+
1
=
Shrink
x
k
-
Γ
g
k
,
Γ
β
Ψ
∗;
λ
k
+
1
=
λ
k
-
γ
β
A
x
k
+
1
+
r
k
+
1
-
y;
end while


## 3. 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 Algorithm1 with the state-of-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(17)
P
S
N
R
=
10
log
10
⁡
25
5
2
M
S
E
,where M
S
E
=
1
/
m
n
∑
i
=
0
m
-
1
∑
j
=
0
n
-
1
I
o
r
i
(
i
,
j
)
-
I
r
e
c
(
i
,
j
)
2, I
o
r
i and I
r
e
c are the original image and reconstructed image, respectively, and m, n are the size of the images.In the following experiments, we choose the measurement matrixΦ as a partial Fourier transform matrix with m rows and n columns and define the sampling ratio as m
/
n. 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 Figure1(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
Comparison of CS-MRI reconstruction results obtained from different algorithms. (a) Original MR image. (b) Variable density sampling pattern with sampling rate 0.2. ((c)–(h)) Reconstructed result from Algorithm1 with contourlet transform, SparseMRI, ICOTA, FICOTA, SpaRSA, and Algorithm 1 with wavelet transform.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure1 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.Table 1
Comparisons of different algorithms.
PSNR (dB)
CPU time (s)
Algorithm1 with Contourlet
48.28
4.52
SparseMRI
36.71
16.17
ICOTA
45.79
15.46
FICOTA
46.93
7.79
SpaRSA
38.46
14.52
Algorithm1 with wavelet
42.27
4.21In the second experiment, we illustrate the reconstruction performance of each algorithm as the sampling rate varies from 0.1 to 0.9.Figure2 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.Figure 2
Comparison of CS-MRI reconstruction results obtained from different algorithms with different sampling rates. The results are average of 50 runs. (a) PSNR versus sampling rate. (b) CPU time versus sampling rate.
(a)
(b)


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



---

*Source: 175231-2015-11-04.xml*