# A Model of Regularization Parameter Determination in Low-Dose X-Ray CT Reconstruction Based on Dictionary Learning

Authors: Cheng Zhang; Tao Zhang; Jian Zheng; Ming Li; Yanfei Lu; Jiali You; Yihui Guan
Journal: Computational and Mathematical Methods in Medicine
 (2015)

---

## Abstract

In recent years, X-ray computed tomography (CT) is becoming widely used to reveal patient’s anatomical information. However, the side effect of radiation, relating to genetic or cancerous diseases, has caused great public concern. The problem is how to minimize radiation dose significantly while maintaining image quality. As a practical application of compressed sensing theory, one category of methods takes total variation (TV) minimization as the sparse constraint, which makes it possible and effective to get a reconstruction image of high quality in the undersampling situation. On the other hand, a preliminary attempt of low-dose CT reconstruction based on dictionary learning seems to be another effective choice. But some critical parameters, such as the regularization parameter, cannot be determined by detecting datasets. In this paper, we propose a reweighted objective function that contributes to a numerical calculation model of the regularization parameter. A number of experiments demonstrate that this strategy performs well with better reconstruction images and saving of a large amount of time.


---

## Body

## 1. Introduction

Nowadays, X-ray computed tomography (CT) is still an important part of biomedical imaging technologies for the reason that the reconstructed image is of high spatial resolution and quality. Nevertheless, it confirms that an overdose of radiation possibly increases the risk of genetic or cancerous diseases, making it urgent to develop creative and effective reconstruction techniques to fit low-dose CT scanning protocol. Obviously, the X-ray flux cannot be reduced much since the signal-to-noise ratio (SNR) of measured data declines with the reduction of dose. Another approach is to decrease the number of projection angles, which will lead to incomplete few-view data. In this case, analytic-based algorithms like FDK [1], which are derived from a continuous imaging model and in need of dense sampled projections, are sensitive to insufficient projection data and arrive at a terrible result. However, algebraic algorithms like the simultaneous algebraic reconstruction technique (SART) [2] solved the problem better by transforming it to a series of linear equations.Recently, Candes et al. [3, 4] have made compressed sensing theory popular in information theory field. This theory indicates that a variety of signals can be represented sparsely in a certain transform domain. Therefore, original signal can be recovered accurately by far fewer samples while there is no need to follow the Shannon/Nyquist sampling theorem. A principle called restricted isometry property (RIP) guarantees the perfect recovery of any sparse signal [5]. This novel theory has been applied to many regions, like information technology [6], signal and image processing [7], inverse filtering [8], and so on. It is said that the data acquisition process with compression is good for enhancing image quality because this method can increase imaging speed and suppress the artifacts caused by patients’ movement [9]. For these benefits, many compressed sensing based algorithms are created to deal with few-view CT reconstruction problem. One major group is based on the total variation (TV), which takes the TV of the image as the sparse constraint. The image is determined by minimizing the TV term with the constraints of the linear projection equations. Sidky and Pan presented an improved TV-based algorithm named adaptive steepest descent projection onto convex sets (ASD-POCS) in circular cone-beam framework [10]. Another similar method called gradient projection Barzilari Borwein (GPBB) has a faster convergence speed [11]. Besides the TV minimization algorithms, dictionary learning is also helpful to sparse representation. During the reconstruction process, the image is divided into many overlapped patches, represented sparsely by overcomplete elements of a particular dictionary. Xu et al. combined statistical iterative reconstruction (SIR) with dictionary learning and got a better reconstruction result than TV-based methods in the low-dose CT condition [12]. According to Xu’s paper, this method is robust to noise and obtains a better reconstructed image with more details than the TV-based methods do. Naturally, there are some parameters relevant to the final result. Some of them, like the sparse level, the scale of the dictionary, and so on, have less change due to different scanning data and then can be empirically selected. However, there is a special parameter changing according to the phantom, the scanning protocol, the noise level, and other factors. This parameter plays an important role in the reconstruction program to balance the data fidelity term and the regularization term while determining its value is time consuming with many attempts. Hence, there is no doubt that providing a model to select a proper value of this parameter according to the scanning data is essential for the algorithm based on dictionary learning, which leads to better result and time saving.This paper is organized as follows. In Section2, the problem of low-dose CT reconstruction is stated and the algorithm based on dictionary learning is reviewed. In Section 3, the model of regularization parameter determination is proposed by function fitting method. In Section 4, a series of experiments are performed and corresponding discussions are given. Finally, there is the conclusion at the end of this paper.


## 2. Notation and Problem Description

### 2.1. Background and Notation Interpretation

According to previous work by Xu et al. [12], SIR is united with dictionary learning to derive the algorithm. SIR assumes that the measured data can be regarded as the Poisson distribution(1)
y
i
~
Poisson
b
i
e
-
l
i
+
r
i
,
i
=
1
,
…
,
I
,where b
=
b
1
,
b
2
,
…
,
b
I
T
∈
R
I
×
1 is the entrance X-ray intensity, y
=
y
1
,
y
2
,
…
,
y
I
T
∈
R
I
×
1 is the exit X-ray intensity, l
=
l
1
,
l
2
,
…
,
l
I
T
∈
R
I
×
1 is the integral of the linear attenuation coefficient with l
i
=
A
μ
i
=
∑
j
=
1
N
2
a
i
j
μ
j, A
=
a
i
j
∈
R
I
×
N
2 is the system matrix, the reconstructed image μ
=
μ
1
,
μ
2
,
…
,
μ
N
2
T is a linear attenuation coefficient distribution, which transforms the initial image of N
×
N pixels to a vector μ
∈
R
N
2
×
1, and r
i represents the read-out noise.The objective function of SIR is as(2)
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
R
μ
,where φ
μ
=
∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
^
i
2 is the data fidelity term, l
^
=
l
^
1
,
l
^
2
,
…
,
l
^
I
T
∈
R
I
×
1 is the measured data of l calculated by l
^
i
=
ln
⁡
b
i
/
y
i
-
r
i, ω
i
=
y
i
-
r
i
2
/
y
i is the statistical weight, and R
μ is the regularization term.The regularization term usually contains prior information of the image, like sparse constraint. When the sparse representation is acquired by dictionary learning theory, we can replaceR
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
+
∑
s
ν
s
α
s
0 in the objective function. Therefore, the reconstruction problem is equivalent to the following minimization:(3)
min
μ
,
α
,
D
⁡
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
∑
s
E
s
μ
-
D
α
s
2
2
+
∑
s
ν
s
α
s
0
,where E
s
=
e
n
j
s
∈
R
N
o
2
×
N
2 is an operator to extract patches with N
0
×
N
0 pixels from the image, D
=
d
1
,
d
2
,
…
,
d
K
∈
R
N
0
2
×
K is the training dictionary whose column d
k
∈
R
N
0
2
×
1 is called an atom of the same size of a patch, α
s
∈
R
K
×
1 has few nonzero entries as a sparse representation of patches by the dictionary basis D, and the variables λ and ν
s are regularization parameters. In this optimization problem, μ, α, and D are all unknown; hence, a practical plan of minimizing the object function is an alternating minimization scheme. The plan divides the primary problem into two recursive steps: update of the dictionary model and update of the image. The final result is acquired by operating the two steps alternately until reaching a stopping criterion.


