# A CT Reconstruction Algorithm Based on Non-Aliasing Contourlet Transform and Compressive Sensing

Authors: Lu-zhen Deng; Peng Feng; Mian-yi Chen; Peng He; Quang-sang Vo; Biao Wei
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## Abstract

Compressive sensing (CS) theory has great potential for reconstructing CT images from sparse-views projection data. Currently, total variation (TV-) based CT reconstruction method is a hot research point in medical CT field, which uses the gradient operator as the sparse representation approach during the iteration process. However, the images reconstructed by this method often suffer the smoothing problem; to improve the quality of reconstructed images, this paper proposed a hybrid reconstruction method combining TV and non-aliasing Contourlet transform (NACT) and using the Split-Bregman method to solve the optimization problem. Finally, the simulation results show that the proposed algorithm can reconstruct high-quality CT images from few-views projection using less iteration numbers, which is more effective in suppressing noise and artefacts than algebraic reconstruction
technique (ART) and TV-based reconstruction method.


---

## Body

## 1. Introduction

Since computed tomography (CT) [1] technique was born in 1973, CT has been widely applied in medical diagnose, industrial nondestructive detection, and so forth. In medical CT field, how to reconstruct high-quality CT images from few-views or sparse-views data is a significant research problem. Recently, compressive sensing (CS) [2] theory has been applied in CT images reconstruction which makes it possible to reconstruct high-quality images from few-views data. In CS theory, CT images can be sparsely represented in an appropriate domain, such as gradient transform and Wavelet transform, and the quality of CT reconstructed images will be improved by some appropriate sparse representations in CT images reconstruction.Contourlet transform [3] is proposed by Do and Vetterli in 2002, which is a sparse representation for 2D images with some properties such as multiresolution, multiscale, and multidirection. Contourlet transform can also get important smooth contour features of the image with few data, but there is frequency aliasing in Contourlet transform. Sharp frequency localization Contourlet transform [4] is firstly proposed by Lu and Do in 2006 and Feng et al. introduced a detailed explanation and construction in 2009 which is named as non-aliasing Contourlet transform (NACT) [5]. NACT which can eliminate the frequency aliasing in Contourlet transform is more efficient in capturing geometrical structure and can represent image sparser than traditional Contourlet transform.To solve the optimization problem in CT images reconstruction based on CS, Goldstein and Osher proposed Split-Bregman [6] method, which is derived from Bregman [7] iteration and can accelerate iteration convergence and produce better reconstruction results. Split-Bregman method uses an intermediate variable to split L
1 regularization and L
2 regularization into two equations, where L
2 and L
1 regularization equation can be solved by steepest descent method and thresholding algorithm, respectively. Based on Split-Bregman method, Vandeghinste et al. proposed Split-Bregman-based sparse-view CT reconstruction approach [8]. Furthermore, an iterative CT reconstruction is proposed using shearlet-based regularization [9]. Chu et al. proposed multienergy CT reconstruction based on low rank and sparsity with the Split-Bregman method (MLRSS) [10]. Chang et al. proposed a few-view reweighted sparsity hunting (FRESH) method for CT images reconstruction [11].In this paper, we propose a CT reconstruction algorithm based on NACT and compressive sensing which tries to explore the sparse capability of NACT in order to reconstruct high-quality CT images. In the following section, the proposed algorithm will be introduced. In the third section, we will analyze the experimental results and discuss relevant issues. In the last section, Section4, we will conclude the paper.


## 2. Theory and Method

### 2.1. CT Reconstruction Theory Based on Compressive Sensing

Theoretically, the mathematical CT model can be expressed as:(1)
A
f
=
p
,
where A is the system matrix, f is the original image, and p is the projection data. Traditional CT reconstruction algorithms such as filtered backprojection (FBP) [12] and algebraic reconstruction technique (ART) [13] cannot reconstruct high quality CT images with the sparse sampling or limited projection data.In 2006, Candes and Donoho put forward the CS theory which makes it possible to get high quality CT images with sparse projection data. The main idea of CS is that a signal can be reconstructed with far less sampled frequency than required by conventional Nyquist-Shannon sampling frequency, if the image has a sparse/compressible representation in a transform domain.Compressive sensing theory can be expressed by the following equation:(2)
min
⁡
∥
y
∥
0
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,
where Φ is a orthogonal transform, Φ
H is the corresponding inverse transform, and f is the CT image to be reconstructed and has a special relationship with Φ
H; that is, f
=
Φ
H
y. p is the projection data of f through matrix A.Inspired by CS theory, Sidky et al. proposed a total variation (TV-) based CT reconstruction algorithm using gradient operator as the sparse representation [14], in which TV is defined as follows:
(3)
∥
f
∥
TV
=
∑
|
∇
f
|
=
∑
s
,
t
(
∇
x
f
)
2
+
(
∇
y
f
)
2
=
∑
s
,
t
(
f
s
,
t
-
f
s
-
1
,
t
)
2
+
(
f
s
,
t
-
f
s
,
t
-
1
)
2
,
where ∇
f represents gradient operator of an image f.


### 2.2. Non-Aliasing Contourlet Transform

