# Approximate Sparsity and Nonlocal Total Variation Based Compressive MR Image Reconstruction

Authors: Chengzhi Deng; Shengqian Wang; Wei Tian; Zhaoming Wu; Saifeng Hu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Recent developments in compressive sensing (CS) show that it is possible to accurately reconstruct the magnetic resonance (MR) image from undersampledk-space data by solving nonsmooth convex optimization problems, which therefore significantly reduce the scanning time. In this paper, we propose a new MR image reconstruction method based on a compound regularization model associated with the nonlocal total variation (NLTV) and the wavelet approximate sparsity. Nonlocal total variation can restore periodic textures and local geometric information better than total variation. The wavelet approximate sparsity achieves more accurate sparse reconstruction than fixed wavelet l
0 and l
1  norm. Furthermore, a variable splitting and augmented Lagrangian algorithm is presented to solve the proposed minimization problem. Experimental results on MR image reconstruction demonstrate that the proposed method outperforms many existing MR image reconstruction methods both in quantitative and in visual quality assessment.


---

## Body

## 1.  Introduction

Magnetic resonance imaging (MRI) is a noninvasive and nonionizing imaging processing. Due to its noninvasive manner and intuitive visualization of both anatomical structure and physiological function, MRI has been widely applied in clinical diagnosis. Imaging speed is important in many MRI applications. However, both scanning and reconstruction speed of MRI will affect the quality of reconstructed image. In spite of advances in hardware and pulse sequences, the speed, at which the data can be collected in MRI, is fundamentally limited by physical and physiological constraints. Therefore many researchers are seeking methods to reduce the amount of acquired data without degrading the image quality [1–3].In recent years, the compressive sensing (CS) framework has been successfully used to reconstruct MR images from highly undersampledk-space data [4–9]. According to CS theory [10, 11], signals/images can be accurately recovered by using significantly fewer measurements than the number of unknowns or than mandated by traditional Nyquist sampling. MR image acquisition can be looked at as a special case of CS where the sampled linear combinations are simply individual Fourier coefficients (k-space samples). Therefore, CS is claimed to be able to make accurate reconstructions from a small subset of k-space data. In compressive sensing MRI (CSMRI), we can reconstruct a MR image with good quality from only a small number of measurements. Therefore, the application of CS to MRI has potential for significant scan time reductions, with benefits for patients and health care economics.Because of the ill-posed nature of the CSMRI reconstruction problem, regularization terms are required for a reasonable solution. In existing CSMRI models, the most popular regularizers arel
0, l
1 sparsity [4, 9, 12] and total variation (TV) [3, 13]. The l
0 sparsity regularized CSMRI model can be understood as a penalized least square with l
0 norm penalty. It is well known that the complexity of this model is proportional with the number of variables. Particularly when the number is large, solving the model generally is intractable. The l
1 regularization problem can be transformed into an equivalent convex quadratic optimization problem and, therefore, can be very efficiently solved. And under some conditions, the resultant solution of l
1 regularization coincides with one of the solutions of l
0 regularization [14]. Nevertheless, while l
1 regularization provides the best convex approximation to l
0 regularization and it is computationally efficient, the l
1 regularization often introduces extra bias in estimation and cannot reconstruct an image with the least measurements when applied to CSMRI [15]. In recent years, the l
q
(
0
<
q
<
1
) regularization [16, 17] was introduced into CSMRI, since l
q regularization can assuredly generate much sparser solutions than l
1 regularization. Although the l
q regularizations achieve better performance, they always fall into local minima. Moreover, which q should yield a best result is also a problem. Trzasko and Manduca [18] proposed a CSMRI paradigm based on homotopic approximation of the l
0 quasinorm. Although this method has no guarantee of achieving a global minimum, it achieves accurate MR image reconstructions at higher undersampling rates than l
1 regularization. And it was faster than those l
q regularization methods. Recently, Chen and Huang [19] accelerated MRI by introducing the wavelet tree structural sparsity into the CSMRI.Despite high effectiveness in CSMRI recovery, sparsity and TV regularizers often suffer from undesirable visual artifacts and staircase effects. To overcome those drawbacks, some hybrid sparsity and TV regularization methods [5–8] are proposed. In [5], Huang et al. proposed a new optimization algorithm for MR image reconstruction method, named fast composite splitting algorithm (FCSA), which is based on the combination of variable and operator splitting techniques. In [8], Yang et al. proposed a variable splitting method (RecPF) to solve hybrid sparsity and TV regularized MR image reconstruction optimization problem. Ma et al. [20] proposed an operator splitting algorithm (TVCMRI) for MR reconstruction. In order to deal with the problem of low and high frequency coefficients measurement, Zhang et al. [6] proposed a new so-called TVWL2-L1 model which measures low frequency coefficients and high frequency coefficients with l
2 norm and l
1 norm. In [7], an experimental study on the choice of CSMRI regularizations was given. Although the classical TV regularization performs well in CSMRI reconstruction while preserving edges, especially for cartoon-like MR images, it is well known that TV regularization is not suitable for images with fine details and it often tends to oversmooth image details and textures. Nonlocal TV regularization extends the classical TV regularization by nonlocal means filter [21] and has been shown to outperform the TV in several inverse problems such as image deonising [22], deconvolution [23], and compressive sensing [24, 25]. In order to improve the signal-to-noise ratio and preserve the fine details of MR images, Gopi et al. [26], Huang and Yang [27], and Liang et al. [28] have proposed nonlocal TV regularization based MR reconstruction and sensitivity encoding reconstruction.In this paper, we proposed a novel compound regularization based compressive MR image reconstruction method, which exploits the nonlocal total variation (NLTV) and the approximate sparsity prior. The approximate sparsity, which is used to replace the traditionall
0 regularizer and l
1 regularizer of compressive MR image reconstruction model, is sparser and much easier to be solved. The NLTV is much better than TV for preserving the sharp edges and meanwhile recovering the local structure details. In order to compound regularization model, we develop an alternative iterative scheme by using the variable splitting and augmented Lagrangian algorithm. Experimental results show that the proposed method can effectively improve the quality of MR image reconstruction. The rest of the paper is organized as follows. In Section 2 we review the compressive sensing and MRI reconstruction. In Section 3 we propose our model and algorithm. The experimental results and conclusions will be shown in Sections 4 and 5, respectively.