### 2.2. Update of the Dictionary Model

During this procedure, the imageμ is supposed to be fixed, meaning that the data fidelity term is a constant. The optimization problem is simplified to the one as(4)
min
D
,
α
⁡
∑
s
E
s
μ
t
-
D
α
s
2
2
+
∑
x
ν
s
α
s
0
,where μ
t is an intermediate image of the last updating step. In the adaptive dictionary based statistical iterative reconstruction (ADSIR), the dictionary is defined dynamically based on the unknown image while the dictionary in the global dictionary based statistical iterative reconstruction (GDSIR) is predefined beforehand [12]. Previous researches have proved that the K-SVD algorithm performs well at training the dictionary [13]. Once the dictionary is determined, the OMP algorithm is used to update the sparse coding [14] with a predetermined sparse level, instead of solving the l
0-norm problem as (4) directly.


### 2.3. Update of the Image

While updating the image, the dictionary and sparse coding remain invariable. In other words, the problem transforms to the form as(5)
min
μ
⁡
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
∑
s
E
s
μ
-
D
α
s
2
2
,where λ is the regularization parameter balancing the data fidelity term ∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
^
i
2 and the regularization term ∑
s
E
s
μ
-
D
α
s
2
2. The regularization term is already a separable quadratic function. By replacing the data fidelity term with a separable paraboloid surrogate [15], the optimization can be iteratively solved by(6)
μ
j
t
+
1
=
μ
j
t
-
∑
i
=
1
I
a
i
j
ω
i
A
μ
t
i
-
l
^
i
+
2
λ
∑
s
∑
n
=
1
N
0
2
e
n
j
s
E
s
μ
t
n
-
D
α
s
n
∑
i
=
1
I
a
i
j
ω
i
∑
k
=
1
N
2
a
i
k
+
2
λ
∑
s
∑
n
=
1
N
0
2
e
n
j
s
∑
k
=
1
N
2
e
n
k
s
+
j
=
1,2
,
…
,
N
2
.


## 2.1. Background and Notation Interpretation

According to previous work by Xu et al. [12], SIR is united with dictionary learning to derive the algorithm. SIR assumes that the measured data can be regarded as the Poisson distribution(1)
y
i
~
Poisson
b
i
e
-
l
i
+
r
i
,
i
=
1
,
…
,
I
,where b
=
b
1
,
b
2
,
…
,
b
I
T
∈
R
I
×
1 is the entrance X-ray intensity, y
=
y
1
,
y
2
,
…
,
y
I
T
∈
R
I
×
1 is the exit X-ray intensity, l
=
l
1
,
l
2
,
…
,
l
I
T
∈
R
I
×
1 is the integral of the linear attenuation coefficient with l
i
=
A
μ
i
=
∑
j
=
1
N
2
a
i
j
μ
j, A
=
a
i
j
∈
R
I
×
N
2 is the system matrix, the reconstructed image μ
=
μ
1
,
μ
2
,
…
,
μ
N
2
T is a linear attenuation coefficient distribution, which transforms the initial image of N
×
N pixels to a vector μ
∈
R
N
2
×
1, and r
i represents the read-out noise.The objective function of SIR is as(2)
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
R
μ
,where φ
μ
=
∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
^
i
2 is the data fidelity term, l
^
=
l
^
1
,
l
^
2
,
…
,
l
^
I
T
∈
R
I
×
1 is the measured data of l calculated by l
^
i
=
ln
⁡
b
i
/
y
i
-
r
i, ω
i
=
y
i
-
r
i
2
/
y
i is the statistical weight, and R
μ is the regularization term.The regularization term usually contains prior information of the image, like sparse constraint. When the sparse representation is acquired by dictionary learning theory, we can replaceR
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
+
∑
s
ν
s
α
s
0 in the objective function. Therefore, the reconstruction problem is equivalent to the following minimization:(3)
min
μ
,
α
,
D
⁡
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
∑
s
E
s
μ
-
D
α
s
2
2
+
∑
s
ν
s
α
s
0
,where E
s
=
e
n
j
s
∈
R
N
o
2
×
N
2 is an operator to extract patches with N
0
×
N
0 pixels from the image, D
=
d
1
,
d
2
,
…
,
d
K
∈
R
N
0
2
×
K is the training dictionary whose column d
k
∈
R
N
0
2
×
1 is called an atom of the same size of a patch, α
s
∈
R
K
×
1 has few nonzero entries as a sparse representation of patches by the dictionary basis D, and the variables λ and ν
s are regularization parameters. In this optimization problem, μ, α, and D are all unknown; hence, a practical plan of minimizing the object function is an alternating minimization scheme. The plan divides the primary problem into two recursive steps: update of the dictionary model and update of the image. The final result is acquired by operating the two steps alternately until reaching a stopping criterion.


## 2.2. Update of the Dictionary Model

During this procedure, the imageμ is supposed to be fixed, meaning that the data fidelity term is a constant. The optimization problem is simplified to the one as(4)
min
D
,
α
⁡
∑
s
E
s
μ
t
-
D
α
s
2
2
+
∑
x
ν
s
α
s
0
,where μ
t is an intermediate image of the last updating step. In the adaptive dictionary based statistical iterative reconstruction (ADSIR), the dictionary is defined dynamically based on the unknown image while the dictionary in the global dictionary based statistical iterative reconstruction (GDSIR) is predefined beforehand [12]. Previous researches have proved that the K-SVD algorithm performs well at training the dictionary [13]. Once the dictionary is determined, the OMP algorithm is used to update the sparse coding [14] with a predetermined sparse level, instead of solving the l
0-norm problem as (4) directly.


## 2.3. Update of the Image

While updating the image, the dictionary and sparse coding remain invariable. In other words, the problem transforms to the form as(5)
min
μ
⁡
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
^
i
2
+
λ
∑
s
E
s
μ
-
D
α
s
2
2
,where λ is the regularization parameter balancing the data fidelity term ∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
^
i
2 and the regularization term ∑
s
E
s
μ
-
D
α
s
2
2. The regularization term is already a separable quadratic function. By replacing the data fidelity term with a separable paraboloid surrogate [15], the optimization can be iteratively solved by(6)
μ
j
t
+
1
=
μ
j
t
-
∑
i
=
1
I
a
i
j
ω
i
A
μ
t
i
-
l
^
i
+
2
λ
∑
s
∑
n
=
1
N
0
2
e
n
j
s
E
s
μ
t
n
-
D
α
s
n
∑
i
=
1
I
a
i
j
ω
i
∑
k
=
1
N
2
a
i
k
+
2
λ
∑
s
∑
n
=
1
N
0
2
e
n
j
s
∑
k
=
1
N
2
e
n
k
s
+
j
=
1,2
,
…
,
N
2
.