The traditional sparse representation, such as gradient operator and Wavelet transform [15], cannot get ideal sparse representation of CT images. In 2002, Ni et al. proposed Contourlet transform [16] which can utilize intrinsic structure information of image to represent images more efficiently compared with Wavelet transform. However, suffering from frequency aliasing, Contourlet transform does not show good performance in image denoising, fusion, and enhancement. In order to solve this problem, a new multiscale analysis method, named non-aliasing Contourlet transform (NACT) was proposed. NACT consists of non-aliasing pyramidal filter banks (NPFB) and directional filter banks (DFB). NPFB contains two different filter banks: L
0
(
ω
)
,
D
0
(
ω
) and L
1
(
ω
)
,
D
1
(
ω
). L
0
(
ω
) and L
1
(
ω
) mean low-pass filters. D
0
(
ω
) and D
1
(
ω
) mean high-pass filters. The relationships of two different filter banks are as follows:
(4)
D
1
2
(
ω
)
+
L
1
2
(
ω
)
4
=
1
,
D
0
2
(
ω
)
+
L
0
2
(
ω
)
=
1
.We assume thatω
p
,
0 and ω
s
,
0 represent pass-band frequency and stop-band frequency of L
0
(
ω
), respectively. Accordingly, ω
p
,
1 and ω
s
,
1 represent pass-band frequency and stop-band frequency of L
1
(
ω
), respectively. In order to eliminate frequency aliasing, the filter parameters should meet (
1
)  ω
s
,
1
<
π
/
2; (
2
)  (
ω
p
,
0
+
ω
s
,
0
)
/
2
=
π
/
2 and (
ω
p
,
1
+
ω
s
,
1
)
/
2
=
π
/
4; (
3
)  ω
s
,
0
≤
π
-
a and ω
s
,
1
≤
(
π
-
a
)
/
2, where a is the maximum width of mixing ingredients in DFB [3].As a sparse representation approach, NACT integrate NPFB and DFB which can decompose image into multidirection and multiresolution. NPFB decomposes image into an approximation subband and several detail subbands with different resolutions; DFB decomposes the detail subbands into directional subbands. The process of decomposition with 3 levels is shown in Figure1. We will use “9-7” filter and “pkva” directional filter bank [17] in the study.Figure 1
Flowchart of 3 levels of decomposition of NACT.


### 2.3. Split Bregman Method

In CS theory,L
0 norm is the most ideal regularization norm, but it is difficult to solve equations and easily interfered by noise in CT reconstruction, so L
0 norm is commonly replaced by L
1 norm. Then the reconstruction problem depicted by (2) can be converted into(5a)
min
⁡
∥
y
∥
0
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,
(5b)
min
⁡
∥
Φ
f
∥
1
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,where y
=
Φ
f, Φ is the sparse transform which is normally used as Wavelet transform, Curvelet transform, gradient operator, and so forth.Furthermore, (5b) can be converted into
(6)
f
=
arg
min
⁡
f
⁡
∥
Φ
f
∥
1
+
λ
∥
A
f
-
p
∥
2
2
,
where λ is penalty function weight.In order to solve (6), Goldstein and Osher proposed Split-Bregman method [6], using an intermediate variable to split L
1 regularization and L
2 regularization into two equations; L
2 regularization equation can be solved by gradient descent method and L
1 regularization equation can be solved by thresholding algorithm. Split-Bregman method contains the following three iteration steps:Step 1.
(7)
f
k
+
1
=
arg
min
⁡
f
⁡
λ
∥
A
f
-
p
∥
2
2
+
μ
∥
d
k
-
Φ
f
-
b
k
∥
2
2
.Step 2.
(8)
d
k
+
1
=
min
⁡
d
⁡
∥
d
∥
1
+
μ
∥
d
k
-
Φ
f
k
+
1
-
b
k
∥
2
2
.Step 3.
(9)
b
k
+
1
=
b
k
+
(
Φ
f
k
+
1
-
d
k
+
1
)
,
where k is the Split-Bregman iteration index, μ is convergence parameter, and d and b are intermediate variables, with which each subproblem can be solved easily.


### 2.4. Proposed Algorithm