## 2. Compressive Sensing and MRI Reconstruction

Compressive sensing [10, 11], as a new sampling and compression theory, is able to reconstruct an unknown signal from a very limited number of samples. It provides a firm theoretical foundation for the accurate reconstruction of MRI from highly undersampled K-space measurements and significantly reduces the MRI scan duration.Supposeu
∈
R
N is a MR image and F
∈
R
M
×
N is a partial Fourier transform; then the sampling measurement b
∈
R
M of MR image u in K-space can be defined as
(1)
b
=
F
u
.
The compressive MR image reconstruction problem is to recover u given the measurement b and the sampling matrix F. Undersampling occurs whenever the number of K-space sample is less than the number of unknowns (
M
<
N
). In that case, the compressive MR image reconstruction is an underdetermined problem.In general, compressive sensing reconstructs the unknownsu from the measurements b by minimizing the l
0 norm of the sparsified image Φ
u, where Φ represents a sparsity transform for the image. In this paper, we choose the orthonormal wavelet transform as the sparsity transform for the image. Then the typical compressive MR image reconstruction is obtained by solving the following constrained optimization problem [4, 9, 12]:
(2)
min
⁡
u
∥
Φ
u
∥
0
s
.
t
.
b
=
F
u
.
However, in terms of computational complexity, the l
0 norm optimization problem (2) is a typical NP-hard problem, and it was difficult to solve. According to the certain condition of the restricted isometric property, the l
0 norm can be replaced by the l
1 norm. Therefore, the optimization problem (2) is relaxed to alternative convex optimization problem as follows:
(3)
min
⁡
u
∥
Φ
u
∥
1
s
.
t
.
b
=
F
u
.
When the measurements b are contaminated with noise, the typical compressive MR image reconstruction problem using l
1 relaxation of the l
0 norm is formulated as the following unconstrained Lagrangian version:
(4)
min
⁡
u
∥
F
u
-
b
∥
2
2
+
α
∥
Φ
u
∥
1
,
where α is a positive parameter.Despite high effectiveness of sparsity regularized compressive MR image reconstruction methods, they often suffer from undesirable visual artifacts such as Gibbs ringing in the result. Due to its desirable ability to preserve edges, total variation (TV) model is successfully used in compressive MR image reconstruction [3, 13]. But the TV regularizer has still some limitations that restrict its performance, which cannot generate good enough results for images with many small structures and often suffers from staircase artifacts. In order to combine the advantages of sparsity-based and TV model and avoid their main drawbacks, a TV regularizer, corresponding to a finite-difference for the sparsifying transform, is typically incorporated into the sparsity regularized compressive MR image reconstruction [5–8]. In this case the optimization problem is written as
(5)
min
⁡
u
∥
F
u
-
b
∥
2
2
+
α
∥
Φ
u
∥
1
+
β
∥
∇
u
∥
1
,
where β is a positive parameter. The TV was defined discretely as ∥
∇
u
∥
1
=
∑
i
(
|
∇
h
u
i
|
+
|
∇
v
u
i
|
), where ∇
h and ∇
v are the horizontal and the vertical gradient operators, respectively. The compound optimization model (5) is based on the fact that the piecewise smooth MR images can be sparsely represented by the wavelet and should have small total variations.


## 3. Proposed Model and Algorithm

As mentioned above, joint TV andl
1 norm minimization model is a useful way to reconstruct MR images. However, they have still some limitations that restrict their performance. l
0 norm needs a combinatorial search for its minimization and its too high sensibility to noise. l
1 problems can be very efficiently solved. But the solution is not sparse, which influences the performance of MRI reconstruction. The TV model can preserve edges, but it tends to flatten inhomogeneous areas, such as textures. To overcome those shortcomings, a novel method is proposed for compressive MR imaging based on the wavelet approximate sparsity and nonlocal total variation (NLTV) regularization, named WasNLTV.


### 3.1. Approximate Sparsity