## 3. Materials and Methods

### 3.1. Effect of the Regularization Parameter

As mentioned above, the regularization parameterλ is of great importance during the update of image (5). We consider the optimizing problem of the form(7)
μ
λ
=
arg
⁡
min
μ
⁡
φ
μ
+
λ
R
μ
,
φ
μ
=
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
i
2
,
R
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
.If there is μ
λ
1
=
min
μ
⁡
φ
μ
+
λ
1
R
μ, μ
λ
2
=
min
μ
⁡
φ
μ
+
λ
2
R
μ, 0
<
λ
1
<
λ
2, then we can get φ
μ
λ
1
≤
φ
μ
λ
2 and R
μ
λ
1
≥
R
μ
λ
2 by an easy derivation of the unequal relations:(8)
φ
μ
λ
1
+
λ
1
R
μ
λ
1
≤
φ
μ
λ
2
+
λ
1
R
μ
λ
2
,
φ
μ
λ
2
+
λ
2
R
μ
λ
2
≤
φ
μ
λ
1
+
λ
2
R
μ
λ
1
.It shows that a smaller λ makes the data fidelity term smaller and the regularization term bigger, which means that the sparse constraint has less effect on the optimizing process and more noise will appear in the final image. On the other hand, a bigger λ weakens the effect of the data fidelity term, generating a loss of some fine details in the image. For example, λ should be increased to suppress the noise increment in the projection domain since the data fidelity term is proportional to the noise standard deviation. In order to get an optimal result, previous work selects a great many values of λ and picks out the best one by comparing the final images. This testing strategy is of great time consuming, making the algorithm based on dictionary learning not friendly to the reconstruction task.


### 3.2. Morozov’s Principle and the Balancing Principle

The research on the choices of regularization parameters in linear inverse problems appears early in 1998 [16]. The original optimizing function of the inverse problem is like(9)
J
f
,
β
=
1
2
T
f
-
z
Y
2
+
β
f
X
2
.Take Y
=
2, X
=
2 as an example. If the noise level of z is known, one efficient tool for selecting the proper regularization parameter β is the well-known Morozov discrepancy principle [16]. To adaptively determine the regularization parameter, a model function is brought in [17]:(10)
m
β
=
b
+
s
t
+
β
.To find a solution of β, (9) is rewritten as(11)
F
β
=
1
2
T
f
β
-
z
2
2
+
β
f
β
2
2
=
φ
β
+
β
F
′
β
with
φ
β
=
1
2
T
f
β
-
z
2
2
,
F
′
β
=
f
β
2
2
.When β
→
∞, it obviously shows that the solution of the minimization problem is f
β
=
0, and then b
=
1
/
2
z
2
2 is obtained easily. The other two variables s and t can be determined by the equations m
β
k
=
F
β
k, m
′
β
k
=
F
′
β
k. To solve the parameter β iteratively, the balancing principle [18] is introduced as(12)
σ
-
1
φ
β
∗
=
β
∗
F
′
β
∗
,where σ
>
1 controls the relative weight of the two terms. Equation (12) can be written as(13)
F
β
∗
=
σ
F
β
∗
-
β
∗
F
′
β
∗
,which is a fixed point iteration. β
k
+
1 is calculated by the formula(14)
F
β
k
+
1
=
σ
F
β
k
-
β
k
F
′
β
k
=
σ
m
β
k
-
β
k
m
′
β
k
.Although the balancing principle behaves well in the inverse problem model, there is no direct way introducing the method to the dictionary based algorithm. The regularization term has a minimum value greater than zero, which leads to the derivation as follows:(15)
for
R
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
≥
δ
>
0
assume
F
λ
=
∑
i
=
1
I
ω
i
2
A
μ
λ
i
-
l
i
2
+
λ
∑
s
E
s
μ
λ
-
D
α
s
2
2
=
φ
λ
+
λ
F
′
λ
F
λ
=
m
λ
=
b
+
s
t
+
λ
when
λ
⟶
∞
,
b
≥
λ
δ
⟶
∞
.Therefore, the strategy that determining the regularization parameter adaptively accords to the last iterative result is not reasonable. We should look for a selecting strategy which can determine the proper value of the regularization parameter by making an analysis of the projection data.


### 3.3. Weight Modification of the Objective Function

In order to find out an applicable selecting model of the regularization parameter, we reconsider the minimizing problem (5) and the updating formula (6). According to the former work, the data fidelity term is replaced with a separable surrogate [15](16)
Q
μ
;
μ
t
=
∑
i
=
1
I
∑
j
=
1
N
2
α
i
j
ω
i
2
×
a
i
j
α
i
j
μ
j
-
μ
j
t
+
A
μ
t
i
-
l
^
i
2
,where (17)
∑
j
=
1
N
2
α
i
j
=
1
∀
i
,
α
i
j
≥
0
.And one convenient choice is(18)
α
i
j
=
a
i
j
∑
j
=
1
N
2
a
i
j
.By making use of the surrogate, (5) becomes a separable form(19)
μ
=
arg
⁡
min
μ
⁡
Q
μ
;
μ
t
+
λ
R
μ
,where Q
μ
;
μ
t and R
μ both can be represented as the sum of a series of quadratic functions, whose variables are μ
j
j
=
1,2
,
…
,
N
2. After some variable exchanges, (19) can be rewritten as(20)
Q
μ
;
μ
t
=
1
2
∑
j
=
1
N
2
p
j
μ
j
-
c
j
t
2
+
C
1
,
R
μ
=
∑
j
=
1
N
2
q
j
μ
j
-
d
j
t
2
+
C
2
with
p
j
=
∑
i
=
1
I
a
i
j
ω
i
∑
k
=
1
N
2
a
i
k
,
q
j
=
∑
s
∑
n
=
1
N
0
2
e
n
j
s
∑
k
=
1
N
2
e
n
k
s
,
c
j
t
=
μ
j
t
-
∑
i
=
1
I
a
i
j
ω
i
A
μ
t
i
-
l
^
i
p
j
,
d
j
t
=
μ
j
t
-
∑
s
∑
n
=
1
N
0
2
e
n
j
s
E
s
μ
t
n
-
D
α
s
n
q
j
,
μ
j
t
+
1
=
p
j
c
j
t
+
2
λ
q
j
d
j
t
p
j
+
2
λ
q
j
+
.From the above, the image updating formula is just the same as (6); the quadratic term coefficient p
j only depends on the system matrix A and the statistical weight  ω. We make a weight modification on the regularization term(21)
R
~
μ
=
∑
j
=
1
n
r
j
q
j
μ
j
-
d
j
t
2
+
C
2
=
1
2
∑
j
=
1
n
p
j
μ
j
-
d
j
t
2
+
C
2
s
.
t
.
r
j
=
1
2
p
j
q
j
.By eliminating the constant term, the image reconstruction process is equivalent to solving the following optimization problem:(22)
μ
=
arg
⁡
min
μ
⁡
Q
μ
;
μ
t
+
λ
R
~
μ
=
arg
⁡
min
μ
⁡
1
2
∑
j
=
1
N
2
p
j
μ
j
-
c
j
t
2
+
1
2
λ
∑
j
=
1
n
p
j
μ
j
-
d
j
t
2
,which can be solved iteratively by (23)
μ
j
t
+
1
=
c
j
t
+
λ
d
j
t
1
+
λ
+
.As shown in (23), λ determines the relative impact on the updating image μ by the data fidelity term and the regularization term, respectively.