According to aforementioned methods, we propose a CT reconstruction algorithm based on NACT and compressive sensing method which can be defined as a constrained form (10) or an unconstrained form (11) as follows:
(10)
min
⁡
∥
f
∥
TV
+
∥
Φ
f
∥
1
.
s
.
t
.
i
i
i
∥
A
f
-
p
∥
2
2
<
σ
2
(11)
iiiiiiiiiiiiiiiiiii
min
⁡
∥
f
∥
TV
+
∥
Φ
f
∥
1
+
λ
∥
A
f
-
p
∥
2
2
.
Applying the Split-Bregman method to (11), we have the following three iteration steps:Step 1.
(12)
f
k
+
1
=
arg
min
⁡
f
⁡
λ
∥
A
f
-
p
∥
2
2
+
γ
∥
d
k
-
∇
f
-
b
k
∥
2
2
+
μ
∥
d
φ
k
-
Φ
f
-
b
φ
k
∥
2
2
.Step 2.
(13)
d
k
+
1
=
min
⁡
d
⁡
∥
d
∥
1
+
γ
∥
d
k
-
∇
f
k
+
1
-
b
k
∥
2
2
d
φ
k
+
1
=
min
⁡
d
φ
⁡
∥
d
φ
∥
1
+
μ
∥
d
φ
k
-
Φ
f
k
+
1
-
b
φ
k
∥
2
2
.Step 3.
(14)
b
k
+
1
=
b
k
+
(
∇
f
k
+
1
-
d
k
+
1
)
b
φ
k
+
1
=
b
φ
k
+
(
Φ
f
k
+
1
-
d
φ
k
+
1
)
,
where γ is convergence parameter and d
φ and b
φ are intermediate variables.The steepest descent method is applied to solve (12). The derivative of (12) is calculated as follows:
(15)
g
[
n
,
m
+
1
]
=
2
λ
A
m
T
(
A
m
f
m
-
p
m
)
-
2
γ
∇
T
(
d
k
-
∇
f
m
-
b
k
)
-
2
μ
Φ
T
(
d
φ
k
-
Φ
f
m
-
b
φ
k
)
f
m
+
1
=
f
m
+
α
g
[
n
,
m
+
1
]
,
where n denotes the iteration index of the steepest descent method, m
=
2
,
…
,
N
data denotes the projection angles, A
m ismth row vector, and system matrix A includes N
data row vector A
m. Accordingly, N
data row vectors p
m compose the projection-data vector p, α is an appropriate step size. The ART method is used to get initial image of iteration. Equation (13) can be explicitly computed as (16) using the shrinkage operator as follows:
(16)
d
k
+
1
=
shrink
(
∇
f
k
+
1
+
b
k
,
1
λ
)
d
φ
k
+
1
=
shrink
(
Φ
f
k
+
1
+
b
φ
k
,
1
μ
)
.We now describe the iterative steps of the proposed algorithm. The iteration process contains two loops, the outside loop operate ART and the inside loop solve the optimization problem which is constrained by TV and NACT. The outside loop is labeled byn and the inside loop is labeled by k. The steps comprising each loop are the DATA-step, which enforces consistency with the projection data; the POS-step, which ensures a nonnegative image. We use f
(
ART-DATA
)
[
n
,
m
] to denote the mth DATA-step subiteration with the nth iteration and f
(
ART-POS
)
[
n
] to denote the POS-step with the nth iteration in the outside loop. We use f
(
NACTTV-DATA
)
[
k
,
m
] to denote the mth DATA-step subiteration with the kth iteration and f
(
NACTTV-POS
)
[
k
] to denote the POS-step with the kth iteration in the inside loop. The steps of the algorithm are as follows:(A)
initialization:(17)
n
=
1
,
f
(
ART-DATA
)
[
n
,
1
]
=
0
;
(B)
data projection iteration, form
=
2
,
…
,
N
data:
(18)
f
(
ART-DATA
)
[
n
,
m
]
=
f
(
ART-DATA
)
[
n
,
m
-
1
]
+
A
m
p
m
-
A
m
·
f
(
ART-DATA
)
[
n
,
m
-
1
]
A
m
·
A
m
;
(C)
positivity constraint:(19)
(
f
i
,
j
)
(
ART-POS
)
[
n
]
=
{
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
≥
0
,
0
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
<
0
;
(D)
initialization of Split-Bregman:(20)
k
=
1
,
d
(
n
)
=
∥
f
(
ART-DATA
)
[
n
,
1
]
-
f
(
ART-POS
)
[
n
]
∥
2
,
f
(
NACTTV-DATA
)
[
k
,
1
]
=
f
(
ART-POS
)
[
n
]
d
x
k
=
∇
x
f
(
ART-POS
)
[
n
]
,
d
y
k
=
∇
y
f
(
ART-POS
)
[
n
]
,
d
φ
k
=
Φ
f
(
ART-POS
)
[
n
]
,
b
x
k
=
b
y
k
=
b
φ
k
=
0
;
(E)
iteration form
=
2
,
…
,
N
data:
(21)
d
p
=
A
m
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
-
p
m
-
1
g
[
k
,
m
-
1
]
=
2
λ
A
m
d
p
-
2
γ
∇
x
T
(
d
x
k
-
∇
x
f
-
b
x
k
)
-
2
γ
∇
y
T
(
d
y
k
-
∇
y
f
-
b
y
k
)
-
2
μ
Φ
T
(
d
φ
k
-
Φ
f
-
b
φ
k
)
|
f
=
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
,
g
^
[
k
,
m
-
1
]
=
g
[
k
,
m
-
1
]
|
g
[
k
,
m
-
1
]
|
,
f
(
NACTTV-DATA
)
[
k
,
m
]
=
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
-
a
d
(
n
)
g
^
[
k
,
m
-
1
]
;
(F)
positivity constraint:(22)
(
f
i
,
j
)
(
NACTTV-POS
)
[
k
+
1
]
=
{
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
≥
0
,
0
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
<
0
;
(G)
updated
x, d
y, d
φ, b
x, b
y, b
φ, increase k, and return to step (E) until k
=
K
NACTTV as follows:
(23)
d
x
k
+
1
=
shrink
(
∇
x
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
x
k
,
1
λ
)
b
y
k
+
1
=
shrink
(
∇
y
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
y
k
,
1
λ
)
d
φ
k
+
1
=
shrink
(
Φ
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
φ
k
,
1
μ
)
b
x
k
+
1
=
b
x
k
+
(
∇
x
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
x
k
+
1
)
b
y
k
+
1
=
b
y
k
+
(
∇
y
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
y
k
+
1
)
b
φ
k
+
1
=
b
φ
k
+
(
Φ
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
φ
k
+
1
)
;
(H)
initialize next loop:(24)
f
(
ART-DATA
)
[
n
+
1,1
]
=
f
(
NACTTV-POS
)
[
K
NACTTV
,
1
]
;increase n and return to step (
B
). The iteration is stopped when ∥
A
f
-
p
∥
2
2
<
σ
2. In our study, we selected λ
=
1000, γ
=
30, μ
=
30, a
=
0.2 and K
NACTTV
=
10, which can strike a good balance in the steepest descent and generate good reconstruction results in the experiments.


## 2.1. CT Reconstruction Theory Based on Compressive Sensing

Theoretically, the mathematical CT model can be expressed as:(1)
A
f
=
p
,
where A is the system matrix, f is the original image, and p is the projection data. Traditional CT reconstruction algorithms such as filtered backprojection (FBP) [12] and algebraic reconstruction technique (ART) [13] cannot reconstruct high quality CT images with the sparse sampling or limited projection data.In 2006, Candes and Donoho put forward the CS theory which makes it possible to get high quality CT images with sparse projection data. The main idea of CS is that a signal can be reconstructed with far less sampled frequency than required by conventional Nyquist-Shannon sampling frequency, if the image has a sparse/compressible representation in a transform domain.Compressive sensing theory can be expressed by the following equation:(2)
min
⁡
∥
y
∥
0
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,
where Φ is a orthogonal transform, Φ
H is the corresponding inverse transform, and f is the CT image to be reconstructed and has a special relationship with Φ
H; that is, f
=
Φ
H
y. p is the projection data of f through matrix A.Inspired by CS theory, Sidky et al. proposed a total variation (TV-) based CT reconstruction algorithm using gradient operator as the sparse representation [14], in which TV is defined as follows:
(3)
∥
f
∥
TV
=
∑
|
∇
f
|
=
∑
s
,
t
(
∇
x
f
)
2
+
(
∇
y
f
)
2
=
∑
s
,
t
(
f
s
,
t
-
f
s
-
1
,
t
)
2
+
(
f
s
,
t
-
f
s
,
t
-
1
)
2
,
where ∇
f represents gradient operator of an image f.


## 2.2. Non-Aliasing Contourlet Transform

The traditional sparse representation, such as gradient operator and Wavelet transform [15], cannot get ideal sparse representation of CT images. In 2002, Ni et al. proposed Contourlet transform [16] which can utilize intrinsic structure information of image to represent images more efficiently compared with Wavelet transform. However, suffering from frequency aliasing, Contourlet transform does not show good performance in image denoising, fusion, and enhancement. In order to solve this problem, a new multiscale analysis method, named non-aliasing Contourlet transform (NACT) was proposed. NACT consists of non-aliasing pyramidal filter banks (NPFB) and directional filter banks (DFB). NPFB contains two different filter banks: L
0
(
ω
)
,
D
0
(
ω
) and L
1
(
ω
)
,
D
1
(
ω
). L
0
(
ω
) and L
1
(
ω
) mean low-pass filters. D
0
(
ω
) and D
1
(
ω
) mean high-pass filters. The relationships of two different filter banks are as follows:
(4)
D
1
2
(
ω
)
+
L
1
2
(
ω
)
4
=
1
,
D
0
2
(
ω
)
+
L
0
2
(
ω
)
=
1
.We assume thatω
p
,
0 and ω
s
,
0 represent pass-band frequency and stop-band frequency of L
0
(
ω
), respectively. Accordingly, ω
p
,
1 and ω
s
,
1 represent pass-band frequency and stop-band frequency of L
1
(
ω
), respectively. In order to eliminate frequency aliasing, the filter parameters should meet (
1
)  ω
s
,
1
<
π
/
2; (
2
)  (
ω
p
,
0
+
ω
s
,
0
)
/
2
=
π
/
2 and (
ω
p
,
1
+
ω
s
,
1
)
/
2
=
π
/
4; (
3
)  ω
s
,
0
≤
π
-
a and ω
s
,
1
≤
(
π
-
a
)
/
2, where a is the maximum width of mixing ingredients in DFB [3].As a sparse representation approach, NACT integrate NPFB and DFB which can decompose image into multidirection and multiresolution. NPFB decomposes image into an approximation subband and several detail subbands with different resolutions; DFB decomposes the detail subbands into directional subbands. The process of decomposition with 3 levels is shown in Figure1. We will use “9-7” filter and “pkva” directional filter bank [17] in the study.Figure 1
Flowchart of 3 levels of decomposition of NACT.


## 2.3. Split Bregman Method

In CS theory,L
0 norm is the most ideal regularization norm, but it is difficult to solve equations and easily interfered by noise in CT reconstruction, so L
0 norm is commonly replaced by L
1 norm. Then the reconstruction problem depicted by (2) can be converted into(5a)
min
⁡
∥
y
∥
0
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,
(5b)
min
⁡
∥
Φ
f
∥
1
s
.
t
.
p
=
A
f
=
A
Φ
H
y
,where y
=
Φ
f, Φ is the sparse transform which is normally used as Wavelet transform, Curvelet transform, gradient operator, and so forth.Furthermore, (5b) can be converted into
(6)
f
=
arg
min
⁡
f
⁡
∥
Φ
f
∥
1
+
λ
∥
A
f
-
p
∥
2
2
,
where λ is penalty function weight.In order to solve (6), Goldstein and Osher proposed Split-Bregman method [6], using an intermediate variable to split L
1 regularization and L
2 regularization into two equations; L
2 regularization equation can be solved by gradient descent method and L
1 regularization equation can be solved by thresholding algorithm. Split-Bregman method contains the following three iteration steps:Step 1.
(7)
f
k
+
1
=
arg
min
⁡
f
⁡
λ
∥
A
f
-
p
∥
2
2
+
μ
∥
d
k
-
Φ
f
-
b
k
∥
2
2
.Step 2.
(8)
d
k
+
1
=
min
⁡
d
⁡
∥
d
∥
1
+
μ
∥
d
k
-
Φ
f
k
+
1
-
b
k
∥
2
2
.Step 3.
(9)
b
k
+
1
=
b
k
+
(
Φ
f
k
+
1
-
d
k
+
1
)
,
where k is the Split-Bregman iteration index, μ is convergence parameter, and d and b are intermediate variables, with which each subproblem can be solved easily.


## 2.4. Proposed Algorithm

According to aforementioned methods, we propose a CT reconstruction algorithm based on NACT and compressive sensing method which can be defined as a constrained form (10) or an unconstrained form (11) as follows:
(10)
min
⁡
∥
f
∥
TV
+
∥
Φ
f
∥
1
.
s
.
t
.
i
i
i
∥
A
f
-
p
∥
2
2
<
σ
2
(11)
iiiiiiiiiiiiiiiiiii
min
⁡
∥
f
∥
TV
+
∥
Φ
f
∥
1
+
λ
∥
A
f
-
p
∥
2
2
.
Applying the Split-Bregman method to (11), we have the following three iteration steps:Step 1.
(12)
f
k
+
1
=
arg
min
⁡
f
⁡
λ
∥
A
f
-
p
∥
2
2
+
γ
∥
d
k
-
∇
f
-
b
k
∥
2
2
+
μ
∥
d
φ
k
-
Φ
f
-
b
φ
k
∥
2
2
.Step 2.
(13)
d
k
+
1
=
min
⁡
d
⁡
∥
d
∥
1
+
γ
∥
d
k
-
∇
f
k
+
1
-
b
k
∥
2
2
d
φ
k
+
1
=
min
⁡
d
φ
⁡
∥
d
φ
∥
1
+
μ
∥
d
φ
k
-
Φ
f
k
+
1
-
b
φ
k
∥
2
2
.Step 3.
(14)
b
k
+
1
=
b
k
+
(
∇
f
k
+
1
-
d
k
+
1
)
b
φ
k
+
1
=
b
φ
k
+
(
Φ
f
k
+
1
-
d
φ
k
+
1
)
,
where γ is convergence parameter and d
φ and b
φ are intermediate variables.The steepest descent method is applied to solve (12). The derivative of (12) is calculated as follows:
(15)
g
[
n
,
m
+
1
]
=
2
λ
A
m
T
(
A
m
f
m
-
p
m
)
-
2
γ
∇
T
(
d
k
-
∇
f
m
-
b
k
)
-
2
μ
Φ
T
(
d
φ
k
-
Φ
f
m
-
b
φ
k
)
f
m
+
1
=
f
m
+
α
g
[
n
,
m
+
1
]
,
where n denotes the iteration index of the steepest descent method, m
=
2
,
…
,
N
data denotes the projection angles, A
m ismth row vector, and system matrix A includes N
data row vector A
m. Accordingly, N
data row vectors p
m compose the projection-data vector p, α is an appropriate step size. The ART method is used to get initial image of iteration. Equation (13) can be explicitly computed as (16) using the shrinkage operator as follows:
(16)
d
k
+
1
=
shrink
(
∇
f
k
+
1
+
b
k
,
1
λ
)
d
φ
k
+
1
=
shrink
(
Φ
f
k
+
1
+
b
φ
k
,
1
μ
)
.We now describe the iterative steps of the proposed algorithm. The iteration process contains two loops, the outside loop operate ART and the inside loop solve the optimization problem which is constrained by TV and NACT. The outside loop is labeled byn and the inside loop is labeled by k. The steps comprising each loop are the DATA-step, which enforces consistency with the projection data; the POS-step, which ensures a nonnegative image. We use f
(
ART-DATA
)
[
n
,
m
] to denote the mth DATA-step subiteration with the nth iteration and f
(
ART-POS
)
[
n
] to denote the POS-step with the nth iteration in the outside loop. We use f
(
NACTTV-DATA
)
[
k
,
m
] to denote the mth DATA-step subiteration with the kth iteration and f
(
NACTTV-POS
)
[
k
] to denote the POS-step with the kth iteration in the inside loop. The steps of the algorithm are as follows:(A)
initialization:(17)
n
=
1
,
f
(
ART-DATA
)
[
n
,
1
]
=
0
;
(B)
data projection iteration, form
=
2
,
…
,
N
data:
(18)
f
(
ART-DATA
)
[
n
,
m
]
=
f
(
ART-DATA
)
[
n
,
m
-
1
]
+
A
m
p
m
-
A
m
·
f
(
ART-DATA
)
[
n
,
m
-
1
]
A
m
·
A
m
;
(C)
positivity constraint:(19)
(
f
i
,
j
)
(
ART-POS
)
[
n
]
=
{
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
≥
0
,
0
(
f
i
,
j
)
(
ART-DATA
)
[
n
,
N
data
]
<
0
;
(D)
initialization of Split-Bregman:(20)
k
=
1
,
d
(
n
)
=
∥
f
(
ART-DATA
)
[
n
,
1
]
-
f
(
ART-POS
)
[
n
]
∥
2
,
f
(
NACTTV-DATA
)
[
k
,
1
]
=
f
(
ART-POS
)
[
n
]
d
x
k
=
∇
x
f
(
ART-POS
)
[
n
]
,
d
y
k
=
∇
y
f
(
ART-POS
)
[
n
]
,
d
φ
k
=
Φ
f
(
ART-POS
)
[
n
]
,
b
x
k
=
b
y
k
=
b
φ
k
=
0
;
(E)
iteration form
=
2
,
…
,
N
data:
(21)
d
p
=
A
m
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
-
p
m
-
1
g
[
k
,
m
-
1
]
=
2
λ
A
m
d
p
-
2
γ
∇
x
T
(
d
x
k
-
∇
x
f
-
b
x
k
)
-
2
γ
∇
y
T
(
d
y
k
-
∇
y
f
-
b
y
k
)
-
2
μ
Φ
T
(
d
φ
k
-
Φ
f
-
b
φ
k
)
|
f
=
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
,
g
^
[
k
,
m
-
1
]
=
g
[
k
,
m
-
1
]
|
g
[
k
,
m
-
1
]
|
,
f
(
NACTTV-DATA
)
[
k
,
m
]
=
f
(
NACTTV-DATA
)
[
k
,
m
-
1
]
-
a
d
(
n
)
g
^
[
k
,
m
-
1
]
;
(F)
positivity constraint:(22)
(
f
i
,
j
)
(
NACTTV-POS
)
[
k
+
1
]
=
{
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
≥
0
,
0
(
f
i
,
j
)
(
NACTTV-DATA
)
[
k
,
N
data
]
<
0
;
(G)
updated
x, d
y, d
φ, b
x, b
y, b
φ, increase k, and return to step (E) until k
=
K
NACTTV as follows:
(23)
d
x
k
+
1
=
shrink
(
∇
x
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
x
k
,
1
λ
)
b
y
k
+
1
=
shrink
(
∇
y
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
y
k
,
1
λ
)
d
φ
k
+
1
=
shrink
(
Φ
f
(
NACTTV-POS
)
[
k
+
1,1
]
+
b
φ
k
,
1
μ
)
b
x
k
+
1
=
b
x
k
+
(
∇
x
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
x
k
+
1
)
b
y
k
+
1
=
b
y
k
+
(
∇
y
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
y
k
+
1
)
b
φ
k
+
1
=
b
φ
k
+
(
Φ
f
(
NACTTV-POS
)
[
k
+
1,1
]
-
d
φ
k
+
1
)
;
(H)
initialize next loop:(24)
f
(
ART-DATA
)
[
n
+
1,1
]
=
f
(
NACTTV-POS
)
[
K
NACTTV
,
1
]
;increase n and return to step (
B
). The iteration is stopped when ∥
A
f
-
p
∥
2
2
<
σ
2. In our study, we selected λ
=
1000, γ
=
30, μ
=
30, a
=
0.2 and K
NACTTV
=
10, which can strike a good balance in the steepest descent and generate good reconstruction results in the experiments.


## 3. Experimental Results

### 3.1. The Image Quality Evaluation

This paper uses the root mean square errors (RMSE) and universal quality index (UQI) [18] to evaluate the quality of the reconstructed images.RMSE is the most widely applied way to evaluate image quality, and RMSE is defined as(25)
RMSE
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
i
,
j
R
)
2
,
where f
i
,
j is the pixel value of original image and f
i
,
j
R is the pixel value of reconstructed image.Wang and Bovic proposed UQI mode which evaluates images distortion problem including correlation distortion, brightness distortion, and contrast distortion. The value of UQI is between −1 and 1. When the reconstructed image is the same as the original image, the value of UQI is 1. UQI is defined as(26)
UQI
=
4
σ
f
-
,
f
-
R
f
-
×
f
-
R
(
σ
f
2
+
σ
f
R
2
)
[
(
f
-
)
2
+
(
f
-
R
)
2
]
,
where
(27)
f
-
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
f
i
,
j
,
f
-
R
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
f
i
,
j
R
,
σ
f
2
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
-
)
2
,
σ
f
R
2
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
R
-
f
-
R
)
2
,
σ
f
,
f
R
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
-
)
(
f
i
,
j
R
-
f
-
R
)
.