The problems of usingl
0 norm in compressive MR imaging (i.e., the need for a combinatorial search for its minimization and its too high sensibility to noise) are both due to the fact that the l
0 norm of a vector is a discontinuous function of that vector. The same as [29, 30], our idea is to approximate this discontinuous function by a continuous one, named approximate sparsity function, which provides smooth measure of l
0 norm and better sparsity than l
1 regularizer.The approximate sparsity function is defined as(6)
ψ
σ
(
x
)
=
2
π
arctan
(
|
x
|
σ
2
)
,
x
∈
R
,
σ
∈
R
+
.
The parameter σ may be used to control the accuracy with which ψ
σ approximate the Kronecker delta. In mathematical terms, we have
(7)
lim
⁡
σ
→
0
ψ
σ
(
x
)
=
{
1
,
x
≠
0
,
0
,
x
=
0
.
Define the continuous multivariate approximate sparsity function Ψ
σ
(
x
) as
(8)
Ψ
σ
(
x
)
=
∑
i
=
1
m
ψ
σ
(
x
i
)
,
x
∈
R
m
×
1
.
It is clear from (7) that Ψ
σ
(
x
) is an indicator of the number of zero-entries in x for small values of σ. Therefore, l
0 norm can be approximate by
(9)
∥
x
∥
0
≈
Ψ
σ
(
x
)
=
∑
i
=
1
m
ψ
σ
(
x
i
)
.
Note that the larger the value of σ, the smoother the Ψ
σ
(
x
) and the worse the approximation to l
0 norm; the smaller the value of l
0 norm, the closer the behavior of Ψ
σ
(
x
) to l
0 norm.


### 3.2. Nonlocal Total Variation

Although the classical TV is surprisingly efficient for preserving edges, it is well known that TV is not suitable for images with fine structures, details, and textures which are very important to MR images. The NLTV is a variational extension of the nonlocal means filter proposed by Wang et al. [30]. NLTV uses the whole image information instead of using adjacent pixel information to calculate the gradients in regularization term. The NLTV has been proven to be more efficient than TV for improving the signal-to-noise ratio, on preserving not only sharp edges, but also fine details and repetitive patterns [26–28]. In this paper, we use the NLTV to replace the TV in compound regularization based compressive MR image reconstruction.LetΩ
⊂
R
2, i
,
j
∈
Ω, u
(
x
) be a real function u
:
Ω
→
R, and let w
(
x
,
y
) be a weight function. For a given image u
(
x
), the weighted graph gradient is ∇
N
L
u
(
x
) if defined as the vector of all directional derivatives ∇
N
L
u
(
x
,
·
) at x:
(10)
∇
N
L
u
(
x
,
y
)
:
=
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
∀
y
∈
Ω
.
The directional derivatives apply to all the nodes y since the weight w
(
x
,
y
) is extended to the whole domain Ω
×
Ω. Let us denote vectors such that p
⃑
=
p
(
x
,
y
)
∈
Ω
×
Ω; the nonlocal graph divergence (
di
v
N
L
p
⃑
)
:
Ω
×
Ω
→
Ω is defined as the adjoint of the nonlocal gradient:
(11)
(
di
v
N
L
p
⃑
)
(
x
)
:
=
∑
y
∈
Ω
(
p
(
x
,
y
)
-
p
(
y
,
x
)
)
w
(
x
,
y
)
.Due to being analogous to classical TV, thel
1 norm is in general more efficient than the l
2 norm for sparse reconstruction. In this paper, we are interested in NLTV. Based on the above definition, the NLTV is defined as follows:
(12)
∥
∇
N
L
u
∥
1
:
=
∑
x
∈
Ω
|
∇
N
L
u
(
x
)
|
∥
∇
N
L
u
∥
1
:
=
∑
x
∈
Ω
∑
y
∈
Ω
(
u
(
x
)
-
u
(
y
)
)
2
w
(
x
,
y
)
.
The weight function w
(
x
,
y
) denotes how much the difference between pixels x and y is penalized in the images, which is calculated by
(13)
w
(
x
,
y
)
=
1
C
u
exp
⁡
(
-
∥
f
u
(
x
)
-
f
u
(
y
)
∥
2
2
h
2
)
,
where f
u
(
x
) and f
u
(
y
) denote a small patch in image u centering at the coordinates x and y, respectively. C
u
=
∑
y
∈
Ω
w
(
x
,
y
) is the normalizing factor. h is a filtering parameter.


### 3.3. The Description of Proposed Model and Algorithm

According to the compressive MR image reconstruction models described in Section2, the proposed WasNLTV model for compressive MR image reconstruction is
(14)
min
⁡
u
∥
F
u
-
b
∥
2
2
+
α
Ψ
(
Φ
u
)
+
β
∥
∇
N
L
u
∥
1
.
It should be noted that the optimization problem in (14), although convex, is very hard to solve owing to nonsmooth terms and its huge dimensionality. To solve the problem in (14), we use the variable splitting and augmented Lagrangian algorithm following closely the methodology introduced in [31]. The core idea is to introduce a set of new variables per regularizer and then exploit the alternating direction method of multipliers (ADMM) to solve the resulting constrained optimization problems.By introducing an intermediate variable vector(
v
1
,
v
2
,
v
3
), the problem (14) can be transformed into an equivalent one; that is,
(15)
min
⁡
u
,
v
1
,
v
2
,
v
3
∥
v
1
-
b
∥
2
2
+
α
Ψ
(
v
2
)
+
β
∥
v
3
∥
1
s
.
t
.
iiiii
F
u
=
v
1
,
Φ
u
=
v
2
,
∇
N
L
u
=
v
3
.
The optimization problem (15) can be written in a compact form as follows:
(16)
min
⁡
u
,
v
g
(
v
)
s
.
t
.
G
u
+
B
v
=
0
, where (17)
v
≡
(
v
1
,
v
2
,
v
3
)
,
g
(
v
)
=
∥
v
1
-
b
∥
2
2
+
α
Ψ
(
v
2
)
+
β
∥
v
3
∥
1
,
G
=
[
F
Φ
∇
N
L
]
,
B
=
[
-
I
0
0
0
-
I
0
0
0
-
I
]
.
The augmented Lagrangian of problem (16) is
(18)
L
(
u
,
v
,
d
)
=
g
(
v
)
+
μ
2
∥
G
u
+
B
v
-
d
∥
2
2
,
where μ
>
0 is a positive constant, d
≡
(
d
1
,
d
2
,
d
3
), and d
/
μ denotes the Lagrangian multipliers associated to the constraint G
u
+
B
v
=
0. The basic idea of the augmented Lagrangian method is to seek a saddle point of L
(
u
,
v
,
d
), which is also the solution of problem (16). By using ADMM algorithm, we solve the problem (16) by iteratively solving the following problems:
(19)
(
u
k
+
1
,
v
k
+
1
)
=
min
⁡
u
,
v
L
(
u
,
v
,
d
)
,
(20)
d
k
+
1
=
d
k
-
G
u
k
+
1
-
B
v
k
+
1
.It is evident that the minimization problem (19) is still hard to solve efficiently in a direct way, since it involves a nonseparable quadratic term and nondifferentiability terms. To solve this problem, a quite useful ADMM algorithm is employed, which alternatively minimizes one variable while fixing the other variables. By using ADMM, the problem (19) can be solved by the following four subproblems with respect to u and v.(1)
u subproblem: by fixing v and d, the optimization problem (19) to be solved is
(21)
u
k
+
1
=
min
⁡
u
μ
2
∥
F
u
-
v
1
k
-
d
1
k
∥
2
2
+
μ
2
∥
Φ
u
-
v
2
k
-
d
2
k
∥
2
2
iiiiiiiiiiii
+
μ
2
∥
∇
N
L
u
-
v
3
k
-
d
3
k
∥
2
2
.It is clear that problem (21) is a quadratic function. By direct computation, we get the Euler-Lagrange equation for (21):
(22)
μ
F
T
(
F
u
-
v
1
k
-
d
1
k
)
+
μ
Φ
T
(
Φ
u
-
v
2
k
-
d
2
k
)
+
μ
di
v
N
L
(
∇
N
L
u
-
v
3
k
-
d
3
k
)
=
0
.
Therefore, the solution of problem (21) is
(23)
u
k
+
1
=
(
F
T
F
+
Φ
T
Φ
+
Δ
N
L
)
-
1
i
i
i
i
i
i
i
i
i
×
(
(
v
3
k
+
+
d
3
k
+
1
)
F
T
(
v
1
k
+
1
+
d
1
k
+
1
)
+
Φ
T
(
v
2
k
+
1
+
d
2
k
+
1
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
+
di
v
N
L
(
v
3
k
+
+
d
3
k
+
1
)
)
.Due to the computational complexity of NLTV, the same as [27], the NLTV regularization in this paper only runs one time.(2)
v
1 subproblem: by fixing v
2, v
3, u
, and d, the optimization problem (19) to be solved is
(24)
v
1
k
+
1
=
min
⁡
v
1
∥
v
1
-
b
∥
2
2
+
μ
2
∥
F
u
k
+
1
-
v
1
-
d
1
k
+
1
∥
2
2
.Clearly, the problem (24) is a quadratic function; its solution is simply
(25)
v
1
k
+
1
=
(
2
b
+
μ
F
u
k
+
1
-
μ
d
1
k
+
1
)
(
2
+
μ
)
.
(3)
v
2 subproblem: by fixing v
1, v
3, u
, and d, the optimization problem (19) to be solved is
(26)
v
2
k
+
1
=
min
⁡
v
2
α
Ψ
(
v
2
)
+
μ
2
∥
Φ
u
k
+
1
-
v
2
-
d
2
k
+
1
∥
2
2
.The same as problem (24), the problem (26) is a quadratic function and its gradient ∇
v
2 is simplified as ∇
v
2
=
μ
(
v
2
-
Φ
u
k
+
1
+
d
2
k
+
1
)
+
2
α
σ
2
/
π
(
σ
2
+
v
2
)
-
2. The steepest descent method is desirable to use to solve (26) iteratively by applying
(27)
v
2
k
+
1
=
v
2
k
-
η
∇
v
2
.
(4)
v
3 subproblem: by fixing v
1, v
2, u
, and d, the optimization problem (19) to be solved is
(28)
v
3
k
+
1
=
min
⁡
v
3
β
∥
v
3
∥
1
+
μ
2
∥
∇
N
L
u
k
+
1
-
v
3
-
d
3
k
+
1
∥
2
2
.Problem (28) is a l
1 norm regularized optimization problem. Its solution is the well-known soft threshold [32]:
(29)
v
3
k
+
1
=
soft
(
∇
N
L
u
k
+
1
-
d
3
k
+
1
,
β
μ
)
,
where soft
(
y
,
τ
)
=
sign
⁡
(
y
)
max
⁡
{
|
y
|
-
τ
,
0
} denotes the component-wise application of soft-threshold function.In conclusion, the ADMM algorithm for optimization problem (16) is shown in Algorithm 1.Algorithm 1:Pseudocode of WasNLTV based compressive MR image reconstruction.
Initialization: set k
=
0, choose α, β, μ, u
0, v
1
0, v
2
0, and v
3
0
repeat:
compute
u
sub-problem: u
k
+
1
=
min
u
⁡
μ
/
2
∥
G
u
+
B
v
-
d
∥
2
2 based on (23)
for  i
=
1
,
…
,
3
computer
v
i  sub-problem: v
i
k
+
1
=
min
v
i
⁡
g
(
v
i
)
+
μ
/
2
∥
G
u
+
B
v
i
-
d
i
∥
2
2 based on (25), (27) and (29)
end for
update Lagrange multipliers:
d
1
k
+
1
=
d
1
k
-
F
u
k
+
1
+
v
1
k
+
1
d
2
k
+
1
=
d
2
k
-
Φ
u
k
+
1
+
v
2
k
+
1
d
3
k
+
1
=
d
3
k
-
∇
NL
u
k
+
1
+
v
3
k
+
1
update  iteration: k
=
k
+
1
untilthe stopping criterion is satisfied.


## 3.1. Approximate Sparsity

The problems of usingl
0 norm in compressive MR imaging (i.e., the need for a combinatorial search for its minimization and its too high sensibility to noise) are both due to the fact that the l
0 norm of a vector is a discontinuous function of that vector. The same as [29, 30], our idea is to approximate this discontinuous function by a continuous one, named approximate sparsity function, which provides smooth measure of l
0 norm and better sparsity than l
1 regularizer.The approximate sparsity function is defined as(6)
ψ
σ
(
x
)
=
2
π
arctan
(
|
x
|
σ
2
)
,
x
∈
R
,
σ
∈
R
+
.
The parameter σ may be used to control the accuracy with which ψ
σ approximate the Kronecker delta. In mathematical terms, we have
(7)
lim
⁡
σ
→
0
ψ
σ
(
x
)
=
{
1
,
x
≠
0
,
0
,
x
=
0
.
Define the continuous multivariate approximate sparsity function Ψ
σ
(
x
) as
(8)
Ψ
σ
(
x
)
=
∑
i
=
1
m
ψ
σ
(
x
i
)
,
x
∈
R
m
×
1
.
It is clear from (7) that Ψ
σ
(
x
) is an indicator of the number of zero-entries in x for small values of σ. Therefore, l
0 norm can be approximate by
(9)
∥
x
∥
0
≈
Ψ
σ
(
x
)
=
∑
i
=
1
m
ψ
σ
(
x
i
)
.
Note that the larger the value of σ, the smoother the Ψ
σ
(
x
) and the worse the approximation to l
0 norm; the smaller the value of l
0 norm, the closer the behavior of Ψ
σ
(
x
) to l
0 norm.


## 3.2. Nonlocal Total Variation

Although the classical TV is surprisingly efficient for preserving edges, it is well known that TV is not suitable for images with fine structures, details, and textures which are very important to MR images. The NLTV is a variational extension of the nonlocal means filter proposed by Wang et al. [30]. NLTV uses the whole image information instead of using adjacent pixel information to calculate the gradients in regularization term. The NLTV has been proven to be more efficient than TV for improving the signal-to-noise ratio, on preserving not only sharp edges, but also fine details and repetitive patterns [26–28]. In this paper, we use the NLTV to replace the TV in compound regularization based compressive MR image reconstruction.LetΩ
⊂
R
2, i
,
j
∈
Ω, u
(
x
) be a real function u
:
Ω
→
R, and let w
(
x
,
y
) be a weight function. For a given image u
(
x
), the weighted graph gradient is ∇
N
L
u
(
x
) if defined as the vector of all directional derivatives ∇
N
L
u
(
x
,
·
) at x:
(10)
∇
N
L
u
(
x
,
y
)
:
=
(
u
(
y
)
-
u
(
x
)
)
w
(
x
,
y
)
,
∀
y
∈
Ω
.
The directional derivatives apply to all the nodes y since the weight w
(
x
,
y
) is extended to the whole domain Ω
×
Ω. Let us denote vectors such that p
⃑
=
p
(
x
,
y
)
∈
Ω
×
Ω; the nonlocal graph divergence (
di
v
N
L
p
⃑
)
:
Ω
×
Ω
→
Ω is defined as the adjoint of the nonlocal gradient:
(11)
(
di
v
N
L
p
⃑
)
(
x
)
:
=
∑
y
∈
Ω
(
p
(
x
,
y
)
-
p
(
y
,
x
)
)
w
(
x
,
y
)
.Due to being analogous to classical TV, thel
1 norm is in general more efficient than the l
2 norm for sparse reconstruction. In this paper, we are interested in NLTV. Based on the above definition, the NLTV is defined as follows:
(12)
∥
∇
N
L
u
∥
1
:
=
∑
x
∈
Ω
|
∇
N
L
u
(
x
)
|
∥
∇
N
L
u
∥
1
:
=
∑
x
∈
Ω
∑
y
∈
Ω
(
u
(
x
)
-
u
(
y
)
)
2
w
(
x
,
y
)
.
The weight function w
(
x
,
y
) denotes how much the difference between pixels x and y is penalized in the images, which is calculated by
(13)
w
(
x
,
y
)
=
1
C
u
exp
⁡
(
-
∥
f
u
(
x
)
-
f
u
(
y
)
∥
2
2
h
2
)
,
where f
u
(
x
) and f
u
(
y
) denote a small patch in image u centering at the coordinates x and y, respectively. C
u
=
∑
y
∈
Ω
w
(
x
,
y
) is the normalizing factor. h is a filtering parameter.


## 3.3. The Description of Proposed Model and Algorithm

According to the compressive MR image reconstruction models described in Section2, the proposed WasNLTV model for compressive MR image reconstruction is
(14)
min
⁡
u
∥
F
u
-
b
∥
2
2
+
α
Ψ
(
Φ
u
)
+
β
∥
∇
N
L
u
∥
1
.
It should be noted that the optimization problem in (14), although convex, is very hard to solve owing to nonsmooth terms and its huge dimensionality. To solve the problem in (14), we use the variable splitting and augmented Lagrangian algorithm following closely the methodology introduced in [31]. The core idea is to introduce a set of new variables per regularizer and then exploit the alternating direction method of multipliers (ADMM) to solve the resulting constrained optimization problems.By introducing an intermediate variable vector(
v
1
,
v
2
,
v
3
), the problem (14) can be transformed into an equivalent one; that is,
(15)
min
⁡
u
,
v
1
,
v
2
,
v
3
∥
v
1
-
b
∥
2
2
+
α
Ψ
(
v
2
)
+
β
∥
v
3
∥
1
s
.
t
.
iiiii
F
u
=
v
1
,
Φ
u
=
v
2
,
∇
N
L
u
=
v
3
.
The optimization problem (15) can be written in a compact form as follows:
(16)
min
⁡
u
,
v
g
(
v
)
s
.
t
.
G
u
+
B
v
=
0
, where (17)
v
≡
(
v
1
,
v
2
,
v
3
)
,
g
(
v
)
=
∥
v
1
-
b
∥
2
2
+
α
Ψ
(
v
2
)
+
β
∥
v
3
∥
1
,
G
=
[
F
Φ
∇
N
L
]
,
B
=
[
-
I
0
0
0
-
I
0
0
0
-
I
]
.
The augmented Lagrangian of problem (16) is
(18)
L
(
u
,
v
,
d
)
=
g
(
v
)
+
μ
2
∥
G
u
+
B
v
-
d
∥
2
2
,
where μ
>
0 is a positive constant, d
≡
(
d
1
,
d
2
,
d
3
), and d
/
μ denotes the Lagrangian multipliers associated to the constraint G
u
+
B
v
=
0. The basic idea of the augmented Lagrangian method is to seek a saddle point of L
(
u
,
v
,
d
), which is also the solution of problem (16). By using ADMM algorithm, we solve the problem (16) by iteratively solving the following problems:
(19)
(
u
k
+
1
,
v
k
+
1
)
=
min
⁡
u
,
v
L
(
u
,
v
,
d
)
,
(20)
d
k
+
1
=
d
k
-
G
u
k
+
1
-
B
v
k
+
1
.It is evident that the minimization problem (19) is still hard to solve efficiently in a direct way, since it involves a nonseparable quadratic term and nondifferentiability terms. To solve this problem, a quite useful ADMM algorithm is employed, which alternatively minimizes one variable while fixing the other variables. By using ADMM, the problem (19) can be solved by the following four subproblems with respect to u and v.(1)
u subproblem: by fixing v and d, the optimization problem (19) to be solved is
(21)
u
k
+
1
=
min
⁡
u
μ
2
∥
F
u
-
v
1
k
-
d
1
k
∥
2
2
+
μ
2
∥
Φ
u
-
v
2
k
-
d
2
k
∥
2
2
iiiiiiiiiiii
+
μ
2
∥
∇
N
L
u
-
v
3
k
-
d
3
k
∥
2
2
.It is clear that problem (21) is a quadratic function. By direct computation, we get the Euler-Lagrange equation for (21):
(22)
μ
F
T
(
F
u
-
v
1
k
-
d
1
k
)
+
μ
Φ
T
(
Φ
u
-
v
2
k
-
d
2
k
)
+
μ
di
v
N
L
(
∇
N
L
u
-
v
3
k
-
d
3
k
)
=
0
.
Therefore, the solution of problem (21) is
(23)
u
k
+
1
=
(
F
T
F
+
Φ
T
Φ
+
Δ
N
L
)
-
1
i
i
i
i
i
i
i
i
i
×
(
(
v
3
k
+
+
d
3
k
+
1
)
F
T
(
v
1
k
+
1
+
d
1
k
+
1
)
+
Φ
T
(
v
2
k
+
1
+
d
2
k
+
1
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
+
di
v
N
L
(
v
3
k
+
+
d
3
k
+
1
)
)
.Due to the computational complexity of NLTV, the same as [27], the NLTV regularization in this paper only runs one time.(2)
v
1 subproblem: by fixing v
2, v
3, u
, and d, the optimization problem (19) to be solved is
(24)
v
1
k
+
1
=
min
⁡
v
1
∥
v
1
-
b
∥
2
2
+
μ
2
∥
F
u
k
+
1
-
v
1
-
d
1
k
+
1
∥
2
2
.Clearly, the problem (24) is a quadratic function; its solution is simply
(25)
v
1
k
+
1
=
(
2
b
+
μ
F
u
k
+
1
-
μ
d
1
k
+
1
)
(
2
+
μ
)
.
(3)
v
2 subproblem: by fixing v
1, v
3, u
, and d, the optimization problem (19) to be solved is
(26)
v
2
k
+
1
=
min
⁡
v
2
α
Ψ
(
v
2
)
+
μ
2
∥
Φ
u
k
+
1
-
v
2
-
d
2
k
+
1
∥
2
2
.The same as problem (24), the problem (26) is a quadratic function and its gradient ∇
v
2 is simplified as ∇
v
2
=
μ
(
v
2
-
Φ
u
k
+
1
+
d
2
k
+
1
)
+
2
α
σ
2
/
π
(
σ
2
+
v
2
)
-
2. The steepest descent method is desirable to use to solve (26) iteratively by applying
(27)
v
2
k
+
1
=
v
2
k
-
η
∇
v
2
.
(4)
v
3 subproblem: by fixing v
1, v
2, u
, and d, the optimization problem (19) to be solved is
(28)
v
3
k
+
1
=
min
⁡
v
3
β
∥
v
3
∥
1
+
μ
2
∥
∇
N
L
u
k
+
1
-
v
3
-
d
3
k
+
1
∥
2
2
.Problem (28) is a l
1 norm regularized optimization problem. Its solution is the well-known soft threshold [32]:
(29)
v
3
k
+
1
=
soft
(
∇
N
L
u
k
+
1
-
d
3
k
+
1
,
β
μ
)
,
where soft
(
y
,
τ
)
=
sign
⁡
(
y
)
max
⁡
{
|
y
|
-
τ
,
0
} denotes the component-wise application of soft-threshold function.In conclusion, the ADMM algorithm for optimization problem (16) is shown in Algorithm 1.Algorithm 1:Pseudocode of WasNLTV based compressive MR image reconstruction.
Initialization: set k
=
0, choose α, β, μ, u
0, v
1
0, v
2
0, and v
3
0
repeat:
compute
u
sub-problem: u
k
+
1
=
min
u
⁡
μ
/
2
∥
G
u
+
B
v
-
d
∥
2
2 based on (23)
for  i
=
1
,
…
,
3
computer
v
i  sub-problem: v
i
k
+
1
=
min
v
i
⁡
g
(
v
i
)
+
μ
/
2
∥
G
u
+
B
v
i
-
d
i
∥
2
2 based on (25), (27) and (29)
end for
update Lagrange multipliers:
d
1
k
+
1
=
d
1
k
-
F
u
k
+
1
+
v
1
k
+
1
d
2
k
+
1
=
d
2
k
-
Φ
u
k
+
1
+
v
2
k
+
1
d
3
k
+
1
=
d
3
k
-
∇
NL
u
k
+
1
+
v
3
k
+
1
update  iteration: k
=
k
+
1
untilthe stopping criterion is satisfied.


## 4. Experimental Results

In this section, a series of experiments on four 2D MR images (named brain, chest, artery, and cardiac) are implemented to evaluate the proposed and existing methods. Figure1 shows the test images. All experiments are conducted on a PC with an Intel Core i7-3520M, 2.90 GHz CPU, in MATLAB environment. The proposed method (named WasNLTV) is compared with the existing methods including TVCMRI [19], RecPF [8], and FCSA [5]. We evaluate the performance of various methods both visually and qualitatively in signal-to-noise ratio (SNR) and root-mean-square error (RMSE) values. The SNR and RMSE are defined as
(30)
SNR
=
10
lo
g
10
∥
u
-
E
(
u
)
∥
2
2
∥
u
-
u
~
∥
2
2
,
RMSE
=
E
(
(
u
-
u
~
)
2
)
,
where u and u
~ denote the original image and the reconstructed image, respectively, and E
(
·
) is the mean function.2D test MR images.
(a)
Chest
(b)
Artery
(c)
Cardiac
(d)
BrainFor fair comparisons, experiment uses the same observation methods with TVCMRI. In theK-space, we randomly obtain more samples in low frequencies and fewer samples in higher frequencies. This sampling scheme is widely used for compressed MR image reconstructions. Suppose a MR image u has N pixels and the partial Fourier transform F in problem (1) consists of M rows of N
×
N matrix corresponding to the full 2D discrete Fourier transform. The M chosen rows correspond to the sampling measurements b. Therefore, the sampling ratio is defined as M
/
N. In the experiments, the Gaussian white noise generated by σ
n
×
randn
(
M
,
1
) in MATLAB is added, where standard deviation σ
n
=
0.01. The regularization parameters α, β, and μ are set as 0.001, 0.035, and 1, respectively. To be fair to compare the reconstruction MR images of various algorithms, all methods run 50 iterations and the Rice wavelet toolbox is used as the wavelet transform.Table1 summarizes the average reconstruction accuracy obtained by using different methods at different sampling ratios on the set of test images. From Table 1, it can be seen that the proposed WasNLTV method attains the highest SNR (dB) in all cases. Figure 2 plots the SNR values with sampling ratios for different images. It can also be seen that the WasNLTV method achieves the larger improvement of SNR values.Table 1
SNR (dB) results of different methods with different sampling ratios.
Image
Samp. ratio
TVCMRI
RecPF
FCSA
WasNLTV
Chest (220 × 220)
10
8.12
8.15
8.31
8.74
20
13.12
13.28
13.75
14.70
30
17.71
18.15
19.19
20.39
40
20.95
21.35
22.59
23.64
50
24.19
24.49
26.12
26.97
Artery (220 × 220)
10
8.17
8.39
12.40
13.35
20
14.78
15.51
20.35
21.20
30
18.57
19.86
24.61
25.43
40
22.02
23.66
28.35
29.33
50
24.74
26.62
32.08
33.60
Cardiac (192 × 192)
10
3.77
1.06
4.34
5.19
20
13.45
13.80
15.24
16.01
30
16.93
17.38
18.25
18.97
40
17.07
17.67
18.62
21.71
50
21.26
21.62
22.52
23.24
Brain (210 × 210)
10
4.09
4.14
4.84
5.23
20
10.40
10.64
16.18
17.58
30
14.37
14.85
21.33
22.13
40
17.65
18.34
25.17
25.93
50
21.57
22.45
28.85
29.75Performance comparisons (sampling rate versus SNR) on different MR images.
(a)
Chest
(b)
Artery
(c)
Cardiac
(d)
BrainTable2 gives the RMSE results of reconstructed MRI after applying different algorithms. From Table 2, it can be seen that WasNLTV method attains the lowest RMSE in all cases. As is known, the lower the RMSE is, the better the reconstructed image is. That is to say the MR images reconstructed by WasNLTV have the best visual quality.Table 2
RMSE results of different methods with different sampling ratios.
Image
Samp. ratio
TVCMRI
RecPF
FCSA
WasNLTV
Chest (220 × 220)
10
38.13
38.04
35.09
19.81
20
10.90
10.60
9.71
8.67
30
6.02
5.73
4.94
3.52
40
3.72
3.58
3.00
2.20
50
2.22
2.27
1.90
1.22
Artery (220 × 220)
10
43.95
43.81
41.03
14.67
20
10.19
9.72
3.77
3.35
30
5.06
4.46
2.18
1.77
40
3.37
2.88
1.52
1.17
50
2.17
1.88
1.09
0.64
Cardiac (192 × 192)
10
18.95
18.81
15.86
14.77
20
8.12
7.89
5.71
4.28
30
2.85
2.81
2.61
2.24
40
2.16
2.20
2.15
1.74
50
1.29
1.57
1.51
1.09
Brain (210 × 210)
10
40.27
40.11
35.87
29.11
20
16.11
15.71
6.48
5.87
30
9.48
9.05
3.63
3.14
40
6.46
6.05
2.73
2.57
50
3.25
3.07
2.41
2.01To illustrate visual quality, reconstructed compressive MR images obtained using different methods with sampling ratios 20% are shown in Figures3, 4, 5, and 6. For better visual comparison, we zoom in a small patch where the edge and texture are much more abundant. From the figures, it can be observed that the WasNLTV always obtains the best visual effects on all MR images. In particular, the edge of organs and tissues obtained by WasNLTV are much more clear and easy to identify.Reconstructed cardiac MR images from 20% sampling.
(a)
TVCMRI
(b)
Cropped TVCMRI
(c)
RecPF
(d)
Cropped RecPF
(e)
FCSA
(f)
Cropped FCSA
(g)
WasNLTV
(h)
Cropped WasNLTVReconstructed brain MR images from 20% sampling.
(a)
TVCMRI
(b)
Cropped TVCMRI
(c)
RecPF
(d)
Cropped RecPF
(e)
FCSA
(f)
Cropped FCSA
(g)
WasNLTV
(h)
Cropped WasNLTVReconstructed chest MR images from 20% sampling.
(a)
TVCMRI
(b)
Cropped TVCMRI
(c)
RecPF
(d)
Cropped RecPF
(e)
FCSA
(f)
Cropped FCSA
(g)
WasNLTV
(h)
Cropped WasNLTVReconstructed artery MR images from 20% sampling.
(a)
TVCMRI
(b)
Cropped TVCMRI
(c)
RecPF
(d)
Cropped RecPF
(e)
FCSA
(f)
Cropped FCSA
(g)
WasNLTV
(h)
Cropped WasNLTVFigure7 gives the performance comparisons between different methods with sampling ratios 20% in terms of the CPU time over the SNR. In general, the computational complexity of NLTV is much higher than TV. In order to reduce the computational complexity of WasNLTV, in the experiment, we perform the NLTV regularization once in some iterations. Despite the higher computational complexity of WasNLTV, the WasNLTV obtains the best reconstruction results on all MR images by achieving the highest SNR in less CPU time.Performance comparisons (CPU-time versus SNR) on different MR images.
(a)
Chest
(b)
Artery
(c)
Cardiac
(d)
Brain


## 5. Conclusions

In this paper, we propose a new compound regularization based compressive sensing MRI reconstruction model, which exploits the NLTV regularization and wavelet approximate sparsity prior. The approximate sparsity prior is used in compressive MR image reconstruction model instead ofl
0 or l
1 norm, which can produce much sparser results. And the optimization problem is much easier to be solved. Because the NLTV takes advantage of the redundancy and self-similarity in a MR image, it can effectively avoid blocky artifacts caused by traditional TV regularization and keep fine edge of organs and tissues. As for the algorithm, we apply the variable splitting and augmented Lagrangian algorithm to solve the compound regularization minimization problem. Experiments on test images demonstrate that the proposed method leads to high SNR measure and more importantly preserves the details and edges of MR images.



---

*Source: 137616-2014-08-28.xml*