### 3.4. Evaluation Model of Regularization Parameter

Before developing the evaluation model, some discussions about the reconstruction result are displayed firstly. Once a value ofλ is selected randomly, by solving (22) iteratively, it infers that the relative error of the data fidelity term is as(24)
δ
=
∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
i
2
∑
i
=
1
I
ω
i
/
2
l
i
2
.The relative error depends on the phantom image, the noise level, the regularization parameter, and so on. When it comes to the situation that the regularization parameter is infinite as λ
→
∞, (23) and (24) will be like the following form:(25)
μ
j
t
+
1
=
d
j
t
+
,
(26)
δ
λ
→
∞
=
∑
i
=
1
I
ω
i
/
2
A
μ
λ
→
∞
∗
i
-
l
i
2
∑
i
=
1
I
ω
i
/
2
l
i
2
s
.
t
.
λ
⟶
∞
.It is naturally derived that the relative error δ
λ
→
∞ increases with the increment of the noise level in projection domain. In addition, it has been mentioned above (in Section 3.1) that λ should be increased with the noise increment. Therefore, the proper λ has a monotonous relation with the parameter δ
λ
→
∞. Since δ
λ
→
∞ can be easily determined by operating the reconstruction algorithm based on dictionary learning once with λ
→
∞, the proper value of λ can be calculated if a reasonable function as λ
∗
=
f
δ
λ
→
∞ can be found.With the help of a series of tests, the relationship betweenλ
∗ and δ
λ
→
∞ is fitted by a piecewise quadratic function as follows:(27)
λ
∗
=
1.74485
δ
G
2
+
0.58883
δ
G
-
6.88253
if
δ
G
>
1.96
λ
∗
=
-
0.21545
δ
G
2
+
1.08602
δ
G
-
0.32634
i
f
δ
G
≤
1.96
w
i
t
h
δ
G
=
1
0
6
δ
λ
→
∞
.Finally, taking ADSIR as an example, the workflow of the developed algorithm is exhibited in Algorithm1. In addition, the ordered subsets convex (OSC) algorithm [19] is utilized as an acceleration of the convergence.Algorithm 1:Workflow of the developed algorithm.
Determination of
λ
Initializeμ
0
,
D
0
,
α
s
0, and t
=
0.
While the stopping criterion is not satisfied, do
(1) Implement the OSC algorithm for acceleration;
(2) Extract patches from the intermediate imageμ
t;
(3) Update the dictionaryD
t
+
1 by K-SVD algorithm;
(4) Update the sparse codingα
s
t
+
1 by OMP algorithm;
(5) Update the imageμ
t
+
1 by (25) with λ
→
∞
,
t
=
t
+
1;
Output the final imageμ
∗
=
μ
t and relative error δ
λ
→
∞ calculated by (26).
When the iteration is stopped, determineλ by (27).
Image Reconstruction
Initializeμ
0
,
D
0
,
α
s
0, and t
=
0, λ is determined in former step.
While the stopping criterion is not satisfied, do
(1) Implement the OSC algorithm for acceleration;
(2) Extract patches from the intermediate imageμ
t;
(3) Update the dictionaryD
t
+
1 by K-SVD algorithm;
(4) Update the sparse codingα
s
t
+
1 by OMP algorithm;
(5) Update the imageμ
t
+
1 by (23), t
=
t
+
1;
Output the final reconstruction.


## 3.1. Effect of the Regularization Parameter

As mentioned above, the regularization parameterλ is of great importance during the update of image (5). We consider the optimizing problem of the form(7)
μ
λ
=
arg
⁡
min
μ
⁡
φ
μ
+
λ
R
μ
,
φ
μ
=
∑
i
=
1
I
ω
i
2
A
μ
i
-
l
i
2
,
R
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
.If there is μ
λ
1
=
min
μ
⁡
φ
μ
+
λ
1
R
μ, μ
λ
2
=
min
μ
⁡
φ
μ
+
λ
2
R
μ, 0
<
λ
1
<
λ
2, then we can get φ
μ
λ
1
≤
φ
μ
λ
2 and R
μ
λ
1
≥
R
μ
λ
2 by an easy derivation of the unequal relations:(8)
φ
μ
λ
1
+
λ
1
R
μ
λ
1
≤
φ
μ
λ
2
+
λ
1
R
μ
λ
2
,
φ
μ
λ
2
+
λ
2
R
μ
λ
2
≤
φ
μ
λ
1
+
λ
2
R
μ
λ
1
.It shows that a smaller λ makes the data fidelity term smaller and the regularization term bigger, which means that the sparse constraint has less effect on the optimizing process and more noise will appear in the final image. On the other hand, a bigger λ weakens the effect of the data fidelity term, generating a loss of some fine details in the image. For example, λ should be increased to suppress the noise increment in the projection domain since the data fidelity term is proportional to the noise standard deviation. In order to get an optimal result, previous work selects a great many values of λ and picks out the best one by comparing the final images. This testing strategy is of great time consuming, making the algorithm based on dictionary learning not friendly to the reconstruction task.


## 3.2. Morozov’s Principle and the Balancing Principle