### 3.2. Numerical Simulation

In this section, a head phantom as shown in Figure2 is used to reconstruct and compare by 3 different methods: ART, ART-TV, and our proposed algorithm (SpBr-NACT). The size of phantom image is 200
×
200. We assume that the CT system was viewed as in a typical pencil-beam geometry, and the scanning range was from 1° to 360° with a θ angular increment; projection angles can be indicated as
(28)
θ
i
=
1
+
360
×
(
i
-
1
)
N
view
,
i
=
1,2
,
…
,
N
view
2
θ
i
=
182
+
360
×
(
i
-
N
view
/
2
)
N
view
,
i
=
N
view
2
,
N
view
2
+
1
,
…
,
N
view
.Figure 2
Head phantom.In the simulation, we reconstruct the head phantom from noise-free and noisy projection data. To obtain noisy projection data, we add 10 dB Gaussian noise into noise-free projection data. Projection numberN
view is 60, and iteration numbers for all reconstruction algorithms are 50. The reconstructed images are shown in Figure 3 and the profile of line 140 in different reconstructed images is plotted in Figure 4.The reconstructed images using three different reconstruction algorithms from the noise-free and noisy data. Top row is for noise-free data and bottom row is for noisy data. (a) and (d) are reconstructed by ART, (b) and (e) are reconstructed by ART-TV, and (c) and (f) are reconstructed by SpBr-NACT method.
(a)
(b)
(c)
(d)
(e)
(f)The profile of line 140 in different reconstructed images. Left column is for noise-free date and right column is for noisy data. (a) and (b) are for ART method; (c) and (d) are for ART-TV method; (e) and (f) are for SpBr-NACT method.
(a)
(b)
(c)
(d)
(e)
(f)From Figures3 and 4, we can see that the reconstructed images using ART and ART-TV methods contain a lot of noise and artifacts, while the reconstructed images using SpBr-NACT method contain less noise and artifacts and have clearer edges.Table1 lists all the RMSE and UQI calculated from reconstructed images. It is obvious that the RMSE of reconstructed images using SpBr-NACT method is much smaller than that of reconstructed images using ART and ART-TV methods; the UQI is much bigger. Thus SpBr-NACT method can reconstruct higher quality images.Table 1
RMSE and UQI of reconstructed images using three different algorithms.
RMSE
UQI
Methods
ART
ART-TV
SpBr-NACT
ART
ART-TV
SpBr-NACT
Noisy-free
0.0502
0.0321
0.0196
0.9869
0.9947
0.9980
Noisy
0.0554
0.0406
0.0318
0.9839
0.9914
0.9948Figure5 plots the change of RMSE and UQI with respect to iteration number. Figure 6 plots the change of RMSE and UQI with respect to projection number N
view. In both figures, ART, ART-TV, and the proposed SpBr-NACT approach are used to reconstruct images from noise-free and noisy data. The blue-solid line is for ART, the green-dashed line is for ART-TV, and the red dashed line is SpBr-NACT. For Figure 5, the projection number is fixed and N
view is 60. For figure 6, the iteration number is fixed and equals 50. From both Figures, it is easy to find that with the increase of projection number or iteration number, SpBr-NACT approach can always get the minimum RMSE and maximum UQI which means that the quality of reconstructed images with SpBr-NACT is better than those with ART and ART-TV. And also we see from Figure 5 when the iteration number is relatively small, 3 methods that almost have the same RMSE and UQI, which implies that our proposed method has no advantage if the iteration step does not converge.The relationship of RMSE and UQI with respect to iteration number. (a) and (b) are the RMSE and UQI of reconstructed images from noise-free data with ART, ART-TV, and SpBr-NACT method, respectively. (c) and (d) are the RMSE and UQI of reconstructed images from noisy data with ART, ART-TV, and SpBr-NACT method, respectively.
(a)
(b)
(c)
(d)The relationship of RMSE and UQI with respect to projection numberN
view. (a) and (b) are the RMSE and UQI of reconstructed images from noise-free data with ART, ART-TV, and SpBr-NACT method, respectively. (c) and (d) are the RMSE and UQI of reconstructed images from noisy data with ART, ART-TV, and SpBr-NACT method, respectively.
(a)
(b)
(c)
(d)


