Journal of Biosciences and Medicines, 2013, 1, 6-9 JBM
http://dx.doi.org/10.4236/jbm.2013.12002 Published Online October 2013 (http://www.scirp.org/journal/jbm/)
Towards n ove l regularization approaches to PET i m age
reconstruction
E. Karali, D. Koutsouris
Department of Electrical and Computer Engineering, National Technical University of Athens, Athens, Greece
Received August 2013
ABSTRACT
The purpose of this study is to introduce a novel em-
pirical iterative algorithm for medical image recon-
struction, under the short name MRP-ISWLS (Me-
dian Root Prior Image Space Weighted Least Squares).
Further, we assess the performance of the new algo-
rithm by comparing it to the simultaneous version of
known MRP algorithms. All algorithms are com-
pared in terms of cross-correlation and CNRs (Con-
trast-to-Noise Ratios). As it turns out, MRP-ISWLS
presents higher CNRs than the known algorithms for
objects of different size. Also MRP-ISWLS has better
noise manipulation.
Keywords: Image Reconstruction; Positron Emission
Tomography (PET); Small Animal Imaging; Median
Root Prior (MRP)
1. INTRODUCTION
Iterative techniques are divided into two main categories:
algebraic and statistical. Statistical techniques are classi-
fied to maximum-likelihood algorithms and least squares
methods. The most famous maximum likelihood tech-
nique is the expectation maximization maximum likeli-
hood (EM-ML) algorithm, which was first presented by
Shepp and Vardi. Image Space Reconstruction Algorithm
(ISRA) is a least square reconstruction method intro-
duced by Daube-Witherspoon and Muehllehner. It shows
better noise manipulation than EM-ML. Another least
squares algorithm is the Weighted Least Square tech-
nique (WLS), due to Anderson et al. WLS assumes that
random independent noise variables present different
standard deviations. The matrix of these standard devia-
tions consists of the expected projection data. WLS ac-
celerates the reconstruction process and results in recon-
structed images of better spatial resolution [1-3].
For the reduction of the noise many regularization
methods have been proposed, which reduce drastically
the noise with a small image resolution reduction. These
methods take into account a priori information for the
radioactivity spatial distribution inside the object under
examination [1]. The success of a regularization method
depends on the mathematical formula of the prior. Me-
dian root prior (MRP) [4] belongs to the most popular
priors. It is derived from a Gaussian distribution with
mean value the median value of reconstructed image
pixels in the vicinity of pixel i. The use of MRP results in
noise component reduction while at the same time it
preserves the edges.
The purpose of this study is on the one hand to intro-
duce a new empirical MRP image reconstruction algo-
rithm, under the short name MRP-ISWLS (Median Root-
Prior-Image Space Weighted Least Squares). MRP-
ISWLS will be an MRP version of ISWLS (Image Space
Weigthed Least Squares), which was introduced in [5].
ISWLS has ISRA properties in noise manipulation and
WLS acceleration of reconstruction process. To assess
the performance of the new iterative reconstruction me-
thod we have used phantom data produced from simu-
lating a prototype small-animal PET system. We com-
pared reconstruction data with MRP-EM-ML, MRP-
ISRA and MRP-WLS following the OSL (One Step Late)
philoshophy [6]. The methods presented here are applied
to 2D sinograms. Moreov e r, the simultaneous version of
the aforementioned algorithms is implemented.
We note that simultaneous versions of reconstruction
algorithms, that is, algorithms where all image pixels are
simultaneously updated in every iteration, are of great
interest because of their ability to be implemented in
parallel computing architectures, which decreases dras-
tically the total reconstruction time [1].
2. THEORY
In general, every iterative method relies on the hypothe-
sis that the projection data y is linearly connected to the
image x of radiopharmaceutical spatial distribution, ac-
cording to the equation:
xy
T
A =
(1)
where
A
is a matrix that characterizes the PET system