The research on the choices of regularization parameters in linear inverse problems appears early in 1998 [16]. The original optimizing function of the inverse problem is like(9)
J
f
,
β
=
1
2
T
f
-
z
Y
2
+
β
f
X
2
.Take Y
=
2, X
=
2 as an example. If the noise level of z is known, one efficient tool for selecting the proper regularization parameter β is the well-known Morozov discrepancy principle [16]. To adaptively determine the regularization parameter, a model function is brought in [17]:(10)
m
β
=
b
+
s
t
+
β
.To find a solution of β, (9) is rewritten as(11)
F
β
=
1
2
T
f
β
-
z
2
2
+
β
f
β
2
2
=
φ
β
+
β
F
′
β
with
φ
β
=
1
2
T
f
β
-
z
2
2
,
F
′
β
=
f
β
2
2
.When β
→
∞, it obviously shows that the solution of the minimization problem is f
β
=
0, and then b
=
1
/
2
z
2
2 is obtained easily. The other two variables s and t can be determined by the equations m
β
k
=
F
β
k, m
′
β
k
=
F
′
β
k. To solve the parameter β iteratively, the balancing principle [18] is introduced as(12)
σ
-
1
φ
β
∗
=
β
∗
F
′
β
∗
,where σ
>
1 controls the relative weight of the two terms. Equation (12) can be written as(13)
F
β
∗
=
σ
F
β
∗
-
β
∗
F
′
β
∗
,which is a fixed point iteration. β
k
+
1 is calculated by the formula(14)
F
β
k
+
1
=
σ
F
β
k
-
β
k
F
′
β
k
=
σ
m
β
k
-
β
k
m
′
β
k
.Although the balancing principle behaves well in the inverse problem model, there is no direct way introducing the method to the dictionary based algorithm. The regularization term has a minimum value greater than zero, which leads to the derivation as follows:(15)
for
R
μ
=
∑
s
E
s
μ
-
D
α
s
2
2
≥
δ
>
0
assume
F
λ
=
∑
i
=
1
I
ω
i
2
A
μ
λ
i
-
l
i
2
+
λ
∑
s
E
s
μ
λ
-
D
α
s
2
2
=
φ
λ
+
λ
F
′
λ
F
λ
=
m
λ
=
b
+
s
t
+
λ
when
λ
⟶
∞
,
b
≥
λ
δ
⟶
∞
.Therefore, the strategy that determining the regularization parameter adaptively accords to the last iterative result is not reasonable. We should look for a selecting strategy which can determine the proper value of the regularization parameter by making an analysis of the projection data.


## 3.3. Weight Modification of the Objective Function

In order to find out an applicable selecting model of the regularization parameter, we reconsider the minimizing problem (5) and the updating formula (6). According to the former work, the data fidelity term is replaced with a separable surrogate [15](16)
Q
μ
;
μ
t
=
∑
i
=
1
I
∑
j
=
1
N
2
α
i
j
ω
i
2
×
a
i
j
α
i
j
μ
j
-
μ
j
t
+
A
μ
t
i
-
l
^
i
2
,where (17)
∑
j
=
1
N
2
α
i
j
=
1
∀
i
,
α
i
j
≥
0
.And one convenient choice is(18)
α
i
j
=
a
i
j
∑
j
=
1
N
2
a
i
j
.By making use of the surrogate, (5) becomes a separable form(19)
μ
=
arg
⁡
min
μ
⁡
Q
μ
;
μ
t
+
λ
R
μ
,where Q
μ
;
μ
t and R
μ both can be represented as the sum of a series of quadratic functions, whose variables are μ
j
j
=
1,2
,
…
,
N
2. After some variable exchanges, (19) can be rewritten as(20)
Q
μ
;
μ
t
=
1
2
∑
j
=
1
N
2
p
j
μ
j
-
c
j
t
2
+
C
1
,
R
μ
=
∑
j
=
1
N
2
q
j
μ
j
-
d
j
t
2
+
C
2
with
p
j
=
∑
i
=
1
I
a
i
j
ω
i
∑
k
=
1
N
2
a
i
k
,
q
j
=
∑
s
∑
n
=
1
N
0
2
e
n
j
s
∑
k
=
1
N
2
e
n
k
s
,
c
j
t
=
μ
j
t
-
∑
i
=
1
I
a
i
j
ω
i
A
μ
t
i
-
l
^
i
p
j
,
d
j
t
=
μ
j
t
-
∑
s
∑
n
=
1
N
0
2
e
n
j
s
E
s
μ
t
n
-
D
α
s
n
q
j
,
μ
j
t
+
1
=
p
j
c
j
t
+
2
λ
q
j
d
j
t
p
j
+
2
λ
q
j
+
.From the above, the image updating formula is just the same as (6); the quadratic term coefficient p
j only depends on the system matrix A and the statistical weight  ω. We make a weight modification on the regularization term(21)
R
~
μ
=
∑
j
=
1
n
r
j
q
j
μ
j
-
d
j
t
2
+
C
2
=
1
2
∑
j
=
1
n
p
j
μ
j
-
d
j
t
2
+
C
2
s
.
t
.
r
j
=
1
2
p
j
q
j
.By eliminating the constant term, the image reconstruction process is equivalent to solving the following optimization problem:(22)
μ
=
arg
⁡
min
μ
⁡
Q
μ
;
μ
t
+
λ
R
~
μ
=
arg
⁡
min
μ
⁡
1
2
∑
j
=
1
N
2
p
j
μ
j
-
c
j
t
2
+
1
2
λ
∑
j
=
1
n
p
j
μ
j
-
d
j
t
2
,which can be solved iteratively by (23)
μ
j
t
+
1
=
c
j
t
+
λ
d
j
t
1
+
λ
+
.As shown in (23), λ determines the relative impact on the updating image μ by the data fidelity term and the regularization term, respectively.


## 3.4. Evaluation Model of Regularization Parameter

Before developing the evaluation model, some discussions about the reconstruction result are displayed firstly. Once a value ofλ is selected randomly, by solving (22) iteratively, it infers that the relative error of the data fidelity term is as(24)
δ
=
∑
i
=
1
I
ω
i
/
2
A
μ
i
-
l
i
2
∑
i
=
1
I
ω
i
/
2
l
i
2
.The relative error depends on the phantom image, the noise level, the regularization parameter, and so on. When it comes to the situation that the regularization parameter is infinite as λ
→
∞, (23) and (24) will be like the following form:(25)
μ
j
t
+
1
=
d
j
t
+
,
(26)
δ
λ
→
∞
=
∑
i
=
1
I
ω
i
/
2
A
μ
λ
→
∞
∗
i
-
l
i
2
∑
i
=
1
I
ω
i
/
2
l
i
2
s
.
t
.
λ
⟶
∞
.It is naturally derived that the relative error δ
λ
→
∞ increases with the increment of the noise level in projection domain. In addition, it has been mentioned above (in Section 3.1) that λ should be increased with the noise increment. Therefore, the proper λ has a monotonous relation with the parameter δ
λ
→
∞. Since δ
λ
→
∞ can be easily determined by operating the reconstruction algorithm based on dictionary learning once with λ
→
∞, the proper value of λ can be calculated if a reasonable function as λ
∗
=
f
δ
λ
→
∞ can be found.With the help of a series of tests, the relationship betweenλ
∗ and δ
λ
→
∞ is fitted by a piecewise quadratic function as follows:(27)
λ
∗
=
1.74485
δ
G
2
+
0.58883
δ
G
-
6.88253
if
δ
G
>
1.96
λ
∗
=
-
0.21545
δ
G
2
+
1.08602
δ
G
-
0.32634
i
f
δ
G
≤
1.96
w
i
t
h
δ
G
=
1
0
6
δ
λ
→
∞
.Finally, taking ADSIR as an example, the workflow of the developed algorithm is exhibited in Algorithm1. In addition, the ordered subsets convex (OSC) algorithm [19] is utilized as an acceleration of the convergence.Algorithm 1:Workflow of the developed algorithm.
Determination of
λ
Initializeμ
0
,
D
0
,
α
s
0, and t
=
0.
While the stopping criterion is not satisfied, do
(1) Implement the OSC algorithm for acceleration;
(2) Extract patches from the intermediate imageμ
t;
(3) Update the dictionaryD
t
+
1 by K-SVD algorithm;
(4) Update the sparse codingα
s
t
+
1 by OMP algorithm;
(5) Update the imageμ
t
+
1 by (25) with λ
→
∞
,
t
=
t
+
1;
Output the final imageμ
∗
=
μ
t and relative error δ
λ
→
∞ calculated by (26).
When the iteration is stopped, determineλ by (27).
Image Reconstruction
Initializeμ
0
,
D
0
,
α
s
0, and t
=
0, λ is determined in former step.
While the stopping criterion is not satisfied, do
(1) Implement the OSC algorithm for acceleration;
(2) Extract patches from the intermediate imageμ
t;
(3) Update the dictionaryD
t
+
1 by K-SVD algorithm;
(4) Update the sparse codingα
s
t
+
1 by OMP algorithm;
(5) Update the imageμ
t
+
1 by (23), t
=
t
+
1;
Output the final reconstruction.