## 3.1. The Image Quality Evaluation

This paper uses the root mean square errors (RMSE) and universal quality index (UQI) [18] to evaluate the quality of the reconstructed images.RMSE is the most widely applied way to evaluate image quality, and RMSE is defined as(25)
RMSE
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
i
,
j
R
)
2
,
where f
i
,
j is the pixel value of original image and f
i
,
j
R is the pixel value of reconstructed image.Wang and Bovic proposed UQI mode which evaluates images distortion problem including correlation distortion, brightness distortion, and contrast distortion. The value of UQI is between −1 and 1. When the reconstructed image is the same as the original image, the value of UQI is 1. UQI is defined as(26)
UQI
=
4
σ
f
-
,
f
-
R
f
-
×
f
-
R
(
σ
f
2
+
σ
f
R
2
)
[
(
f
-
)
2
+
(
f
-
R
)
2
]
,
where
(27)
f
-
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
f
i
,
j
,
f
-
R
=
1
M
×
N
∑
0
≤
i
<
N
∑
0
≤
j
<
M
f
i
,
j
R
,
σ
f
2
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
-
)
2
,
σ
f
R
2
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
R
-
f
-
R
)
2
,
σ
f
,
f
R
=
1
M
×
N
-
1
∑
0
≤
i
<
N
∑
0
≤
j
<
M
(
f
i
,
j
-
f
-
)
(
f
i
,
j
R
-
f
-
R
)
.