E. Karali, D. Koutsouris / Journal of Biosciences and Medicines 1 (2013) 6-9
7
being used for data acquisition. In bibliography this ma-
trix is called system or probability matrix and it projects
image data to sinogram domain (the term sinogram refers
to the projection data matrix) [1]. Every element αij of
the system matrix
A
represents the probability an an-
nihilation event emitted in image pixel
i
to be detected
in LORj (Line of Response). The significance of the
probability matrix lies on the valuable information re-
lated to the data acquisition process, that it can contain
(e.g. number of detector rings, number of detector ele-
ments in every ring, ring d iameter, diameter of transaxial
field of view, detector size, image size, spatial and angu-
lar sampling).
2.1. The ISWLS Algorithm
Consider an image discretized into
N
pixels and the
measured data
y
collected by M detector tubes. The
updating scheme of the ISWLS algorithm in the kth it-
eration is:
( )
2
1
1
2
1
11
M
ij j
j
kk
ii
MN k
iji ji
ji
ay
xx
a ax
=
−
−
′′
′
= =
=
∑
∑∑ (2)
2.2. The Median Root Prior (MRP)
It is derived from a Gaussian distribution with mean
value the median value of reconstructed image pixels in
the vicinity of pixel
i
.
Suppose that:
M
M
i
x
b
iexf 2
2
)(
)(
−
−
≈
(3)
where
),(ixmedMi =
, the median value of recon-
structed image pixels in the vicinity of pixel i. Then:
),(
),()(
)( ixmed
ixmedx
b
x
xf
xu
i
ii
i
i
i
−
=
∂
∂
=
(4)
The term b
∈
[0,2] determines the degree of smoo-
thing in reconstruction images. If b = 0 no prior is ap-
plied. Big values of b cau se over-smoothing, while small
values of b result in images with high resolution but with
increased noise.
2.3. MRP Algorithms
According to the one-step-late philosophy, where the
prior is applied to the previous radiopharmaceutical dis-
tribution estimation, we can extract an empirical iterative
formula for the ISWLS algorithm in combination with
MRP prior. The new empirical iterative algorithm has
updating scheme:
2
11
11 2
1
''
11 '1
( ,)
1( ,)
i
i
M
ij j
kj
ki
ikk
MN
ik
iji ji
kji
ay
x
x
xmed xi
ba ax
med xi
−=
−−
−
−= =
=−
+

∑
∑∑
(5)
The MRP-EM-ML updating scheme in kth iteration is:
1
11
11
''
1'1
( ,)
1( ,)
i
i
kMij j
ki
ikk
N
jk
iij i
ki
ay
x
xxmed xiax
bmed xi
−
−−
=−
−=
=−

+


∑∑
(6)
The MRP-ISRA updating scheme in kth iteration is:
11
11 1
111
( ,)
1( ,)
ii
i
M
ij j
kj
ki
ikkMN k
iiji j
kji
ay
x
xxmed xia ax
bmed xi
′
−=
−− −
′
−′
= =
=−
+
∑
∑∑
(7)
The MRP-WLS updating s cheme in kth iteration is
2
1
11 2
11
''
1'1
( ,)
1( ,)
j
i
i
kMij
ki
ikk
N
j
ik
ij i
ki
ay
x
xxmed xi
bax
med xi
−
−−
=−
−=
=−
+

∑
∑
(8)
3. RESULT S
For the evaluation of the iterative reconstruction methods
presented in Section 2, projection data of a Derenzo-type
phantom have been used. The Derenzo-type phantom
consists in sets of rods filled with F18, with diameters 4.8,
4, 3.2, 2.4, 1.6, and 1.2 mm, and the same separation
between surfaces in the corresponding sets. The rods
were surrounded by plastic (polyethylene). Data were
produced using Monte Carlo simulation of a small-ani-
mal PET scanner, consisting of two detector heads.
Further, 18 × 106 coincidence events were collected.
Projection data was binned to a 2D s inogram, 55 pixels ×
170 pixels in size, which means that data from 55 TORs
(Tube of Response) per rotation angle were collected and
170 totally angular samples were used. Since the two
detector heads rotate from 0˚ to 180˚ the angular step
size was 1.0647˚.
The system matrix was derived from an analytical
method and calculated once before reconstruction. Each
element
ij
a
was computed as the area of intersection Eij,
of TORj (Tube-of-Response) with image pixel
i
. The
calculated system matrix is a sparse matrix. It consists of
zero-valued elements in majority that have no contribu-
tion during iterative reconstruction process. So, only the
non-zero elements were stored, resulting in significant
reduction in system matrix size and consequently in re-
quired storage. The reconstructed 2D images were 128
pixels × 128 pixels in size, thus the system matrix con-
sisted of 55 × 170 × 128 × 128 elements with 4.33%

E. Karali, D. Koutsouris / Journal of Bios ci ences and Medicine s 1 (2013) 6-9
8
sparsity.
The initial image estimate for all MRP algorithms was:
1
, i1,2,,N
i
M
j
j
o
y
xN
=
= =
∑

(9)
where yj is the value of the
th
j
sinogram element and N
represents the total number of image pixels (
N
= 128 ×
128 in this implementation). The value of b was 0.00 1.
Figu re 1 shows the reconstructed transaxial images
with MRP-EMML, MRP-ISRA, MRP-WLS and MRP-
ISWLS after 1, 10 and 50 iterations.
In Figure 2, cross-correlation coefficient
c
of every
iterative method is plotted versus the number of itera-
tions. The cross-correlation coefficient
c
was calcu-
lated according to the equation:
()( )
() ()
11
22
11 11
NN
ij ij
ij
NN NN
ij ij
ij ij
IreconIrecon IrealIreal
c
Irecon IreconIrealIreal
= =
= == =
−−
=
−−
∑∑
∑∑ ∑∑
(10)
where
reconI
and
realI
are the reconstructed image
and the true phantom activity image mean values, re-
spectively. Cross-correlation coefficient is a similarity
Figure 1. Reconstructed images with: (a) MRP-EMML, (b)
MRP-ISRA; (c) MRP-WLS, and d) MRP-ISWLS, after 1, 10
and 50 iterations respectively.
Figure 2. Cross-correlation coefficient versus the number of
iterations for MRP-EMML, MRP-ISRA, MRP-WLS, and
MRP-ISWLS.
measure between reconstructed and real radio distribu-
tion image. Its values are in the range [−1,1]. Value
1= c
corresponds to fully correlated images.
Except for the cross-correlation coefficient that shows
the average performance of the reconstruction methods,
local contrast-to-noise ratios (CNR) [7] for rods with
different diameters were calculated. CNRs for 4.8, 3.2,
and 1.6 mm rods diameter were computed, using squared
regions-of-interest (ROIs), 4.55, 3.85 and 2.15 mm in
size, respectively. The ROIs were placed inside the cor-
responding objects. The number of selected ROIs was
equal to the number of same sized objects. ROIs of the
same sizes were positioned in three different background
areas, each. CNRROI was defined as:
ROI
ROIROI
Backg
Backgobj
ROI
RR
CNR σ
−
=
, (11)
where
ROI
obj
R
is the mean value of reconstructed objects
in the corresponding ROIs, and
ROI
Backg
R
is the mean
value of the three background ROIs in each case. Further,
ROI
Backg
σ
is the standard deviation of background values
in the corresponding ROIs. The graphs in Figure 3 illus-
trate the variation of CNRROI with respect to the number
of iterations for the three different object diameters.
4. DISCUSSIO N
MRP-ISWLS presents higher cross-correlation values
than MRP-EM-ML and MRP-ISRA. Although it shows
the same high values of cross-correlation coefficient as
MRP-WLS during the first 50 iterations it has better
noise manipulation. In our research cross-correlation
coefficient reaches up to 0.75, and not to 1. This is be-
cause extra corrections prior or during the reconstruction
should be made, like attenuation, scatter and random
b
a
c
d
1iter 10ite r50iter