## 4. Experimental Results and Discussion

To make the evaluation of the regularization parameter possible, the developed algorithm improves ADSIR with a modified weight of the regularization term while the weight is adaptive to the data fidelity term. So the proposed algorithm is named adaptive weight regularized ADSIR (AWR-ADSIR). In this section, a series of reconstruction experiments are exhibited to validate that the regularization selecting principle in AWR-ADSIR is practical. The simulation numerical phantoms are Shepp-Logan phantom, human head slice image, and human abdomen slice image. The Shepp-Logan phantom is a numeric phantom with pixels intensities ranging from 0 to 1. The sample images of human head slice and human abdomen slice are the FBP reconstruction results based on full-sampling scanning data, which are obtained from our collaborator. All of these phantom images are of256
×
256 pixels presented in Figure 1. The scanning data are simulated as an undersampling situation with different noise levels. Firstly, different regularization parameters are selected to demonstrate that the one chosen by the algorithm leads to the best reconstruction result. Secondly, by comparing the quality of images reconstructed by diverse algorithms, which are SART, GPBB, ADSIR, and AWR-ADSIR, it can be confirmed that AWR-ADSIR is of remarkable performance among these algorithms. Finally, the selecting principle is used in the GDSIR model, proving that it also works. All the algorithms above are coded in MATLAB and run on a dual-core PC with 3.10 GHz Intel Core i5-2400 and 4 GB RAM.Figure 1
From left to right are Shepp-Logan phantom, human head slice image, and human abdomen slice image. The whole windows are[
0,1
], [
-
1000,1436
] HU, and [
-
1000,837
] HU while the display windows are [
0.15,0.45
], [
400,1000
] HU, and [
-
160,400
] HU, respectively.


### 4.1. Comparison of Different Regularization Parameters

In the following experiments, all the parameters except the regularization parameterλ keep invariant for the same phantom with the same projection noise level. Three values of λ are selected, of which one is calculated by the proposed algorithm, another one is multiplied by 0.1, and the third one is multiplied by 10. The distance from the X-ray source to the center point of the phantom is twice the length of the image edge. The iteration of the algorithm is stopped when the relative error e
r
r
δ
=
δ
t
-
δ
t
-
1
/
δ
t is less than a stopping value (δ is calculated by (24)).To compare the difference between different selections of the regularization parameter, the human abdomen slice image is tested as an example. The projection data are simulated by 180 views of 2° step length over a 360° range, and 512 detector elements are distributed in fan-beam geometry covering the phantom. The noise levels added to the projection data are 0.0% and 0.1% Gaussian random noise. The results are displayed in Figures2 and 3. For the reason that biomedical images are often observed by a proper window to find more details, the images are displayed with a window [
-
160,400
] HU. The difference between the reconstructed image and the phantom image is displayed by a window [
-
90,90
] HU.Figure 2
The results of human abdomen slice simulation study. From top to bottom, the noise levels are 0.0% and 0.1% in turn. From left to right, the regularization parameters areλ
∗, 0.1
λ
∗, and 10
λ
∗. The display window is [
-
160,400
] HU.Figure 3
The difference between the reconstructed image and the original image (OI) of the human abdomen slice image. From top to bottom, the noise levels are 0.0% and 0.1% in turn. From left to right, the regularization parameters areλ
∗, 0.1
λ
∗, and 10
λ
∗. The display window is [
-
90,90
] HU.There are two criterions to evaluate the reconstructed image. One is the normalized mean absolute deviation (NMAD), defined as(28)
N
M
A
D
%
=
∑
i
,
j
μ
i
j
-
μ
i
j
truth
∑
i
,
j
μ
i
j
t
r
u
t
h
×
100
.The other one is the signal-to-noise ratio (SNR), defined as(29)
SNR
=
10
lg
⁡
∑
i
,
j
μ
i
j
truth
2
∑
i
,
j
μ
i
j
-
μ
i
j
truth
2
.The values of the two criterions are presented in Table1. Comparing the results with the same noise level of the λ situation and the 10
λ situation, the NMAD of the λ situation is smaller and the SNR of the λ situation is larger mostly, which proves that the image reconstructed by choosing the regularization parameter as λ is more close to the sample image. When it comes to the 0.1
λ situation, although the values of the two criterions are a little better, there are some artifacts appearing in the reconstructed image, leading to a decline of the image quality. In the middle column of Figures 2 and 3, some horizontal line artifacts appear in the images (regions D, E, and F in Figure 2 and the ellipse regions in Figure 3). It seems that there are more horizontal artifacts in the 0.0% noise level image. The probable reason is the noise added to the projection data since the noise covers the inconspicuous artifact in the 0.1% noise level image. So what is the reason for the fact that the NMAD and SNR of the 0.1
λ situation are better? One reasonable explanation might be the smoothing effect of the dictionary learning algorithm. When the iterative image is updated by (6) or (23), the sparse constraint added by dictionary learning method smoothes not only the noise but also the margin details. This effect becomes more significant when the regularization parameter is becoming larger. Therefore, the difference between the reconstructed image and the sample image in the margin area becomes more obvious, which is displayed in the left and right columns of Figure 3. By discovering this effect making the criterion worse, future work should be devoted to improving the reconstruction algorithm based on dictionary learning in order to smooth the noise and preserve the margin details meanwhile.Table 1
Quantitative evaluation of the results with different regularization parameters.
Noise level
δ
G (10
6
δ
λ
→
∞)
λ
NMAD (%)
SNR (dB)
0.0%
1.7094
0.9005
λ
1.5819
35.5246
0.1
λ
1.2406
37.6730
10
λ
2.0147
32.9517
0.1%
2.5269
5.7462
λ
2.0371
33.1220
0.1
λ
1.8742
34.5110
10
λ
2.1561
32.4527


### 4.2. Comparison of Different Reconstruction Algorithms