## 3.2. Numerical Simulation

In this section, a head phantom as shown in Figure2 is used to reconstruct and compare by 3 different methods: ART, ART-TV, and our proposed algorithm (SpBr-NACT). The size of phantom image is 200
×
200. We assume that the CT system was viewed as in a typical pencil-beam geometry, and the scanning range was from 1° to 360° with a θ angular increment; projection angles can be indicated as
(28)
θ
i
=
1
+
360
×
(
i
-
1
)
N
view
,
i
=
1,2
,
…
,
N
view
2
θ
i
=
182
+
360
×
(
i
-
N
view
/
2
)
N
view
,
i
=
N
view
2
,
N
view
2
+
1
,
…
,
N
view
.Figure 2
Head phantom.In the simulation, we reconstruct the head phantom from noise-free and noisy projection data. To obtain noisy projection data, we add 10 dB Gaussian noise into noise-free projection data. Projection numberN
view is 60, and iteration numbers for all reconstruction algorithms are 50. The reconstructed images are shown in Figure 3 and the profile of line 140 in different reconstructed images is plotted in Figure 4.The reconstructed images using three different reconstruction algorithms from the noise-free and noisy data. Top row is for noise-free data and bottom row is for noisy data. (a) and (d) are reconstructed by ART, (b) and (e) are reconstructed by ART-TV, and (c) and (f) are reconstructed by SpBr-NACT method.
(a)
(b)
(c)
(d)
(e)
(f)The profile of line 140 in different reconstructed images. Left column is for noise-free date and right column is for noisy data. (a) and (b) are for ART method; (c) and (d) are for ART-TV method; (e) and (f) are for SpBr-NACT method.
(a)
(b)
(c)
(d)
(e)
(f)From Figures3 and 4, we can see that the reconstructed images using ART and ART-TV methods contain a lot of noise and artifacts, while the reconstructed images using SpBr-NACT method contain less noise and artifacts and have clearer edges.Table1 lists all the RMSE and UQI calculated from reconstructed images. It is obvious that the RMSE of reconstructed images using SpBr-NACT method is much smaller than that of reconstructed images using ART and ART-TV methods; the UQI is much bigger. Thus SpBr-NACT method can reconstruct higher quality images.Table 1
RMSE and UQI of reconstructed images using three different algorithms.
RMSE
UQI
Methods
ART
ART-TV
SpBr-NACT
ART
ART-TV
SpBr-NACT
Noisy-free
0.0502
0.0321
0.0196
0.9869
0.9947
0.9980
Noisy
0.0554
0.0406
0.0318
0.9839
0.9914
0.9948Figure5 plots the change of RMSE and UQI with respect to iteration number. Figure 6 plots the change of RMSE and UQI with respect to projection number N
view. In both figures, ART, ART-TV, and the proposed SpBr-NACT approach are used to reconstruct images from noise-free and noisy data. The blue-solid line is for ART, the green-dashed line is for ART-TV, and the red dashed line is SpBr-NACT. For Figure 5, the projection number is fixed and N
view is 60. For figure 6, the iteration number is fixed and equals 50. From both Figures, it is easy to find that with the increase of projection number or iteration number, SpBr-NACT approach can always get the minimum RMSE and maximum UQI which means that the quality of reconstructed images with SpBr-NACT is better than those with ART and ART-TV. And also we see from Figure 5 when the iteration number is relatively small, 3 methods that almost have the same RMSE and UQI, which implies that our proposed method has no advantage if the iteration step does not converge.The relationship of RMSE and UQI with respect to iteration number. (a) and (b) are the RMSE and UQI of reconstructed images from noise-free data with ART, ART-TV, and SpBr-NACT method, respectively. (c) and (d) are the RMSE and UQI of reconstructed images from noisy data with ART, ART-TV, and SpBr-NACT method, respectively.
(a)
(b)
(c)
(d)The relationship of RMSE and UQI with respect to projection numberN
view. (a) and (b) are the RMSE and UQI of reconstructed images from noise-free data with ART, ART-TV, and SpBr-NACT method, respectively. (c) and (d) are the RMSE and UQI of reconstructed images from noisy data with ART, ART-TV, and SpBr-NACT method, respectively.
(a)
(b)
(c)
(d)


## 4. Conclusion

In this study, we proposed a CT reconstruction algorithm based on NACT and compressive sensing. The experimental results demonstrate that the proposed method can reconstruct high-quality images from few-views data and has a potential for reducing the radiation dose in clinical application. In the further research, we will try to explore more directional information from NACT so as to improve the performance of SpBr-NACT algorithm, especially when the projection number is far more below what we setup in the current experiment.



---

*Source: 753615-2014-06-30.xml*