E. Karali, D. Koutsouris / Journal of Biosciences and Medicines 1 (2013) 6-9
9
Figure 3. CNRs versus iterations for: (a) 4.8 mm; (b) 3.2 mm
and (c) 1.6 mm object diameter.
corrections. Despite the fact that these corrections have
not been made the final result of the average perform-
ance of MRP-ISWLS is not altered.
As illustrated in Figu re 3 MRP-ISWLS presents high
CNR ratios from the first iterations, higher than MRP-
EM-ML and MRP-ISRA. Although it shows similar per-
formance to MRP-WLS, its CNR ratios do not degrade
after 50 iterations but tend to be stabilized. So, MRP-
ISWLS presents a better noise manipulation than MRP-
WLS.
In order to determine a satisfactory value of b, various
values were applied to MRP-ISWLS. These values were
0.001, 0.01, 0.1, 1, 1.5. The value b = 0.001 presented
higher CNR values in comparison to ISWLS’s results.
MRP-ISRA and MRP_ISWLS are slower than MRP-
EM-ML and MRP-WLS during the first 9 iterations (79
s/iteration). Their reconstruction speed is gradually im-
proving with increasing number of iterations. The reason
for slow reconstruction process during the first iterations
lies in the time needed for backprojection computations
in the first iteration.
5. CONCLUSION
In this work, different simultaneous iterative reconstruc-
tion schemes were applied to data acquired from a simu-
lation of a small-animal PET scanner. A new iterative
scheme was introduced, namely MRP-ISWLS.
6. ACKNOWLEDGEMENTS
This research has been co-financed by the European Union (European
Social Fund-ESF) and Greek national funds through the Operational
Program “Education and Lifelong Learning” of the National Strategic
Reference Framework (NSRF)-Research Funding Program: THALES:
Reinforcement of the interdisciplinary and/or inter-institutional re-
search and innovation.
REFERENCES
[1] Phelps M.E. (2004) PET Molecular Imaging and its Ap-
plications. Chapter 1, Springer-Verlag, New York.
[2] Daube-Witherspoon, M.E. and Muehllehner, G. (1986)
An iterative image space reconstruction algorithm suita-
ble for volume ECT. IEEE Transactions on Medical Im-
aging, 5, 61-66.
[3] Anderson, M.M., Mair, B.A., Rao, M. and Wu, C.H.
(1997) Weighted least-squares reconstruction methods for
positron emission tomography. IEEE Transactions on
Medical Imaging, 16, 159-165.
http://dx.doi.org/10.1109/42.563661
[4] Alenius, S. Ruotsalainen, U. and Astola J. (1997) Using
local median as the location of the prior distribution in it-
erative emission tomography image reconstruction. Pro-
ceedings of IEEE Medical Image Conference, 45, 3097-
3104.
[5] Karali, E., Pavlopoulos, S., Lambropoulou, S. and Kout-
souris, D. (2011). ISWLS: A novel algorithm for image
reconstruction in PET. IEEE Transactions on Information
Technology in Biomedicine, 15, 381-386.
http://dx.doi.org/10.1109/TITB.2010.2104161
[6] Green, P. J. (1990) On use of the EM algorithm for penal-
ized likelihood estimation. Journal of Royal Statistical
Society Series B , 52, 443-452.
[7] Cherry, S.R., Sorenson, J.A. and Phelps M.E. (2003)
Physics in Nuclear Medicine. Chapter 15, SAUNDE RS-
Elsevier, Philadelphia.