To present the advantage of AWR-ADSIR, the images reconstructed by AWR-ADSIR and some other algorithms (SART, GPBB, and ADSIR) are compared with the same projection data, initial conditions, and stopping criterions. The testing examples are the Shepp-Logan phantom and the human head slice image. The Shepp-Logan phantom is simulated by 120 views of 3° step length over a 360° range, and 512 detector elements are distributed in fan-beam geometry with three different noise levels while the head slice sample is simulated by 180 views of 2° step length over a 360° range. The reconstructed results are displayed in the four figures (Figure4 to Figure 7).Figure 4
Reconstructed images from low-dose projection data of the Shepp-Logan phantom. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
0.15,0.45
].With the comparative results calculated by (28) and (29) presented in Tables 2 and 3, the quality of all the reconstructed images is decreased with the noise level increasing. Among these four algorithms, SART generates the worst results. GPBB behaves well when the noise level is very low but when the noise level is beyond 0.1%, the quality of reconstructed image degenerates quickly. The regularization parameter in ADSIR is empirically selected according to [12]. The fact that the image quality and comparing criterions of ADSIR and AWR-ADSIR are mostly the same proves that the parameter selecting model in AWR-ADSIR is practical and efficient. In addition, the marginal details appearing in Figures 5 and 7 of algorithms ADSIR and AWR-ADSIR indicate the smoothing effect on the margins again. Since the high-contrast edge is smoothed by the algorithm, the structural boundaries appear in the difference image in Figure 7. Obviously, when the empirical regularization parameter is unknown, the ADSIR algorithm costs a large amount of time to determine the proper value by repeatedly operating the iterative process (usually more than ten times) while the AWR-ADSIR algorithm only operates the iterative process twice as Algorithm 1 shows. The model indeed reduces the time consumption to determine the value of the regularization parameter.Table 2
Comparing criterions of the results reconstructed by different algorithms (Shepp-Logan).
Algorithm
Criterion
Noise level
0.0%
0.1%
0.2%
SART
NMAD (%)
2.1510
2.9273
4.1032
SNR (dB)
33.5448
31.5220
28.7997
GPBB
NMAD (%)
1.3472
1.7826
3.0473
SNR (dB)
36.5443
34.6140
30.5004
ADSIR
NMAD (%)
0.8110
0.9553
1.1823
SNR (dB)
35.7740
34.8516
33.3413
AWR-ADSIR
NMAD (%)
0.8223
1.0639
1.1549
SNR (dB)
35.5354
34.5480
33.7801Table 3
Comparing criterions of the results reconstructed by different algorithms (head).
Algorithm
Criterion
Noise level
0.0%
0.1%
0.2%
SART
NMAD (%)
1.0904
2.2017
4.0009
SNR (dB)
36.8554
31.6575
26.4913
GPBB
NMAD (%)
0.4837
1.1631
2.7798
SNR (dB)
45.0405
36.7565
29.4104
ADSIR
NMAD (%)
0.5981
0.9370
1.1693
SNR (dB)
39.0515
34.9548
32.8972
AWR-ADSIR
NMAD (%)
0.6765
0.9438
1.1572
SNR (dB)
37.7225
34.9687
32.9703Figure 5
The difference between the reconstructed image and the original image (OI) of the Shepp-Logan phantom. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
-
0.05,0.05
].Figure 6
The results of human head slice simulation study. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
400,1000
] HU.Figure 7
The difference between the reconstructed image and the original image (OI) of the human head slice image. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
-
100
,
100
] HU.


### 4.3. Adaptive Weight Regularized GDSIR

The last experiment is applying the parameter selecting strategy into the GDSIR model. The dictionary displayed in Figure8 is learned from the overlapping patches of the original image in Figure 1. The reconstruction process of AWR-GDSIR is almost the same as AWR-ADSIR except that the dictionary has been constructed in advance and dose not change during the reconstruction process. The result shown in Figure 8 indicates that the proposed model is also suitable for the general dictionary condition.Figure 8
The image reconstructed by AWR-GDSIR and the global dictionary. The dictionary is displayed in window[
-
1,1
], which is constructed training based on the patches extracted from original image. The image is reconstructed by 0.0% noise level projection data displayed in window [
400,1000
] HU.


## 4.1. Comparison of Different Regularization Parameters

In the following experiments, all the parameters except the regularization parameterλ keep invariant for the same phantom with the same projection noise level. Three values of λ are selected, of which one is calculated by the proposed algorithm, another one is multiplied by 0.1, and the third one is multiplied by 10. The distance from the X-ray source to the center point of the phantom is twice the length of the image edge. The iteration of the algorithm is stopped when the relative error e
r
r
δ
=
δ
t
-
δ
t
-
1
/
δ
t is less than a stopping value (δ is calculated by (24)).To compare the difference between different selections of the regularization parameter, the human abdomen slice image is tested as an example. The projection data are simulated by 180 views of 2° step length over a 360° range, and 512 detector elements are distributed in fan-beam geometry covering the phantom. The noise levels added to the projection data are 0.0% and 0.1% Gaussian random noise. The results are displayed in Figures2 and 3. For the reason that biomedical images are often observed by a proper window to find more details, the images are displayed with a window [
-
160,400
] HU. The difference between the reconstructed image and the phantom image is displayed by a window [
-
90,90
] HU.Figure 2
The results of human abdomen slice simulation study. From top to bottom, the noise levels are 0.0% and 0.1% in turn. From left to right, the regularization parameters areλ
∗, 0.1
λ
∗, and 10
λ
∗. The display window is [
-
160,400
] HU.Figure 3
The difference between the reconstructed image and the original image (OI) of the human abdomen slice image. From top to bottom, the noise levels are 0.0% and 0.1% in turn. From left to right, the regularization parameters areλ
∗, 0.1
λ
∗, and 10
λ
∗. The display window is [
-
90,90
] HU.There are two criterions to evaluate the reconstructed image. One is the normalized mean absolute deviation (NMAD), defined as(28)
N
M
A
D
%
=
∑
i
,
j
μ
i
j
-
μ
i
j
truth
∑
i
,
j
μ
i
j
t
r
u
t
h
×
100
.The other one is the signal-to-noise ratio (SNR), defined as(29)
SNR
=
10
lg
⁡
∑
i
,
j
μ
i
j
truth
2
∑
i
,
j
μ
i
j
-
μ
i
j
truth
2
.The values of the two criterions are presented in Table1. Comparing the results with the same noise level of the λ situation and the 10
λ situation, the NMAD of the λ situation is smaller and the SNR of the λ situation is larger mostly, which proves that the image reconstructed by choosing the regularization parameter as λ is more close to the sample image. When it comes to the 0.1
λ situation, although the values of the two criterions are a little better, there are some artifacts appearing in the reconstructed image, leading to a decline of the image quality. In the middle column of Figures 2 and 3, some horizontal line artifacts appear in the images (regions D, E, and F in Figure 2 and the ellipse regions in Figure 3). It seems that there are more horizontal artifacts in the 0.0% noise level image. The probable reason is the noise added to the projection data since the noise covers the inconspicuous artifact in the 0.1% noise level image. So what is the reason for the fact that the NMAD and SNR of the 0.1
λ situation are better? One reasonable explanation might be the smoothing effect of the dictionary learning algorithm. When the iterative image is updated by (6) or (23), the sparse constraint added by dictionary learning method smoothes not only the noise but also the margin details. This effect becomes more significant when the regularization parameter is becoming larger. Therefore, the difference between the reconstructed image and the sample image in the margin area becomes more obvious, which is displayed in the left and right columns of Figure 3. By discovering this effect making the criterion worse, future work should be devoted to improving the reconstruction algorithm based on dictionary learning in order to smooth the noise and preserve the margin details meanwhile.Table 1
Quantitative evaluation of the results with different regularization parameters.
Noise level
δ
G (10
6
δ
λ
→
∞)
λ
NMAD (%)
SNR (dB)
0.0%
1.7094
0.9005
λ
1.5819
35.5246
0.1
λ
1.2406
37.6730
10
λ
2.0147
32.9517
0.1%
2.5269
5.7462
λ
2.0371
33.1220
0.1
λ
1.8742
34.5110
10
λ
2.1561
32.4527


