Towards novel regularization approaches to PET image reconstruction

doi:10.4236/jbm.2013.12002

Paper Menu >>

Journal Menu >>

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/)

OPEN ACCESS

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

Email: [email protected]

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:

(1)

where

is a matrix that characterizes the PET system

E. Karali, D. Koutsouris / Journal of Biosciences and Medicines 1 (2013) 6-9

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

represents the probability an an-

nihilation event emitted in image pixel

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

pixels and the

measured data

collected by M detector tubes. The

updating scheme of the ISWLS algorithm in the kth it-

eration is:

( )

ij j

MN k

iji ji

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

Suppose that:

iexf 2

)(

−

≈

(3)

where

),(ixmedMi

, the median value of recon-

structed image pixels in the vicinity of pixel i. Then:

),(

),()(

)( ixmed

ixmedx

−

∂

(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:

11 2

11 '1

( ,)

1( ,)

ij j

ikk

iji ji

kji

xmed xi

ba ax

med xi

−=

−−

−

−= =

=−

+



∑

∑∑

(5)

The MRP-EM-ML updating scheme in kth iteration is:

1'1

( ,)

1( ,)

kMij j

ikk

iij i

xxmed xiax

bmed xi

−

−−

=−

−=

=−







∑∑

(6)

The MRP-ISRA updating scheme in kth iteration is:

11 1

111

( ,)

1( ,)

ij j

ikkMN k

iiji j

kji

xxmed xia ax

bmed xi

′

−=

−− −

′

−′

= =

=−

∑

∑∑

(7)

The MRP-WLS updating s cheme in kth iteration is

11 2

1'1

( ,)

1( ,)

kMij

ikk

ij i

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

was computed as the area of intersection Eij,

of TORj (Tube-of-Response) with image pixel

. 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

sparsity.

The initial image estimate for all MRP algorithms was:

, i1,2,,N

= =

∑



(9)

where yj is the value of the

sinogram element and N

represents the total number of image pixels (

= 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

of every

iterative method is plotted versus the number of itera-

tions. The cross-correlation coefficient

was calcu-

lated according to the equation:

()( )

() ()

11 11

ij ij

NN NN

ij ij

IreconIrecon IrealIreal

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

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

CNR σ

−

, (11)

where

ROI

obj

is the mean value of reconstructed objects

in the corresponding ROIs, and

ROI

Backg

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

1iter 10ite r50iter

E. Karali, D. Koutsouris / Journal of Biosciences and Medicines 1 (2013) 6-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.