## 4.2. Comparison of Different Reconstruction Algorithms

To present the advantage of AWR-ADSIR, the images reconstructed by AWR-ADSIR and some other algorithms (SART, GPBB, and ADSIR) are compared with the same projection data, initial conditions, and stopping criterions. The testing examples are the Shepp-Logan phantom and the human head slice image. The Shepp-Logan phantom is simulated by 120 views of 3° step length over a 360° range, and 512 detector elements are distributed in fan-beam geometry with three different noise levels while the head slice sample is simulated by 180 views of 2° step length over a 360° range. The reconstructed results are displayed in the four figures (Figure4 to Figure 7).Figure 4
Reconstructed images from low-dose projection data of the Shepp-Logan phantom. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
0.15,0.45
].With the comparative results calculated by (28) and (29) presented in Tables 2 and 3, the quality of all the reconstructed images is decreased with the noise level increasing. Among these four algorithms, SART generates the worst results. GPBB behaves well when the noise level is very low but when the noise level is beyond 0.1%, the quality of reconstructed image degenerates quickly. The regularization parameter in ADSIR is empirically selected according to [12]. The fact that the image quality and comparing criterions of ADSIR and AWR-ADSIR are mostly the same proves that the parameter selecting model in AWR-ADSIR is practical and efficient. In addition, the marginal details appearing in Figures 5 and 7 of algorithms ADSIR and AWR-ADSIR indicate the smoothing effect on the margins again. Since the high-contrast edge is smoothed by the algorithm, the structural boundaries appear in the difference image in Figure 7. Obviously, when the empirical regularization parameter is unknown, the ADSIR algorithm costs a large amount of time to determine the proper value by repeatedly operating the iterative process (usually more than ten times) while the AWR-ADSIR algorithm only operates the iterative process twice as Algorithm 1 shows. The model indeed reduces the time consumption to determine the value of the regularization parameter.Table 2
Comparing criterions of the results reconstructed by different algorithms (Shepp-Logan).
Algorithm
Criterion
Noise level
0.0%
0.1%
0.2%
SART
NMAD (%)
2.1510
2.9273
4.1032
SNR (dB)
33.5448
31.5220
28.7997
GPBB
NMAD (%)
1.3472
1.7826
3.0473
SNR (dB)
36.5443
34.6140
30.5004
ADSIR
NMAD (%)
0.8110
0.9553
1.1823
SNR (dB)
35.7740
34.8516
33.3413
AWR-ADSIR
NMAD (%)
0.8223
1.0639
1.1549
SNR (dB)
35.5354
34.5480
33.7801Table 3
Comparing criterions of the results reconstructed by different algorithms (head).
Algorithm
Criterion
Noise level
0.0%
0.1%
0.2%
SART
NMAD (%)
1.0904
2.2017
4.0009
SNR (dB)
36.8554
31.6575
26.4913
GPBB
NMAD (%)
0.4837
1.1631
2.7798
SNR (dB)
45.0405
36.7565
29.4104
ADSIR
NMAD (%)
0.5981
0.9370
1.1693
SNR (dB)
39.0515
34.9548
32.8972
AWR-ADSIR
NMAD (%)
0.6765
0.9438
1.1572
SNR (dB)
37.7225
34.9687
32.9703Figure 5
The difference between the reconstructed image and the original image (OI) of the Shepp-Logan phantom. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
-
0.05,0.05
].Figure 6
The results of human head slice simulation study. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
400,1000
] HU.Figure 7
The difference between the reconstructed image and the original image (OI) of the human head slice image. From top to bottom, the noise levels are 0.0%, 0.1%, and 0.2% in turn. From left to right, the reconstruction algorithms are SART, GPBB, ADSIR, and AWR-ADSIR, respectively. The display window is[
-
100
,
100
] HU.


## 4.3. Adaptive Weight Regularized GDSIR

The last experiment is applying the parameter selecting strategy into the GDSIR model. The dictionary displayed in Figure8 is learned from the overlapping patches of the original image in Figure 1. The reconstruction process of AWR-GDSIR is almost the same as AWR-ADSIR except that the dictionary has been constructed in advance and dose not change during the reconstruction process. The result shown in Figure 8 indicates that the proposed model is also suitable for the general dictionary condition.Figure 8
The image reconstructed by AWR-GDSIR and the global dictionary. The dictionary is displayed in window[
-
1,1
], which is constructed training based on the patches extracted from original image. The image is reconstructed by 0.0% noise level projection data displayed in window [
400,1000
] HU.


## 5. Conclusion

In most optimization problems, the determination of the regularization parameter is still a problem. In this paper, aiming to determine the regularization parameter of the algorithm based on dictionary learning, one model function, whose independent variableδ
λ
→
∞ can be calculated by the known projection data, is proposed depending on some modification on the objective function. When compared to some other algorithms, the images reconstructed by AWR-ADSIR and ADSIR are of similar quality, better than the one reconstructed by SART, and the proposed algorithm is much more robust to noise than GPBB. This indicates that the modification of the objective function does not degrade the performance of ADSIR. What is more, the parameter selection model is demonstrated to be rational by the fact that the image quality of the λ situation is better than the ones of 0.1
λ and 10
λ situations. However, when some other parameters (like the scale of the dictionary, the scale of the patch, and so on) change, the model function might result in some difference. However, it still works to look for the function relationship between these two parameters since the monotonous relation between δ
λ
→
∞ and λ remains.By validating the proposed selecting principle, the smoothing effect on image margins is discovered. Our future work will focus on improving the dictionary learning method, with the expectation to maintain the smoothing effect on noise regions and preserve marginal information better.



---

*Source: 831790-2015-10-04.xml*