ATTENUATION CORRECTION METHODS FOR

POSITRON EMISSION TOMOGRAPHY

By

CHEN-HSIEN WU

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT

OF THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1998

Dedicate to my parents,

Mr. Min-Ming Wu and Ms. Jau-Jr Wu,

and my sister,

Hwa-Yu.

ACKNOWLEDGEMENTS

I would like to thank my advisor, mentor and best friend, Professor John M.

M. Anderson for his endless patience and support. His open-mindedness has made

the time I have spent working for him become one of the most precious experiences

in my life.

Thanks go to Professor Bernard Mair for providing me lots of information and

suggestions which have strong influence on my research.

Thanks go to Professor Murali Rao, who has helped me find the solutions to

numerous problems that I have encountered while working on this dissertation.

Thanks go to Professor Jose Principe and Professor William Edmonson for

offering precious opinions which make this dissertation better.

I would like to thank Dr. Chuan Wang, who advised me on heading in the right

direction. Thanks to Pei-Fang, Hsin-Ho, Byron, Kali, Lisa, Chris, Raymond, Ray,

and Anke for their priceless support in helping me to walk through this journey

and overcome all difficulties.

Best thanks go to my family and ancestors.

TABLE OF CONTENTS

ACKNOWLEDGEMENTS .......................... iii

LIST OF FIGURES .............................. vi

LIST OF TABLES ................................ x

ABSTRACT .................................. xi

1 INTRODUCTION ............................. 1

1.1 Overview: Positron Emission Tomography ............... 2

1.2 Emission Data .............. ...................... 5

1.3 Clinical and Research Applications .................... 7

1.4 Advantages and Disadvantages ........................ 8

2 BACKGROUND .............................. 10

2.1 Basic Physics of PET ......................... 10

2.2 Poisson M odel ............................. 14

2.3 Errors in Positron Emission Tomography ................ 15

2.3.1 Background on Attenuation .... .. 15

2.3.2 Accidental Coincidences . .... .. 20

2.3.3 Scatter . .. .. 22

2.3.4 Other Errors . ... .. 24

2.4 Transmission Tomography in PET ... .. 24

2.5 Modified Poisson Model . ... .. 26

2.6 Existing Attenuation Correction Methods .... .. 29

2.6.1 Standard Correction Method .... .. 29

2.6.2 Reprojection Methods ... ................... .. 30

2.6.3 Simultaneous Estimation Methods .... .. 31

3 MAXIMUM LIKELIHOOD METHOD .... .. 33

3.1 Maximum Likelihood Estimator . ... .. 33

3.2 Concavity of The Log Likelihood Function .... .. 35

3.2.1 Example to Illustrate Concavity of the Log Likelihood Func-

tion . .. 37

3.3 Optimization Procedure . ... .. 38

3.3.1 x-step . .. .. 39

3.3.2 si-step . .. .. 40

3.3.3 Monotonically Increasing Value of the Log Likelihood Func-

tion in x-step . ... .. 41

3.3.4 Property of the ML Estimator .... .. 44

4 LEAST SQUARES REPROJECTION METHOD ... ..46

4.1 LS Algorithm for Estimating the Attenuation Density ...... .. 46

4.2 LS Reprojection Method . ... .. 47

5 SIM ULATION . ... .. 49

5.1 Synthetic Data . ... .. 50

5.1.1 ML Method Implementation of Algorithms ... .. 51

5.1.2 EMML Algorithm with Corrected Data ... .. 51

5.1.3 LS Reprojection Method . ... .. 51

5.1.4 Non-smoothed Typical Transmission Scan Data 52

5.1.5 Smoothed Typical Transmission Scan Data ... ..63

5.1.6 Non-smoothed Short Transmission Scan Data ........ 71

5.1.7 Smoothed Short Transmission Scan Data ... .. 77

5.2 Real Data . ... .. 85

5.2.1 Overview of Implementation ... .. 86

5.2.2 Results . .. .. 87

6 CONCLUSIONS . ... .. 92

6.1 Future W ork . ... .. 93

A EMML Algorithm for PET . ... .. 95

REFERENCES . ... .. 97

BIOGRAPHICAL SKETCH . ... .. 101

LIST OF FIGURES

1.1 Quest 250 from GE Medical Systems (courtesy of GE Medical Sys-

tem s) . . 3

1.2 A positron travels a short distance and annihilates with a nearby

electron to produce a pair of photons that fly off at 180-degree

angle. These photons are detected by a ring of detectors 6

2.1 Annihilation Process . ... .. 11

2.2 Scanditronix PC4096-15WB detector module: 16 crystals, 1 light

guide, and 2 dual cathode PMTs .... .. 13

2.3 Attenuation due to absorption . ... .. 16

2.4 Attenuation due to scattering . ... .. 18

2.5 Basic concept of attenuation . ... .. 19

2.6 Accidental coincidence in PET . ... .. 21

2.7 Scatter in PET . ... .. 23

2.8 Blank scan . ... .. 25

2.9 Transmission scan . ... .. 27

5.1 Reconstruction of cardiac emission phantom using non-smoothed

typical transmission scan data (a) phantom (b) 1st global iteration

(c) 2nd global iteration (d) 3rd global iteration ... .. 53

5.2 Reconstruction of cardiac attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) 1st

global iteration (c) 2nd global iteration (d) 3rd global iteration 54

5.3 Reconstruction of cardiac emission phantom using non-smoothed

typical transmission scan data (a) emission phantom (b) EMML

algorithm with corrected data (c) LS reprojection method (d) ML

m ethod . .. .. 56

5.4 Line plot of the 50th row of the reconstruction of cardiac emission

phantom using non-smoothed typical transmission scan data 57

5.5 Reconstruction of cardiac attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) ML

method (c) LS reprojection method .... .. 59

5.6 ML method (a) emission phantom (b) emission image (c) attenua-

tion phantom (d) attenuation image .... .. 60

5.7 Reconstruction of chest emission phantom using non-smoothed typ-

ical transmission scan data (a) emission phantom (b) EMML algo-

rithm with corrected data (c) LS reprojection method (d) ML method 61

5.8 Reconstruction of chest attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) ML

method (c) LS reprojection method .... .. 62

5.9 Reconstruction of cardiac emission phantom using smoothed typical

transmission scan data (a) emission phantom (b) EMML algorithm

with corrected data (c) LS reprojection method (d) ML method 64

5.10 Line plot of the 50h row of reconstruction of cardiac emission phan-

tom using the smoothed transmission scan data ... ..65

5.11 Reconstruction of cardiac attenuation phantom using the smoothed

typical transmission scan data (a) attenuation phantom (b) ML

method (c) LS reprojection method ... .. 67

5.12 ML method (a) emission phantom (b) emission image (c) attenua-

tion phantom (d) attenuation image .... .. 68

5.13 Emission image of the chest phantom using smoothed typical trans-

mission scan data (a) emission phantom (b) EMML algorithm with

corrected data (c) LS reprojection method (d) ML method 69

5.14 Reconstruction of chest attenuation phantom using smoothed typi-

cal transmission scan data (a) attenuation phantom (b) ML method

(c) LS reprojection method . ... .. 70

5.15 Reconstruction of cardiac emission phantom using the non-smoothed

short transmission scan data (a) emission phantom (b) EMML al-

gorithm with corrected data (c) LS reprojection method (d) ML

m ethod . .. ..73

5.16 Reconstruction of cardiac attenuation phantom using non-smoothed

short transmission scan data (a) attenuation phantom (b) ML method

(c) LS reprojection method . ... .. 74

5.17 ML method (a) emission phantom (b) emission image (c) attenua-

tion phantom (d) attenuation image .... .. 75

5.18 Reconstruction of chest emission phantom using non-smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm

with corrected data (c) LS reprojection method (d) ML method 76

5.19 Reconstruction of chest phantom using non-smoothed short trans-

mission scan data (a) attenuation phantom (b) ML method (c) LS

reprojection method . ... .. 78

5.20 Reconstruction of cardiac emission phantom using smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm

with corrected data (c) LS reprojection method (d) ML method 79

5.21 Reconstruction of cardiac attenuation phantom using smoothed

short transmission scan data (a) attenuation phantom (b) ML method

(c) LS reprojection method . ... .. 81

5.22 ML method (a) emission phantom (b) emission image (c) attenua-

tion phantom (d) attenuation image .... .. 82

5.23 Reconstruction of chest emission phantom using smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm

with corrected data (c) LS reprojection method (d) ML method 83

5.24 Reconstruction of chest attenuation phantom using smoothed short

transmission scan data (a) attenuation phantom (b) ML method (c)

LS reprojection method . ... .. 84

5.25 Emission images using high quality smoothed transmission scan

data (a) Filtered back-projection method (b) EMML algorithm

with corrected data (c) LS reprojection method (d) ML method 88

5.26 Attenuation images using high quality smoothed transmission scan

data (a) ML method (b) LS reprojection method ... 90

5.27 Emission image using high quality non-smoothed transmission scan

data (a) EMML algorithm with corrected data (b) ML method. 91

LIST OF TABLES

5.1 PSNR (dB) of emission images using typical transmission scan data 71

5.2 PSNR (dB) of reconstruction of cardiac emission phantom using

short transmission scan data . ... .. 85

Abstract of Thesis Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

ATTENUATION CORRECTION METHODS FOR

POSITRON EMISSION TOMOGRAPHY

By

Chen-Hsien Wu

May 1998

Chairman: Dr. John M M Anderson

Major Department: Electrical and Computer Engineering

Positron emission tomography (PET) is an imaging modality that provides bio-

chemical information about the human body. Among the sources of error in PET,

attenuation is the one that impacts image quality the most. The standard method

for correcting attenuation consists of the following two step procedure. First, a

blank scan and transmission scan are used to estimate the survival probabilities.

Then, the data for each pair of detectors is divided by the corresponding survival

probability estimate. This corrected data and a reconstruction algorithm, such as

the filtered back-projection algorithm, are then used to generate an estimate of

the emission intensity, which is displayed as an image.

In this dissertation, two iterative methods are proposed for correcting atten-

uation errors in PET. The first method is a maximum likelihood (ML) method

that simultaneously estimates the emission intensity and attenuation density. The

second one is a least squares (LS) reprojection method, where the survival prob-

abilities are estimated by reprojecting a LS estimate of the attenuation density.

The emission intensity is estimated using the expectation maximization maximum

likelihood (EMML) algorithm with a probability matrix that has been modified by

the survival probability estimates. Both methods are guaranteed to produce non-

negative estimates of the emission intensity and the attenuation density, provided

that the initial estimates are nonnegative.

From simulation studies, it appears that the major advantage of the proposed

methods over the standard method with the filtered back-projection algorithm

is that they (1) have less noise, (2) do not have artifacts outside the region of

interest, and (3) they provide high quality images even for short transmission

scans (less than 10 minutes).

CHAPTER 1

INTRODUCTION

For both clinical and research applications, medical imaging modalities are

used to obtain important information on the functionality and anatomy of the

human body. Imaging modalities such as X-ray computed tomography and mag-

netic resonance imaging provide detailed images of the anatomical structure of

the body. However, in many applications, it is more beneficial to have a modality

that determines functional information such as blood flow rates, oxygen utiliza-

tion rates, and glucose utilization rates within the region-of-interest (ROI) (e.g.,

brain, heart) [28]. For example, glucose utilization rates are used in oncology to

diagnose tumors and to determine the efficacy of a treatment.

Functional information on the human body can be obtained using (1) positron

emission tomography (PET), (2) single-photon emission computed tomography

(SPECT) [16], and (3) functional magnetic resonance imaging [51]. Among them,

PET is probably the most flexible imaging modality. For example, PET can

image which parts of the brain are functioning while a subject is performing a

speaking tasking. This is not possible in functional magnetic resonance imaging

because the scanners are too loud [55]. An advantage of PET over SPECT is

that the radionuclides in PET can be incorporated into molecules (or analogs) of

biological interest [56].

In this dissertation, we propose two methods for correcting a certain error in

PET called attenuation. The outline of the dissertation is as follows. An overview

of PET is provided in the remainder of Chapter 1. In Chapter 2, we discuss the

physics behind PET, introduce a mathematical model for the data, and review

existing attenuation correction methods. We introduce the proposed maximum

likelihood and least-squares reprojection methods in Chapter 3 and Chapter 4,

respectively. In Chapter 5, we investigate the performance of these methods using

synthetic and real data. Finally, we discuss our conclusions and future research

direction in Chapter 6.

1.1 Overview: Positron Emission Tomography

PET is exceptional in its capability to quantitatively determine functional

information within a ROI. Since George Hevesy won a Nobel Prize for his pio-

neering research on radionuclides in 1943, numerous researchers have investigated

methods to obtain functional information on the human body. However, the first

PET machine was not built until 1975 by Michael Phelps, Edward Hoffman, and

their colleagues [42]. Since then, scientists and engineers have made significant

advances in the hardware, radionuclides, and reconstruction algorithms used for

PET. In Fig. 1.1 is a PET scanner of the Quest series from the GE Medical

Systems.

We first introduce the major components of a PET scan, which are as follows:

1. A compound is labeled with a radionuclide.

Figure 1.1. Quest 250 from GE Medical Systems (courtesy of GE Medical Sys-

tems)

2. The labeled compound or radiopharmaceutical is introduced into the sub-

ject's body.

3. The data are collected.

4. An image is produced that represents the concentration of the radiophar-

maceutical within the ROI.

Radionuclides in PET are (1) incorporated into compounds that are naturally

used by the body, (2) emit positrons, and (3) have short half-lives. A radionuclide's

half-life must be sufficiently long so that there is enough data to get quality

images [38]. The most common radionuclides in PET include "C, 'O50, and

3N, which are used to label glucose, water and amino acid, respectively. The

most widely used radiopharmaceutical is SF-labeled deoxyglucose (FDG), which

is used to study glucose metabolism. Because of the short half-lives of positron-

emitting radionuclides (2 minutes, 10 minutes, 20 minutes, and 110 minutes for

"11C, O150, 13N, and SF, respectively), they need to be produced on-site using a

cyclotron. The cyclotron creates radionuclides by smashing the atomic shell of

certain molecules with protons. Cyclotrons are expensive and costly to maintain

because they require at least two operators.

SPECT uses radionuclides, such as 99mTc, 201TI, and 1231, which have half-

lives that range from 10 hours to several days. Because of their long half-lives,

radionuclides in SPECT can be made in advance and shipped to hospitals so an

on-site cyclotron is unnecessary. However, long half-lives mean that patients are

exposed to radioactivity for long periods.

5

After being administered to a subject, the labeled compound will localize in the

ROI. For example, FDG is absorbed preferentially by the brain because glucose

is a source of energy [11]. Another example is the use of 150 based radiopharma-

ceuticals to index injured tissues by monitoring the oxygen metabolism of tissues

[52].

1.2 Emission Data

As mentioned previously, the radiopharmaceutical localizes in the ROI. When

the radionuclide decays, it emits positrons which travel a short distance (3-4 mm)

before annihilating with nearby electrons. The areas with a high concentration of

the radiopharmaceutical emit more photons than the areas of low concentration.

When a positron annihilates with an electron, two high-energy (511 keV) photons

are produced that fly off at an almost 180-degree angle from each other. This

annihilation process is shown in Fig. 1.2. The photons go through the subject and

are registered by rings of detectors surrounding the subject. When two detectors

detect a pair of photons within a short period of time, these photons are said to

be in coincidence (see Fig. 1.2). For each pair of detectors, there is a count that

corresponds to the number of observed coincidences. The emission data is the

collection of these counts.

From the emission data, we can estimate the concentration of the radiophar-

maceutical using a reconstruction algorithm. A computer is then used to display

images where the intensity is proportional to the estimated concentration of the

radiopharmaceutical. In the next section, some clinical and research applications

of PET will be presented.

Figure 1.2. A positron travels a short distance and annihilates with a nearby

electron to produce a pair of photons that fly off at 180-degree angle. These

photons are detected by a ring of detectors

1.3 Clinical and Research Applications

PET is widely used in pathology and oncology because it can quantitatively

measure blood flow, and utilization rates of glucose, C02, and fatty acid ammo-

nia. PET provides functional information about tissues or organs, which cannot

be obtained by modalities such as magnetic resonance imaging and computed to-

mography [52]. It is also a powerful tool for investigating brain functions such

as emotion, thinking, vision, and speech [5, 30]. Since the major energy source

of the brain is glucose, FDG is used to map the brain's energy consumption. In

other words, if a specific region of the brain has increased functionality, then the

regional blood flow increases to provide the necessary energy. As a result, the con-

centration of the radionuclide used in FDG increases in that region. Using PET

to map the concentration of the radionuclide, we can observe the higher function-

ality in that region of the brain. Based on the concepts discussed above, scientists

have used PET images to assess the brains of normal children and children under

abuse. They found that the children who were abused had slower developments

in brain functions such as speech and vision as compared to normal children [7].

PET has clinical applications as well. For example, PET has the ability to

measure the size of infarcts in the heart after a coronary attack, to diagnose dis-

eases such as Alzheimer's [45] and epileptic seizures [36], and to survey the effects

of treatments such as radiation therapy, chemotherapy, and hormone therapy. In

addition, whole-body PET is used to spot metastases, which is very difficult for

other modalities because of the small sizes and the random occurrences of metas-

tases inside cancer patients' bodies. The effect of the chemotherapy on cancer

patients can be evaluated by visualizing the sizes and locations of metastases on

PET images [1].

PET is also used in psychological studies. For instance, it was used in the

early 1990s several times by courts in California to determine whether neural

malfunctions in the brain caused a suspect's violent behavior. Another interesting

experiment was the monitoring of six bank officials' brain functions when they

were exposed to a video of an armed bank robbery. This was done to find out the

connection between stress and visual re-experience [22].

1.4 Advantages and Disadvantages

PET is a better imaging modality than SPECT in the sense of offering better

efficiency, sensitivity, and accuracy as compared to SPECT [11]. Also, unlike

SPECT, the radionuclides used in PET are easily incorporated into compounds,

such as glucose, that are naturally used by the body. Using different radionuclides,

PET is able to image different biochemical activities while functional magnetic

resonance imaging only can image hydrogen activity. Additionally, PET can track

chemical changes that do not affect blood flow enough for functional magnetic

resonance imaging to measure.

Unfortunately, PET has some major disadvantages such as its high cost, high

complexity, and the requirement of highly skilled personnel. The cost of a PET

scanner is similar to a magnetic resonance imaging scanner. However, the cy-

clotron used to generate the radionuclides is also very expensive. Additionally,

at least three skilled professionals one to maintain the scanner, one to run the

9

cyclotron and another one to produce radionuclides are needed to operate a

PET laboratory. For SPECT, it is not necessary to create the radionuclides on

site because the radionuclides used in SPECT can be stored. Even though its

image quality is not as good as that of PET, SPECT is more widely available in

practice because it is significantly cheaper.

CHAPTER 2

BACKGROUND

In this chapter, we review the basic physics of PET, introduce the Poisson

model, and discuss sources of errors in the data. Since the focus of this dissertation

is attenuation correction, we also review existing attenuation correction methods.

2.1 Basic Physics of PET

As discussed in Chapter 1, the subject is given a radiopharmaceutical that

emits positrons. When a positron annihilates with a nearby electron, two photons

are produced that depart at an angle near 180 degrees. The region of interest

(ROI) is surrounded by rings of detectors that detect the photons. We will use

the term tube to refer to a pair of detectors. When a pair of photons is incident

on a tube within a small time interval, a counter is increased (see Fig. 2.1). The

emission data are the counts associated with the various tubes.

The detectors are shielded from radiation coming from outside the region of

interest by heavy lead rings at the front and back of the gantry. In order to

switch between the slice-collimated mode and the three-dimensional mode, most

commercial PET scanners contain spherical thin interplane rings called septa.

When the scanner is in slice-collimated mode, the septa only allows in-plane tubes

to detect photons. When a scanner operates in full three-dimensional mode, the

Coincidence

Timing -

Circuit

Figure 2.1. Annihilation Process

septa are withdrawn to permit photons to be registered by in-plane tubes and

cross-plane tubes [37].

A detector module consists of an array of crystals coupled to several photomul-

tiplier tubes. In Fig. 2.2, a detector module of the Scanditronix PC4096-15WB

whole-body PET scanner is shown. A detector ring contains several hundred de-

tector modules (e.g., 384 for the Siemens-CTI ECAT EXACT 921 scanner). When

a crystal is hit by a 511 keV photon, it interacts with the photon and generates nu-

merous light photons. Then, the light guide connected to the crystal transfers the

light photons to the nearest photomultiplier tube (PMT). The PMTs collect light

photons and convert them into electrical signals that are stored in a computer.

The outputs of the PMTs are decoded to determine which detector module was

hit by a photon. The most widely used crystal material is Bismuth-Germanate

because it has good properties when it interacts with 511 keV photons such as a

high light photon producing rate, short decay time, and a fast energy rise time

[26].

Once the data is collected, algorithms are used to reconstruct the PET im-

ages. The most popular method is the filtered back-projection algorithm [45]. It

is extremely fast and easy to implement, but the images are noisy, have nega-

tive values, and contain artifacts outside the ROI. Alternative approaches include

statistically based methods such as maximum likelihood [28, 46], maximum a

posterior [20, 21, 30, 31], least squares [15, 40, 24, 46, 49, 50], and weighted least-

squares [3, 12, 17] methods. An advantage of statistically based methods is that

they account for the errors in the data and do not have artifacts outside the ROI.

A disadvantage of them is that they are computationally demanding.

PMT

\ Light Guide

Crystals

Figure 2.2. Scanditronix PC4096-15WB detector module: 16 crystals, 1 light

guide, and 2 dual cathode PMTs

2.2 Poisson Model

Since the data collection in PET is actually a photon counting process, it is

reasonable to assume that the emission process is a Poisson process with unknown

mean [38]. Based on this assumption, Shepp and Vardi proposed a model called

the Poisson model [46] in 1982. Their model assumes that the data is error free

and obeys a Poisson process.

In the Poisson model, the ROI is discretized into a grid of J voxels. Let Nj

denote the number of emissions in voxel j, where Nj is assumed to be Poisson with

unknown mean xj = E[Nj]. Let Y, denote the number of coincidences detected in

the ith tube, where Yi is assumed to be Poisson with mean di = E[Yi] and I is the

number of tubes. It is also assumed that (1) {Nj}j' are independent, (2) {Yj}j

are independent, and (3) the probability that an emission in voxel j is detected in

the i1th tube, Pij, is known for all i,j. Under these assumptions, it can be shown

that the mean of the ith tube is given by

J

di = E xjPij (2.1)

j=1

= (Px)i, i = 1,2,...,I, (2.2)

where the (i,j)th element of P is Pij and (Pax), denotes the ith element of the

vector Pax.

Given the data, Y1, Y2, ..., yi, where y, is an observation of Yi, the goal is to

estimate x = [x1 x2 ... xj]T. In this way, an estimate for the number of photons

emitted in each voxel is obtained. Throughout the remainder of this dissertation,

x will be referred to as the emission intensity.

Shepp and Vardi used the expectation maximization algorithm [16] to compute

a ML estimate of the emission intensity. This algorithm, which is referred to as

the expectation maximization maximum likelihood algorithm, is briefly discussed

in the Appendix.

In the next section, we discuss sources of errors in PET.

2.3 Errors in Positron Emission Tomography

In reality, PET data are corrupted by errors. In this section, we go through

most of the important sources of error. Of all the errors, the one that has the

most effect on image quality is attenuation [23]. Since we focus on developing new

algorithms to correct attenuation, we provide an in depth discussion of attenua-

tion.

2.3.1 Background on Attenuation

Attenuation is the main error source in PET emission image reconstruction

[34]. The 511 keV photons produced by the annihilations are attenuated by tissue,

bone, and air. Different mediums have different attenuation effects. Before going

into the details of attenuation correction for PET, we introduce the physics behind

attenuation.

The main sources of attenuation in PET are absorption and Compton scatter-

ing. Photons are rarely absorbed because they carry 511 keV of energy. Hence,

Figure 2.3. Attenuation due to absorption

attenuation due to absorption is negligible. Fig. 2.3 is a depiction of photon ab-

sorption. Compton scattering occurs when the path of the photon is deflected by

an outer shell electron [37]. Photons in PET are more likely to undergo Compton

scattering than absorption. Most of the photons that undergo Compton scatter-

ing are not detected. Fig. 2.4 shows a photon undergoing Compton scattering.

The total effect of photon absorption and Compton scattering is referred to as

attenuation.

The survival probability along a line L, a, is the probability that a photon is

not attenuated along that line. It can be shown that the survival probability is

equal to

a = e- fL pd, (2.3)

where it is the attenuation density [29].

Since eq.(2.3) is a fundamental equation in attenuation, we provide an expla-

nation of it. Consider Fig. 2.5, which shows a line source emitting photons toward

an object with attenuation density pt. We assume that N photons incident on a

point within the object travel a short distance, Al, and only N + AN photons

survive. Obviously, AN is a negative number. We may use the following equation

to describe the above process:

AN 1

AN 1 (2.4)

Figure 2.4. Attenuation due to scattering

19

Al

-HH-

/line L

/ ...

Source Ni Nout Detector

N N+AN

Tissue

AN 1

N Ai

Figure 2.5. Basic concept of attenuation

As Al and AN go to zero, we get the differential equation

1 dN (2.5)

which can be solved by integrating along the line L. Doing so, we have that

/Nout dN

jNt = -I fPdl (2.6)

,N, N J

where Nin is the number of photons that go into the object and Nout is the number

of photons that come out of the object. From the above equation, we get

InNoit In Ni. fL J dl. (2.7)

It follows from eq. (2.7) that

Nut =e-fLdl (2.8)

Finally, eq. (2.3) follows by noting that the left side of (2.8) is the probability

that a photon along the line L is not attenuated. In the next several sections, we

discuss some other errors in PET.

2.3.2 Accidental Coincidences

If two photons from different annihilations are registered as a coincidence by

a tube, then an error occurs. This kind of error is referred to as an accidental

coincidence [22]. Fig. 2.6, shows an example of an accident coincidence. Two

annihilations happen inside the ROI. However, in both annihilations, a photon

flies out of the detector ring. The remaining photons are detected by a tube

within the coincidence window and incorrectly recorded as a coincidence.

Out of plane

Figure 2.6. Accidental coincidence in PET

The accidental coincidence rate is proportional to the concentration of the

radionuclide in the ROI. Usually, the accidental coincidences are corrected by

subtracting an estimate of them from the observed photon counts. Accidental

coincidences are estimated either by using an estimate of the single events to

compute them or using a delayed coincidence circuitry to measure them [22].

This subtraction method may produce negative counts in the corrected data,

which means that it is no longer Poisson. Consequently, several methods have

been proposed to resolve this problem [2, 28, 41]. A discussion can be found in

Anderson et al. [4] and Ollinger et al. [37]. In the next section, another error

source, scatter, is reviewed.

2.3.3 Scatter

In a scattered event, one or both of the 511 keV photons created in an an-

nihilation are scattered but they still intersect a detector ring. As a result, the

annihilation event is registered by the wrong tube. An example of an error due

to scatter is shown in Fig. 2.7.

Since photons lose energy after they are scattered, it is very easy to distinguish

them from unscattered photons by measuring the energy they pass to the detector

crystals [8]. This technique can be used to correct the scattered events by setting

a threshold to reject photons with energy sufficiently less than 511 keV. Scatter

between different slices in whole-body PET is minimized by the septa (see Section

2.1).

Figure 2.7. Scatter in PET

2.3.4 Other Errors

In reality, positrons travel a short distance (3-4 mm) before annihilating with

an electron and photons shoot off at an angle less than 180. Also, photons can

penetrate the crystal they hit and leak energy to a neighboring tube [26]. Another

type of error, called detector inefficiency, occurs when a detector fails to detect

an incident photon [17]. Fortunately, most of the time, these errors are small and

usually can be ignored or corrected.

In order to correct for attenuation, certain scans called blank scans and trans-

mission scans are used to estimate the survival probability along each tube. In

the next section, we provide a brief discussion of blank scans and transmission

scans.

2.4 Transmission Tomography in PET

In order to correct for attenuation, a blank scan and transmission scan are

performed. For the blank scan, a scan is done without a subject. The external

photon source, as seen in Fig. 2.8, rotates along the interior of the detector ring.

Let bi denote the number of observed coincidences in tube i during the blank scan,

i = 1,2,..., I. The blank scan lasts about an hour to more than 2 hours. It is

performed at the beginning of every day.

After positioning the subject in the scanner, a scan is performed, which is

referred to as a transmission scan. Like the blank scan, the transmission scan uses

a rotating external source of photons (see Fig. 2.9). The photon counts collected

in a transmission scan are the ones that survive the attenuation effect of the ROI.

Figure 2.8. Blank scan

Let ri denote the number of observed coincidences during the transmission scan

for tube i, i = 1,2,..., I. A transmission scan needs to be performed each time

a new patient is put in the PET scanner. The duration of a transmission scan

ranges from about 2 to 25 minutes. For the safety and comfort of the patient, the

transmission scans are kept as short as possible. (patient should not be exposed

to the radioactive source too long).

Intuitively, it is reasonable to estimate the survival probability by taking the

ratio of the transmission scan count per unit time and the blank scan count per

unit time. Hence, an estimate of the survival probability in tube i, &j, is

&j -- T (2.9)

where Tb and Tr are the durations of the blank scan and the transmission scan,

respectively. The data can then be corrected for attenuation by dividing the data

by the estimated survival probabilities (i.e., ML). This approach is simple and

intuitively appealing, but it does not work well if the transmission scan is too

short. One way to resolve that problem is to modify the Poisson model instead of

directly correcting the data. In next section, we introduce the modified Poisson

model.

2.5 Modified Poisson Model

In this section, we discuss an extension of the Poisson model for the case

of attenuation errors only [28]. Before discussing this model, we introduce the

concept of thinning a Poisson process.

Figure 2.9. Transmission scan

Result 1: Let N be a Poisson random variable with mean x. Suppose a random

variable X is created by thinning N with probability p. Then X is Poisson

with mean px.

Proof:

00

P(X = k) = P(X = kIN = m)P(N = m)

m=k

0(rn)pk(1 p)rn-kP(N =rm)

oo

()kl_ )mpk(_p-k_

=!(P- = m)

m=k

Se-p (l p)kk-kxm

m rk (-

k! m=k ( -k)!

Replacing m k with n, yields

P(X = k)

e-epk- (1 p)nx )+k

k!/- n!

n=O

-e-x(xp)k -,((1 --p)x)"

- e-X(Xp)ke(l-p)x

k!

k-! xp)'

~ k!

which shows that X is a Poisson random variable with mean px.

We can apply Result 1 to model the effect of attenuation in PET (see also

[41]). Let yij be the attenuation density within the j1'h voxel, where it is assumed

that the attenuation density within each voxel is constant. Let L be an I x J

matrix whose (i,j)th element, Lij, is the intersection length of the ith tube and

the Jth voxel. Then, it follows from eq.(2.3) that the probability that a photon

pair is detected by tube i is approximately equal to

-(Lit), (2.10)

O~i = C-

where (LpA), = Zl Lij pj, i = 1,2,... ,JI. Because of attenuation, the data in

tube i, Yi, is created by thinning the error free data in tube i by the survival

probability aj. By Result 1, it follows that Y, is Poisson with unknown mean

E[Y,] = ai(Px)i. (2.11)

With this modified Poisson model, several statistical reconstruction algorithms

have been proposed.

2.6 Existing Attenuation Correction Methods

In this section, we review existing attenuation correction methods including

those based on the modified Poisson model.

2.6.1 Standard Correction Method

The standard way to compensate for the attenuation is to divide the emission

data by the estimate of the survival probabilities. Specifically, the corrected data

for the ith tube, j, is given by

= --,i= ,2,...,J, (2.12)

oti

where &, is given by eq. (2.9).

This method is easy to implement and works fairly well for long transmission

scans (20-25 minutes), for short transmission scans (less than 10 minutes), it

produces noisy emission images. The conventional way this problem is addressed

is to smooth the transmission sinogram. Smoothing reduces the noise, but it also

lessens the quantitative accuracy of the images [34].

2.6.2 Reprojection Methods

Given the blank scan data and the transmission scan data, an estimate of the

attenuation density is obtained [36, 41]. Then, the estimate, it, is "reprojected"

to obtain an estimate of the survival probabilities:

= e- =1L3 ,1 = 1,2, ..., J. (2.13)

The estimated survival probabilities are used to correct the data or modify the

probability matrix, P. In addition to improving the accuracy of the estimated

survival probabilities, the transmission images provide information that can be

used to position patients, estimate scattering errors, and to provide anatomical

landmarks that cannot be seen in emission images [36].

The reprojection method provides a better estimate of the survival probabil-

ities. However, the counting noise from the transmission image introduces addi-

tional error into the emission reconstruction. In regions such as the thorax, where

the attenuation medium is heterogeneous, this problem cannot be resolved by

simply smoothing the estimated survival probabilities because it introduces bias

[13]. Since reprojection methods estimate the attenuation density, they are more

computationally demanding than the standard correction method. In the next

section, we discuss methods that simultaneously estimate the attenuation density

for survival probabilities and the emission intensity.

2.6.3 Simultaneous Estimation Methods

In [41], Politte and Snyder proposed constrained maximum likelihood recon-

struction algorithms by modifying the probability matrix under the assumption of

known survival probabilities. Their methods are not practical since the survival

probabilities must be estimated. The transmission data provides noisy estimate

of the survival probabilities because the count is pretty low in real applications.

Since the emission data are more reliable than the transmission data, in the sense

of having higher data counts, simultaneous estimation methods may provide a

more accurate estimate of the emission intensity.

In [11], a simultaneous estimation method was proposed to estimate the atten-

uation density and the emission intensities for SPECT using the Cyclic Subgradi-

ent Projection (CSP) algorithm. After the Poisson model was proposed, Fessler,

Clinthorne and Rogers proposed a joint maximum likelihood method that esti-

mates the emission intensity and survival probabilities [13, 18]. Lu proposed a

MAP algorithm to estimate the emission intensity and the attenuation density for

SPECT [32]. However, prior information about the attenuation medium must be

known.

32

To simultaneously estimate the attenuation density and the emission intensity,

we propose a new method. In the next chapter, we discuss this proposed method.

CHAPTER 3

MAXIMUM LIKELIHOOD METHOD

In this chapter, we present an attenuation correction method that jointly esti-

mates the attenuation density and emission intensity using a maximum likelihood

(ML) method. Other ML methods have been proposed [13, 18], but they have

the disadvantage of requiring more unknowns than data points. By contrast,

the proposed method has fewer unknowns than data points for typical scanners.

Specifically, it estimates 2J unknowns from I data points, where J and I are the

number of the voxels and the number of the tubes, respectively. A common choice

for the number of voxels is J = 16384 (i.e., J = 128 x 128) and, for the Siemens-

CTI ECAT EXACT 921 PET scanner, the number of the tubes is I = 73536.

3.1 Maximum Likelihood Estimator

As discussed in Section 2.2, the emission data for the ith tube, yi, is assumed

to be an observation of a random variable, Yi, where Yi is Poisson with mean

(Px)ie-(L) (recall, (Px), = xjPij and (Li)i = EJ- 1,jLij). The emis-

sion intensity, ax, and attenuation density, pi, are assumed to be deterministic but

unknown. Given the emission data, yi, y2,... ,yI, blank scan data, bl, b2, ., b,

and transmission scan data, ri,r2,... ,rj, the goal is to estimate the emission

intensity and attenuation density.

It is well known that in many estimation problems ML estimators are consis-

tent and asymptotically unbiased [44]. Thus, we propose to estimate the unknown

parameters using a ML method. In other words, we estimate the emission intensity

and attenuation density by maximizing the log likelihood function, f(ax,/ I):

[a* A*] =rmaa f (x,) ,

where f(x, it) = log P[Yi = y,Y2 = Y2,...,YI = yiIx, A], and x* and i* are ML

estimates of x and p1, respectively. The notation x > 0 means that all of the

elements of x are nonnegative. Under the model assumptions, the log likelihood

function is given by

I

f(x,p.) = logntP[Yi=YjX', ] (3.1)

i=1

Se-[(PX),e-(L) [(pe-(L)]y (3.2)

=logit----L-^-^ ---J (3.2)

= logfly

I

= E[-(pX)ie-(Lit)' + y, log(Px) yi(Lp)i log(yi!)] (3.3)

i=1

Eq. (3.1) follows from the assumption that the data is independent. Note, since

the last term in eq. (3.3) does not contain x or p., it can be ignored in the

maximization procedure. Thus, the above optimization problem is equivalent to

the following one:

(OP) [x* IA*] argwa q(,1),

where

I

O(x, A) = y_[-(Px),e-)(LL) + y, log(Px). yi(LP)i] (3.4)

t=1

The function q(x,p) will be referred to as the log likelihood function for the

remainder of this dissertation.

The problem (OP) is a difficult constrained optimization problem that must

be solved in order to obtain the ML estimate of the attenuation density and the

emission intensity. In the next section, we discuss some properties of the log

likelihood function before presenting our method for solving (OP).

3.2 Concavity of The Log Likelihood Function

Unfortunately, the log likelihood function, S(ax, ft), is not concave over the

A

feasible region, Sx,u = {(Ax, I)Ix > 0, At > 0}. However, for a fixed nonzero

vector zo, q(x, zo) and q(zo,t1z) are concave over the sets Sx = {axax > 0} and

S11 {AIAt > 0}, respectively. These results are proved below.

Result 2: Let zo be a nonzero J x 1 vector. Then, O(x, zo) is concave

over Sx.

Proof:

Our proof is based on the well known fact that a function is concave

if its Hessian is negative semidefinite [6]. Consider l(x, z0), where

Zo > 0 is a fixed J x 1 vector. It follows from eq. (3.4) that (xB, zo)

is given by

I

O(X, zo) = _[-(Px)ie-(Lz)' + y, log(Px), yi(Lz)i] ,

i=1

and its first order partial derivative with respect to xj equals

00(X, Z") Y-P-

_(a -0)- y-3j e-(Lzo)pij] (3.5)

axj j =1(p )

The Hessian of (x, zo,), which we denote by H(x), is a J x J matrix

whose (j, k)th element is given by

H()jk A0 2(X,,') (3.6)

I YiPijPik

= kPi)- ,k=l,2,...,J. (3.7)

Consider the scalar mTH(x)m, where m and x are arbitrary nonzero

vectors:

T J I y2PijmPikmk (3.8)

mTH(x)m = -SEEE ; (3.))

k=1 j=1 i=1 (PA)?

yi(Pm)?

= (px)4 (3.9)

Since for all i, y2 > 0, (Pm)2 > 0, and (Pax), > 0, it follows that

mTH(x)m < 0, which means that the Hessian is negative semidefi-

nite. Thus, O(x, zo) is concave over Sax.

Result 3: Let Zo be a nonzero J x 1 vector. Then, q(zo, Ms) is concave

over the set Sit

Proof:

Similar to the proof for Result 2, we need to prove that the Hessian of

S(zo, Iz) is negative semidefinite for zo, a nonzero fixed J x 1 vector.

From eq. (3.4), it follows that (zo, p) is given by

I

(zo, ,) = Z[-(Pzo),e-(LI)' + y log(Pz), yi(LI),]j ,

i=1

and its first order partial derivative with respect to /,j equals

9 I A I

Z( (^ = ^{Pzo)ie-(LI)'Lij 3yiLij (3.10)

i/ =1 i=1

The Hessian of (Zo, ft), which we denote by H(Is), is a J x J matrix

whose (j, k)t" element is given by

H()jk = ( A)-- (3.11)

I

S- (Pzo)iLjLike-(LA)', j, k = 1,2,...,J.(3.12)

i=1

Consider the scalar mTH(pj)m, where m and 14 are arbitrary nonzero

vectors:

J J I

mTH()m = Ej j -(Pzo)LjmjLikmke-(L)' (3.13)

k=l j=l i=l

-L(Pzo)i(Lm)e-(L) (3.14)

i=1

For all i, (Pzo)i > 0, e-(L>) > 0, and (Lm)? > 0, hence it fol-

lows that mTH(pf)m < 0, which means that the Hessian is negative

semidefinite. Thus, 4(Zo, t) is concave over S1.

3.2.1 Example to Illustrate Concavity of the Log Likelihood Function

To get a sense of the "shape" of the log likelihood function, we consider the

single voxel case. For this case, the log likelihood function becomes

q(x, P) = -pxe-" + y log(px) ylip (3.15)

It is straightforward to show that the Hessian of (x, Ii), which we denote by H,

is given by

__v ple-1p

ple-t. --pxl2e-1o

Consider the scalar mTHm, where m is an arbitrary nonzero vector:

mTHm = --rn2 pxl2e-' m2 + 2ple-mlMrni2. (3.16)

X2

From eq. (3.16), it can be seen that there is no way to guarantee that mTHm is

less than or equal to zero. For p = 1, 1 = 1, x = 40,/p = 0.049, y = 2, mi = 0.2167

and m2 = 0.01,the log likelihood function is not concave (i.e., mTHm > 0).

3.3 Optimization Procedure

We exploit the concavity results discussed in last section and propose the

following optimization procedure to obtain the ML estimates (i.e., solve (OP)):

(Step 1) Given ik-' and xk-1, determine xk by maximizing O(X,pik-l) over Sx.

In other words,

(PI) xk ,k-i1

(Step 2) Given fk-1 and xk, determine ik by maximizing O(Xk, L) over Sif. That

is,

Ak

(P2) a (k

(Step 3) Iterate between Step 1 and Step 2.

Step 1 and Step 2 are referred to as the x-step and the p-step, respectively.

The vectors xk and ,1k are the kth estimate of x and pL, respectively, where an

iteration consists of both steps 1 and 2. Unfortunately, closed form solutions for

x k and ,uk are not available, which means they must be determined using iterative

algorithms. In the next two sections, we develop algorithms for obtaining xk and

k.

A .

3.3.1 x-step

Given Ak-1l > 0, we know that O(X,/1k-l) is concave over Sa from Result 2.

Under the conditions in problem (P1), i is a solution to (P1) if and only if the

Kuhn-Tucker conditions [54] are satisfied:

90(x, -1)

it" I)x=x_ =0 if xj 0 ,j= 1, 2,...,J,

axj

oxj lx.<0O if xj=0, j=l,2,...,J.

ax3

We can rewrite the above two Kuhn-Tucker conditions as

I Yi i

(Cl) ( ) -E(L=-1)i if xj 0, j=1,2,...,J,

(C2) y i

i=1 i=1

The first condition (Cl) suggests the following fixed point iteration for solving

(P1):

(Xk)m z yP,

(,k)m+ (k) j = 1,2,..., J (3.17)

(zB )j+ =(x )jE i -(Luk-'),p'P 3~ i"''

ji F=1 e-La-)p

where (x k)7 is the mth estimate of the jth voxel of aXk and (ak)P = (a k-1) for all

j. Another possible iteration is

,--, -(L jAk-1), .

(Xk)m+ = (Xk)m. i Pi j = 1,2,...,J. (3.18)

2^---1 (P(Ck)>),

However, in our simulations, this iteration "blew up" (i.e., the iterates became

unbounded). Consequently, we use the iteration in eq. (3.17) to estimate xk.

The logic behind the fixed point iteration in eq.(3.17) is discussed below. Sup-

pose (aXk) converges to (ak)- 0 then from eq.(3.17) we have

kio kr 1-=1 / Plykixi)

%I at Pij

(^r () =^ T (pak) ) (3.19)

which implies that condition (Cl) is satisfied. In other words, we know that

the fixed point iteration in eq.(3.17) will provide a solution satisfying (Cl) if it

converges to a positive vector. Additional motivation for the technique used to

derive the algorithm in eq. (3.17) is that the EMML algorithm [28, 46] and the

ISRA algorithm [15] can be derived using it. Applying the alternating projection

method [10] to solve (P1) leads to the same iterative equation as eq. (3.17). This

iteration is guaranteed to converge by convergence results in Byrne's paper [10].

3.3.2 uI-step

From Result 3, we know that (xrk, /) is concave over S/z for xk > 0. Thus,

it is a solution to (P2) if and only if the following Kuhn-Tucker conditions hold:

I I

(C3) Z(PaXk),e-(L)'Lj = 3yiLij if /j 0, j = 1,2,...,J,

i=1 i=1

I I

(C4) Z(Pxk)ie-(L)'Lj < yiLij if pj = 0,j = 1,2,...,J.

i=1 i=1

Using the same approach that leads to the iteration in eq. (3.17), we obtain the

following fixed point iteration for solving problem (P2):

I= :l ) ,e-( L(1k ),2(prk)in3ij

1)T+1 = (l k-- t=, j = 1,2,...,J, (3.20)

,Ei=1 yiLij

where (Jtk)T is the m1th estimate of the jt" voxel of itk and (Itk) = (lk-1 )% for all

j. Another possible iterative equation is

E yiLij=

k)rn+i = (,-k)m (Zf- = (P, j= 1,2,..., J (3.21)

In our simulations, this iteration "blew up". Consequently, we estimate Ak using

the iteration in eq. (3.20).

In the next section, we investigate the behavior of the log likelihood function

as a function of the iterates.

3.3.3 Monotonically Increasing Value of the Log Likelihood Function in x-step

From the log likelihood function given in eq.(3.4) and the iterative equation in

eq.(3.17), we can prove that the log likelihood function increases monotonically

with successive iteration in the x-step. Unfortunately, this result cannot be proven

for the ji-step.

Result 5: Let z0 be a fixed nonzero J x 1 vector. Then,

0((xk)m+1, Z.) 0((Xk),, zo) > 0.

Proof:

From eq.(3.4), it is easy to see that

-((Xk)+, z) (k)m, z.)

S (p(k m+ + (o (p(Xk)m+l,)

= El-- ) )"+ (P(X)m),], +-4- E y' o (g k

,=1 i=l (P' V"' ) I

where ai = -(Lz,),.

Hence, using the iterative equation in eq.(3.17), it follows that the

second term of the above derivation can be written as

y log (p(k)m+)i

log (p(Xk)m).

I 1 Yi Pi(

where wj- = ,P "

Let 3ij = (p (a-)) =1 =

is concave, we have

g J

Yi log ij

i1 j=l

I J P.i(Xk)m+l

-- lyog ] "" )J

i=1 j=1 Pz~~

I i P p1-)

= /y'YTlog[ m

= igE [/(p( (k)-),

;=1 J= 1 j=1'J)'

I J p..({kIm

= E Y 1og -[ ( P k k]

i= 1 j= 1- (p (i ), J J '

1. Using the fact that the log function

I J

> EY 1E ij log.

=1 j=l

(3.22)

Multiplying the right side of eq.(3.22) by ( i'" the value of right

side of eq.(3.22) does not change since =' 1. Then, it follows

that

I J

EY /32j log wk

i=1 j=1

I ,ZjEi, logW

- >jy0(^) PZ.7) -

= E^E-.A)j,1j

J I 1 Yi Pi

= E xj(EZaP,) (P(21k)m), log wk

j=l i=1 ZI=1 I Pij

J I~

-=i x(Z )w l=ogw

j=l V=1

J I

E EXk(E -)(

j=l i=1

J I yiPij

= I i= (p(k)Zx)i

I J

S Z 7P=(1 k)?

i=1 .1=l

I I

= Ey y-ET i(PXk))".

i=1 i=1

Since EI Y, log EJ I k EfI y, EJ I A, log we have

Fi=l Yilg fiwj=I Z=I YEj=I Wjkw

I- log (p(k)m+l) I I ),

=ilog (p(Xk)m). j ai(P(

Then, we get

0((Xk)m+, Z.o) -((k)m, Z.)

-- EL1 y F=l[i(P(xk)m)i + a (P(xk)m), a,(P(Xk)m+l),]

Fro eq E1) a(p(Xk)m+l).

From eq.(3.4), we know that

I

i-- (p(--11)

I

=E

i=1

yi (P((Xk)m)

(P(xk)m).

I

= EY1i

i=1

(3.24)

From eq.(3.24), we obtain

I I

ZYi -- Zao(p@(k)m+l)i 0

1=1 i=1

(3.25)

Thus, we get the desired result.

((,k)m+lZo) ((Xk)m, Zo) > 0

In the next section, an interesting property of the ML estimator is presented.

(3.23)

3.3.4 Property of the ML Estimator

Result 4: Su]

?pose there exists a solution (x*, i*) to (OP). Then,

Ei=i

= e(L ) < 1 (3.

E l e- )., -

(Pa*)e-(L*)L < 1, j 1,2,... ,J. (3.

't=1 yiLij

26)

27)

Proof:

Consider the log likelihood function

I

(t,/I) = [-(Px)ie-(LA) + y, log(Px)i y.(Ljz)]

i=1

Suppose that there exists an optimal solution (ax*,/ *) and let t be

a nonnegative scalar, and x and 1A be nonnegative vectors. Then,

(x* + tax, i* + t1A) < (x*, A*). From the above relationship, we get

Ei {Yi log[ (P*'),]+-^ tyi(LA)i

-[(Pa*), + t(PX),]e-[(L/r*)'+t(Lt&)j] + (Px *)e-(L/I)} < 0

Since (Px*)ie-(LA) > 0 for all i, we have

-If orP t(PX)jl

i= loY g[ (PX ), I

L~i\ \t ig[ (Pa;*)

- ty,(LI), [(Pa*), + t(Px),]e-[(Ljr)+t(LA)j]} < 0

(3.28)

Dividing eq.(3.28) by t, we get

{ylog[l+(p*.)]-yi(Lz),-[( 7 +(Px)i]e-[(L#*)i+t(LI)']i < 0.

=(PrX)i 0

(3.29)

Now letting t -- 0 and using the L'Hoptial's rule, we obtain

{,(Px) yY(Lt)i +(Px)(L)e-(L (PX).Ie-(L'*)} < 0.

t=i (P*^x)'

(3.30)

Since eq.(3.30) must hold for any ji > 0 and x > 0, we consider the

case where x = ej (ej is a J x 1 vector of zeros except for the jh

element which is a one) and pt = 0. For this case, eq. (3.30) becomes

pj e-(L*)Pj < 0 j=1,2,...,J. (3.31)

,ti (P *i.- '

Similarly, for x = 0 and 1t = ej, eq. (3.30) yields

I[-yiL3 + (Pa,*);e-(Lw')Li] < 0, j = 1,2,..., J. (3.32)

i=1

Note, eq. (3.31) and (3.32) are only necessary conditions for (x*,/ *)

to be a maximum of (x, 1A) over the feasible region. Since x*, i* > 0,

it follows from eq. (3.31) and eq. (3.32) that

g(L .),<1, j=1,2,..., J, (3.33)

FfI (px*)ie-(LA*) iij

S=<1, j= 1,2,...,J. (3.34)

E.i=l yiLij

In the next chapter, we introduce a new least squares reprojection attenuation

correction method.

CHAPTER 4

LEAST SQUARES REPROJECTION METHOD

In this chapter, we introduce a reprojection method, where the attenuation

density is estimated using a least squares (LS) algorithm. The algorithm is an ex-

tension of one proposed by Anderson, et. al. [4] for X-ray computed tomography.

4.1 LS Algorithm for Estimating the Attenuation Density

As discussed in Section 2.4 and 2.5, let bi denote the number of photons de-

tected in tube i during the blank scan. It is assumed that bi is an observation of

a random variable, Bi, where B, is Poisson with mean Ui > 0. It is also assumed

that the number of photons detected in tube i during the transmission scan, ri, is

an observation of a random variable, Ri, where R, is Poisson with mean Uia,(r)

(recall ai == e_(L)i). Given the observed blank scan data, bl, b2, bi, and the

observed transmission data, ri, r2,..., rI, the goal is to estimate the attenuation

density, i.

It is straightforward to show that b6 and ri are the ML estimates of Ui and

Uiai(!), respectively. Consequently, bi, Ui and r, P Uiaj(r). These approxi-

mations lead to the following one: ai, ()" Taking the logarithm on both sides

of this relationship and multiplying the result by negative one, we get (Lit), z di,

where di = -log(i(1k)). Exploiting this observation, we propose to estimate the

attenuation density by minimizing the 2-norm distance between d and (Lpf):

(P3) p jzp ,

where p(At) = E'=,(di-(L'IA)I)2, tt* is the LS estimate of t, and d = [dj, d2, ..., dr].

It is straightforward to show that p(jt) is convex over Skt.

A vector jt is a solution to (P3) if and only if it satisfies the following Kuhn-

Tucker conditions:

I

(C5) ZLj(d, (Lit)i) = 0 if j $0, j = 2,...,J,

t==l

I

(C6) -L,j(di (L)i) 0 if j = 0, j = 1,2,...,J.

i=1

As before (see Chapter 3), the first condition (C5) motivates the following fixed

point iteration:

k + 1 k i = 1 *

2j c4--3 = -1,2,...,J (4.1)

,,ji E .i (L d kjL jj j .

where yj is the kth estimate of pzj. This iteration has exactly the same form as

the iterative image space reconstruction algorithm for PET [15]. Since it has been

proven that the iterative image space reconstruction algorithm converges [40], the

proposed LS algorithm converges as well. In the next section, we discuss how the

LS estimate of the attenuation density is used to correct for attenuation.

4.2 LS Reprojection Method

Given the blank scan data and the transmission scan data, we estimate the

attenuation density using eq. (4.1). Let A be the resulting estimate. Then, the

survival probabilities are estimated by reprojecting f:

ai = e- ElLjA3 ,i = 1,2,...,j. (4.2)

Next, the estimated survival probabilities are used to create a modified probability

matrix, P, where Pij = diPij. The emission intensity is then estimated by using

48

the EMML algorithm with the modified probability matrix:

EIl_ yPi,,

SP '' ( j) l I 1 2,..., J (43)

$=1 Pij

It is straightforward to show that this iteration is equivalent to

k+1 = k. l (Pak)i

Xj X4 ,' J = 1,2,'...,J. (4.4)

2-r=1 "irij

In the next chapter, simulations are used to compare the performance of the

ML method, LS reprojection method, and EMML algorithm with corrected data.

CHAPTER 5

SIMULATION

Using synthetic data and real data, we assess the proposed methods and com-

pare them to the EMML algorithm with corrected data. Dr. John Votaw, As-

sociate Professor of Radiology at the Emory University Hospital provided the

real data. This data was collected on the Siemens-CTI ECAT EXACT 921 PET

scanner. The probabilities, Pij, and line length, Lij, were determined from the

diameter of the detector ring, number of detectors, the size of ROI (the ROI is

always a square centered in the middle of the detector ring), and the number of

voxels. Specially, the probabilities are calculated using the angle-of-view method

[46]. It should be noted that the probabilities and line lengths are calculated un-

der the assumption that the detectors are contiguous. The number of detectors

for the simulated scanner and the Siemens-CTI ECAT EXACT 921 scanner are

300 and 384, respectively. The computational expense of an iteration for the ML

method, LS reprojection method, and the EMML algorithm with corrected data

are essentially the same. Consequently, the computational complexity of the al-

gorithms can be compared by considering the total number of iterations required

by each one.

5.1 Synthetic Data

For the synthetic data, we use a cardiac emission phantom (see Fig. 5.3 (a)),

chest emission phantom (see Fig. 5.7 (a)), cardiac attenuation phantom (see Fig.

5.5 (a)), and chest attenuation phantom (see Fig. 5.8 (a)). All of the phantoms are

of size 100 x 100, and the voxel size for the cardiac phantoms and chest phantoms

are 0.422 cm and 0.437 cm, respectively. The lungs, spine, and soft tissue region

in the cardiac attenuation phantom have densities 0.025 cm-1, 0.16 cm-1, and

0.096 cm-1, respectively. The low density regions (light color region) and the

high density regions (dark color region) in the chest attenuation phantom have

densities 0.035 cm-1 and 0.095cm-1, respectively. The diameter of the detector

ring for the cardiac phantoms and chest phantoms is 59.68 cm and 61.8 cm,

respectively. The ROI for the cardiac phantoms and chest phantoms is 42.2 cm

and 43.7 cm, respectively. The number of tubes is I = 44850.

We generate the emission scan, transmission scan, and blank scan data in the

following manner. The observed coincidence for the 'th tube, yi, is the observation

of a Poisson random variable with mean (P)ia), (recall a; = e-(LI)i). The

emission phantoms and attenuation phantoms specify x and it, respectively. The

observed count for the blank scan for the 2th tube, bi, and the transmission scan

for the ith tube, ri, are Poisson random variables with means U, and UiQj(3),

respectively. We specify the total count of the transmission scan to be R. and use

Ui R ( ) = 1,2,.. I.

.r"

5.1.1 ML Method Implementation of Algorithms

The term local iteration refers to the iterations used to determine x k and ik

in Step 1 (x-step) and Step 2 (y-step), respectively (see (P1) and (P2) in Section

3.3). The term global iteration refers to a single execution of both the x-step and

the j-step. For the simulations using synthetic data, we run 8 local iterations for

the x-step, 5 local iterations for the p-step, and 2 global iterations. Thus, a total

of 21 iterations and 26 iterations are used to obtain emission intensity (i.e., x2)

and attenuation density (i.e., At2), respectively.

The ML method is initialized using (x0)j = i and ()j = 0.0214,

except that e-(LMO) is replaced with a9 = (). As state in Section 3.3.1 and

Section 3.3.2, (xk) = (xk-1)j and (k) = (k-1)1.

5.1.2 EMML Algorithm with Corrected Data

For the EMML algorithm with corrected data, the iterative equation is given

below:

?+ = X k (p i (5.1)

i= (Pl k)

where y = -. It is initialized as follows: x= for all j. We run 20

iterations of the EMML algorithm with corrected data.

5.1.3 LS Reprojection Method

To initialize the iteration in eq. (4.1), we use a uniform initial estimate: j =

0.0214 for all j. We run 20 iterations to get the LS estimate of t, A. Then, we

run 20 iterations of the algorithm shown in eq. (4.4). Thus, the LS reprojection

method requires 40 iterations to get an emission image.

The total emission count, typical transmission scan count, and short transmis-

sion scan count are 1.5 million, 0.5 million and 0.125 million, respectively. The

blank scan,duration, the typical transmission scan duration, and the short trans-

mission scan duration are 60 minutes, 20 minutes, and 5 minutes, respectively. We

simulate the ML method, LS reprojection method, and the EMML algorithm with

corrected data using typical transmission scan data with and without smoothing

and the short transmission scan data with and without smoothing. In the next

section, we show the results using non-smoothed typical transmission scan data.

5.1.4 Non-smoothed Typical Transmission Scan Data

We present the results using a non-smoothed typical transmission scan data in

this section. The emission image of the cardiac emission phantom using the ML

method after the 1st global iteration, 2nd global iteration, and 3rd global iteration

are shown in Fig. 5.1 (b), Fig. 5.1 (c), and Fig. 5.1 (d), respectively. It can be

seen that the emission images improve with increasing iteration. However, after

3 global iterations, the emission images become noisy and less accurate, which

indicates that a stopping rule needs to be found. This problem is a common one

for statistically based algorithms.

The reconstructions of the cardiac attenuation phantom using the ML method

after the 1st global iteration, 2nd global iteration, and 3rd global iteration are shown

in Fig. 5.2 (b), Fig. 5.2 (c), and Fig. 5.2 (d), respectively. In the attenuation

images, the lungs are not visible, and only the contour and the spine can be seen.

Also, there are artifacts outside the ROI.

20 40 60 80 100

(a)

20 40 60 80 100

20 40 60 80 1O00

20 40 60 80 100

(b)

20 40 60 80 100

Figure 5.1. Reconstruction of cardiac emission phantom using non-smoothed typ-

ical transmission scan data (a) phantom (b) 1st global iteration (c) 2nd global

iteration (d) 3rd global iteration

100lo aj

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

100l

Figure 5.2. Reconstruction of cardiac attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) 1st global iteration

(c) 2nd global iteration (d) 3rd global iteration

.I

For the cardiac emission phantom, the results of the EMML algorithm with

corrected data, the LS reprojection method, and the ML method are shown in

Fig. 5.3 (a), Fig. 5.3 (b), Fig. 5.3 (c), and Fig. 5.3(d), respectively. From

this figure, it can be seen that the ML method yields an image with better con-

trast and resolution than the EMML algorithm with corrected data and the LS

reprojection method (note the lungs and the heart at the top center of the cardiac

phantom). It also can be seen that the EMML algorithm with corrected data

does not successfully resolve either the lungs or the heart. In the emission image

of the LS reprojection method, the spine can be seen clearly. However, the LS

reprojection method also provides incorrect information as seen by the fact that

the lungs have a higher activity concentration than the soft tissues around them

(see the phantom in Fig. 5.3 (a)).

In Fig. 5.4, we show the line plots of the 50th row of the emission images in Fig.

5.3. We observe that both the LS reprojection method and the EMML algorithm

with corrected data are more biased than the ML method (note the plot between

52 and 62, which is inside the right lung area). From this line plot, we find that

the ML method not only produces emission images with better contrast but also

gives a more accurate estimate of the emission intensity than other methods.

The EMML algorithm with corrected data does not reconstruct attenuation

images, so we only show the attenuation images of the ML method and the LS

reprojection method in Fig. 5.5 (b) and Fig. 5.5 (c), respectively. The ideal

attenuation phantom is shown in Fig. 5.5 (a). We can see that the attenuation

density reconstruction of the ML method in Fig. 5.5 (b) has a lot of artifacts and

1UU -''

20 40 60 80 100

100' --

20 40 60 80 100

20 40 60 80 100

20 40 60 80 100

,di

20 40 60 80 100

Figure 5.3. Reconstruction of cardiac emission phantom using non-smoothed typ-

ical transmission scan data (a) emission phantom (b) EMML algorithm with cor-

rected data (c) LS reprojection method (d) ML method

I

^; t' v.'. ,

." *. *" '-- .

** -

-- 'j::. ^

x 104

I I I I

SSolid: Phantom

Dash-Dash: ML Metho

Dash-Dot: LS Reproje

Dot: EM-Standard Cor

S.* .. .

.......... ................

?... ::. I

: I :

2. I

!! i ? : :! *" :

i i i/ V! ^ /- :: '. '- :

* o : l ;..

I: "*- ;; \

\,^^ / N ~: 2;.j** S -'

i V i i i,./ \ /. .

,, N.

I 1\ / / ..I I ,

d

action Method

Tection Method

I-

2$.:

I.

0 10 20 30 40 50 60 70 80 90 100

Figure 5.4. Line plot of the 50th row of the reconstruction of cardiac emission

phantom using non-smoothed typical transmission scan data

1.4

0.8

0.6

0.4

0.2

n

does not resolve the lungs, spine, and soft tissue. The LS reprojection method

resolves lungs and soft tissues but also has artifacts.

We also show the emission images of the ML method initialized with a LS

estimate of the attenuation density. Since the yt-step is supposed to determine a

better estimate of /t, we want to check if a better initial it will correct the artifact

problem of the attenuation images and provide better emission images for the ML

method. In Fig. 5.6 (b), we can observe that spine is clearer, but the lung area is

incorrect. However, the attenuation image in Fig. 5.6 (d) is better than Fig. 5.5

(b) (resolve lungs, soft tissue but no spine).

The results of the EMML algorithm with corrected data, the LS reprojection

method and the ML method using the chest phantom are shown in Fig. 5.7 (b),

5.7 (c), and 5.7 (d), respectively. From this figure, we can observe that the ML

method and the LS reprojection method produce comparable emission images.

Both images have better contrast and resolution than the EMML algorithm with

corrected data. Note, the emission image produced by the EM algorithm with

corrected data is very noisy.

We present the attenuation images of the ML method and the LS reprojection

method for the chest phantom in Fig. 5.8 (b) and Fig. 5.8 (c), respectively. The

ideal attenuation phantom is shown in Fig. 5.8 (a). For the chest phantom, the

ML method provides better contrast and more details (big and small spots in the

center) than the LS reprojection method. The LS reprojection method fails to

resolve any information.

Ip

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

Figure 5.5. Reconstruction of cardiac attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) ML method (c) LS

reprojection method

1UU ---- ---- ---- ---- ----

20 40 60 80 100

(a)

20

40-

604

80

100 .

20 40 60 80 100

1o0 o jL-

20 40 60 80 100

(b)

20

40

60 -

80

100 -...

20 40 60 80 100

Figure 5.6. ML method (a) emission phantom (b) emission image (c) attenuation

phantom (d) attenuation image

.

20 40 60 80 100

20 40 60 80 100

100, ---

20 40 60 80 100

(b)

20 40 60 80 100

Figure 5.7. Reconstruction of chest emission phantom using non-smoothed typical

transmission scan data (a) emission phantom (b) EMML algorithm with corrected

data (c) LS reprojection method (d) ML method

J A

1UU -

20 40 60 80 100

20 40 60 80 100

(b)

1UU --

20 40 60 80 100

(c)

Figure 5.8. Reconstruction of chest attenuation phantom using non-smoothed

typical transmission scan data (a) attenuation phantom (b) ML method (c) LS

reprojection method

In practice, the transmission scan data is smoothed even though it is ad hoc to

do so [34]. In next section, we show the results of the simulations using smoothed

transmission scan data.

5.1.5 Smoothed Typical Transmission Scan Data

Since the survival probability estimated by the ratio of transmission scan data

and blank scan data is usually noisy and biased, a popular way to address this

problem is to smooth the transmission scan data [13, 34, 41]. We use a 1 x 3 box

car filter, [1 1 '] to smooth the transmission scan data.

We show the emission images of the EMML algorithm with corrected data,

the LS reprojection method and the ML method for the cardiac phantom using

smoothed typical transmission scan data in Fig. 5.9 (b), 5.9 (c), and Fig. 5.9 (d),

respectively. Note that the LS reprojection method resolves more details of the

lungs, heart and the spine than the ML method. Both methods produce images

much better than the EMML algorithm with corrected data.

This study demonstrates that smoothing the transmission scan data improves

the quality of the emission images. In order to compare the accuracy of the

emission image, in Fig. 5.10, we show the line plots of the 50th row of the images

in Fig. 5.9. We can see that the line plot of the LS reprojection method and

the EMML algorithm with corrected data both provide worse estimates of the

emission intensity than the ML method.

We also compare the attenuation images. We show the result of the ML method

and the LS reprojection method in Fig. 5.11 (b) and Fig. 5.11 (c), respectively.

20 40 60 80 100

20 40 60 80 100

., ." ?. | I' ""

S' <-" *-" N' ^ r.

20 40 60 80 100

(b)

Iv|W -

20 40 60 80 100

Figure 5.9. Reconstruction of cardiac emission phantom using smoothed typical

transmission scan data (a) emission phantom (b) EMML algorithm with corrected

data (c) LS reprojection method (d) ML method

65

x 104

4.5iiiii

Solid: Phantom

4 Dash-Dash: ML Method :

Dash-Dot: LS Reprojection Method

Dot: EM-Standard Correction Method

3.5-

3

2.5

2-

1.5

2 "- .: .- .

0.5 ., :

"I A' '" :" "

1... I .

..,.J ,: ." I x.. :... .,.- i

0.5 .. t ./ -,. \: /:* ". s / : *. : :

.- .' "*.. '. .. ". '1

0. .-

0

0 10 20 30 40 50 60 70 80 90 100

Figure 5.10. Line plot of the 50th row of reconstruction of cardiac emission phan-

tom using the smoothed transmission scan data

The ML method still produces an attenuation image with heavy artifacts while

the LS method provides an image with better contrast and quality. Compared to

the attenuation images using non-smoothed transmission scan data in Fig. 5.5,

both attenuation images show more features.

The emission and attenuation image of the ML method using an initial esti-

mate for IA was obtained from the LS reprojection method are shown in Fig. 5.12

(b), Fig. 5.12 (d), respectively. The ML method produces emission image with

better contrast in the heart area in than Fig. 5.9 (b). The attenuation image

provides more information than both attenuation images in Fig. 5.11 (note the

black spine at the center below).

The emission images of the EMML algorithm with corrected data, the LS

reprojection method and the ML method of the chest phantom using smoothed

typical transmission scan data are shown in Fig. 5.13 (b), 5.13 (c), and 5.13 (d),

respectively. We can see that the LS reprojection method and the ML method

produce comparable images. Both methods provide better contrast and the image

quality than the EMML algorithm with corrected data.

The attenuation images of the ML method and the LS reprojection method

using the typical smoothed transmission data for chest attenuation phantom are

shown in Fig. 5.14 (b) and Fig. 5.14 (c), respectively. The ML method resolves

lungs, a heart and soft tissue while the LS method fails to provide any information.

The attenuation image of the ML method is smoother than Fig. 5.8 (b) but the

contrast and quality do not improve.

20 40 60 80 100

(a)

20 40 60 80 100

(b)

20 40 60 80 100

(c)

Figure 5.11. Reconstruction of cardiac attenuation phantom using the smoothed

typical transmission scan data (a) attenuation phantom (b) ML method (c) LS

reprojection method

I

*4 <"<**

60

20 40 60 80 100

20 40 60 80 100

20 40 60 80 100

1001' -

20 40 60

-I (I]

20 40 60

Figure 5.12. ML method (a) emission phantom (b) emission image (c) attenuation

phantom (d) attenuation image

4*an

I

.I4

80 100

80 100

100. .

20 40 60 80 100

20 40 60 80 100

20 40 60 80 100

20 40 60 80 100

Figure 5.13. Emission image of the chest phantom using smoothed typical trans-

mission scan data (a) emission phantom (b) EMML algorithm with corrected data

(c) LS reprojection method (d) ML method

20

40-

60

80-

100 .

20 40 60 80 100

(a)

20 40 60 80 100

(b)

1Ug"-L-

20 40 60 80 100

(c)

Figure 5.14. Reconstruction of chest attenuation phantom using smoothed typical

transmission scan data (a) attenuation phantom (b) ML method (c) LS reprojec-

tion method

In order to compare the methods numerically, we compute their peak-signal-

to-noise-ratio (PSNR) [43] for the different images. PSNR is defined by

X2k

PSNR(dB) = 10logo10 xpe, (5.2)

rose

1 E XJ j2X i

where Xpeak is the peak value of the phantom, mse = 2j=l(j Jj)2, x= is

the jth voxel in phantom, and ij is an estimate of xj. We compare the PSNR

of the methods using the typical transmission scan data in Table 5.1. We can

dB cardiac phantom chest phantom

Algorithm Standard ML LS-Rep. Standard ML LS-Rep.

Non-smooth -7.5 1112.5111 -5.92 -6.3 12.73 -4.8

Smooth -15.5 8.62 -5.62 -6.7 12.78 -4.56

Table 5.1. PSNR (dB) of emission images using typical transmission scan data

observe that the emission images of the ML method have the highest PSNR and

the EMML algorithm with corrected data has the lowest PSNR.

Another application of the proposed methods is the case of short duration

transmission scan data. In the next two sections, we show the results using a 5

minutes short transmission scan data with and without smoothing.

5.1.6 Non-smoothed Short Transmission Scan Data

In this section, we show the results using a 5 minutes transmission scan. The

emission images of the EMML algorithm with corrected data, the LS reprojection

method and the ML method of cardiac emission phantom using non-smoothed

typical transmission scan data are shown in Fig. 5.15 (b), 5.15 (c), and 5.15 (d),

respectively. From this figure, we know that the ML method resolves the heart

and lungs better than the the EMML algorithm with corrected data and the LS

reprojection method. The emission image of the LS reprojection method produces

incorrect information in the lungs region.

We also compare the attenuation images. We show the result of the ML method

and the LS reprojection method in Fig. 5.16 (b) and Fig. 5.16 (c), respectively.

We can see that the attenuation image of ML method in Fig. 5.16 (b) has a lot of

artifacts and only resolves spine (black spot at the center below). The attenuation

image of the LS reprojection method in Fig. 5.16 (c) shows the basic contours

better than the ML method but also suffers from the artifacts. Both methods

provide worse attenuation images when using short transmission scan data.

The emission and attenuation image of the ML method initialized using the

LS estimate of ft are shown in Fig. 5.17 (b), Fig. 5.17 (d), respectively. The ML

method resolves the heart area better than Fig. 5.15 (d) but produces the wrong

information of the lungs. The attenuation image provides more information and

less artifacts than Fig. 5.16 (c).

The results of the EMML algorithm with corrected data, the LS reprojection

method and the ML method of chest emission phantom using non-smoothed short

transmission scan data are shown in Fig. 5.18 (b), 5.18 (c), and 5.18 (d), respec-

tively. From this figure, we can observe that the ML method that produces

image with better contrast, resolution than the EMML algorithm with corrected

data and the LS reprojection method. The EMML algorithm with corrected data

and LS reprojection method both produce image with wrong information. The

1UU .. .

20 40 60 80 100

20 40 60 80 100

20 40 60 80 1O00

1 UUl- -

20 40 60 80 100

20 40 60 80 100

Figure 5.15. Reconstruction of cardiac emission phantom using the non-smoothed

short transmission scan data (a) emission phantom (b) EMML algorithm with

corrected data (c) LS reprojection method (d) ML method

* ,' ,,a. *

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

Figure 5.16. Reconstruction of cardiac attenuation phantom using non-smoothed

short transmission scan data (a) attenuation phantom (b) ML method (c) LS

reprojection method

0

0

0.

0

20 40 60 80 100

(a)

40

604 1^l^^y

80

100 --.

20 40 60 80 100

(c)

20 40 60 80 100

(b)

20 40 60 80 100

(d)

Figure 5.17. ML method (a) emission phantom (b) emission image (c) attenuation

phantom (d) attenuation image

I-

20 40 60 80 100

20

40

60

80

100 ',.. .

20 40 60 80 100

..^ '*:

20 40 60 80 100

(b)

20 40 60 80 100

Figure 5.18. Reconstruction of chest emission phantom using non-smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm with corrected

data (c) LS reprojection method (d) ML method

emission image of the LS reprojection method is smoother than the one of the

EMML algorithm with corrected data.

We present the attenuation images of the ML method and the LS reprojection

method in Fig. 5.19 (b) and Fig. 5.19 (c), respectively. The LS method does not

provide too much information while the ML method provides a noisy attenuation

image. However, the ML method resolves the high density region at the center.

In the next section, we show the results of the simulations using smoothed

short transmission scan data.

5.1.7 Smoothed Short Transmission Scan Data

In this section, we show the results with short transmission scan smoothed

by the same 1 x 3 box car filter, [1 1 1]. The emission images of the EMML

algorithm with corrected data, the LS reprojection method nd the ML method of

cardiac phantom using smoothed short transmission scan data are shown in Fig.

5.20 (b), 5.20 (c), and 5.20 (d), respectively. Again, we observe that We can see

that the LS reprojection method provides the comparable emission image with

the ML method in the sense of the contrast and image quality. Note that the

LS reprojection method resolves more details of the lungs, heart and spine but

also suffers from quantitative inaccuracy. Especially, the LS reprojection method

provides clearer heart area than other methods. The ML method and the LS

reprojection method both methods produce images much better than the EMML

algorithm with corrected data.

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

Figure 5.19. Reconstruction of chest phantom using non-smoothed short trans-

mission scan data (a) attenuation phantom (b) ML method (c) LS reprojection

method

4 ^An

20 40 60 80 100

20

40-

60 d

80

100 .

20 40 60 80 100

Wk

20 40 60 80 100

(b)

0

0

0.

0

20 40 60 80 100

(d)

Figure 5.20. Reconstruction of cardiac emission phantom using smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm with corrected

data (c) LS reprojection method (d) ML method

II

We show the attenuation images of the ML method and the LS reprojection

method in Fig. 5.21 (b) and Fig. 5.21 (c), respectively. The attenuation image

of the LS method provides the shape of the attenuation map while the image of

the ML method suffers from heavy artifacts.

The emission and attenuation image of the ML method initialized with fs

are shown in Fig. 5.22 (b), Fig. 5.22 (d), respectively. The ML method almost

produces the best emission images for the all short transmission scan experiments.

In particularly, it provides a very clear reconstruction of the heart area (see Fig.

5.9 (b)).

The emission images of the EMML algorithm with corrected data, the LS

reprojection method and the ML method of the chest emission phantom using

smoothed short transmission scan data are shown in Fig. 5.23 (b), 5.23 (c), and

5.23 (d), respectively. From this figure, we can observe that the ML method

produces an image with more details than the other images (Note, the small high

density spot area at the left below).

We present the attenuation images of the ML method and the LS reprojection

method in Fig. 5.24 (b) and Fig. 5.24 (c), respectively. We can see that the

ML method produces a more accurate attenuation image than the LS reprojection

method fails to resolve anything (note the two spots at the center). Compared to

the attenuation images using non-smoothed transmission scan data, attenuation

images here have more details and features.

The comparisons of the PSNR for the emission images using the short trans-

mission scan data are shown in Table 5.2. Again, the ML method has the largest

PSNR. In the next section, we show the results using the real data.

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

Figure 5.21. Reconstruction of cardiac attenuation phantom using smoothed short

transmission scan data (a) attenuation phantom (b) ML method (c) LS reprojec-

tion method

20 40 60 80 100

(a)

20 40 60 80 100

(c)

20 40 60 80 100

(b)

20 40 60 80 100

(d)

Figure 5.22. ML method (a) emission phantom (b) emission image (c) attenuation

phantom (d) attenuation image

U-

I

I

20 40 60 80 100

(a)

20 40 60 80 100

40 -

20 40 60 80 100

20 40 60 80 100

Figure 5.23. Reconstruction of chest emission phantom using smoothed short

transmission scan data (a) emission phantom (b) EMML algorithm with corrected

data (c) LS reprojection method (d) ML method

20 40 60 80 100

(a)

20 40 60 80 100

(b)

20 40 60 80 100

(c)

Figure 5.24. Reconstruction of chest attenuation phantom using smoothed short

transmission scan data (a) attenuation phantom (b) ML method (c) LS reprojec-

tion method

i tV ,

dB cardiac phantom chest phantom

Algorithm Standard 11 ML 11 LS-Rep. Standard ML LS-Rep.

Non-smooth -13.9 22.13 -9 -5.8 10.33 9.28

Smooth -13.62 19.85 -11 -15.2 9.02 -8.22

Table 5.2. PSNR (dB) of reconstruction of cardiac emission phantom using short

transmission scan data

5.2 Real Data

We obtained real data from Dr. John Votaw, Associate Professor of Radiology

at the Emory University Hospital. To create the data, a thorax phantom with

a hot heart, two lungs, a spine, and soft tissue was scanned using the Siemens-

CTI ECAT EXACT 921 PET scanner. This scanner is comprised of 24 rings of

diameter 82.5 cm, where each ring contains 384 detectors of size 6.35 x 6.35 mm2.

In addition, the diameter of the scanner's patient port is 56.2 cm and its axial

field of view is 16.2 cm. The field of view is a square where each side equals to

45 cm. There are J = 128 x 128 voxels in the field of view. In two-dimensional

mode, virtual rings are created by the cross-ring tubes. Since the ECAT EXACT

scanner consists of 24 rings and 23 virtual rings, it provides 47 planes of data.

Additional details on the scanner can be found in the paper of Wienhard et. al

[53].

The duration of the emission scan, transmission scan, and blank scan are 25

minutes, 25 minutes, and 60 minutes, respectively. Associated with each plane

are the sinograms for the emission data, transmission data, and the accidental

coincidence data from the emission and transmission scan. The sinograms are of

size 160 x 192 corresponding to 192 angles with 160 bins per angle. The number of

tube is I = 73536. In the next section, we discuss details of the implementation.

5.2.1 Overview of Implementation

All the images shown are for the 10th plane. This plane was chosen because

it contains most of the major features of the thorax phantom. The 10th plane

contains 0.288 million emission counts (i.e., E=1i yi = 2.88 x 105) and 13.82 million

blank scan counts. The number of counts for the low quality transmission scan

is 0.615 million (i.e., EL1 r, = 6.15 x 105 ), while the high quality transmission

scan contains 3.67 million counts (i.e., Eil ri = 3.67 x 106). It should be noted

that the transmission scans were automatically smoothed using the scanner's data

processing algorithm. The emission and transmission data have been corrected

for detector inefficiency and accidental coincidences using standard procedures,

however, no scatter correction has been done.

For the ML method, we run 2 global iterations, 7 and 5 iterations for the

x-step and jz-step, respectively. In other words, the ML method requires total 19

iterations and 24 iterations to get an emission image and an attenuation image,

respectively. The initial estimate for x-step is (x)j = fi or all where

a6 = l The initial estimate for ji-step is (p)j = 0.07 for all j. The initial

bi rr "

estimate for the survival probability is a = where rb = 60 minutes and

Tr = 25 minutes.

Twenty iterations of the EMML algorithm with corrected data method is initial-

ized using x = for all j, where = 4. For the LS reprojection method,

to initialize eq. (4.1), we use a uniform initial estimate: f0 = 0.07 for all j. We

run 20 iterations to get he LS estimate of ft, A. Then, we run 20 iterations of

the algorithm shown in eq. (4.4). Thus, LS reprojection method requires total 40

iterations to get an emission image.

The transmission scan is a high quality scan because it was performed long

after the emission scan was done (no emission counts will be mis-registered as

transmission counts).

5.2.2 Results

The emission images of the filtered back-projection method, the EMML algo-

rithm with corrected data, the LS reprojection method, and the ML method are

shown in Fig. 5.25 (a), Fig. 5.25 (b), Fig. 5.25 (c), and Fig. 5.25 (d), respec-

tively. From this figure, we observe that the EMML algorithm with corrected

data produces an image with clearer lungs boundaries of the other methods. Also,

the ML method produces comparable but smoother image than the EMML algo-

rithm with corrected data. The filtered back-projection method provides a noisy

emission image with lots of artifacts everywhere.

Except the emission images of the filtered back-projection method, all emission

images have a "ring" artifact outside the phantom (see Fig. 5.25 (c) and Fig. 5.25

(d)), which is not visible in the image using the EMML algorithm with corrected

data because the limited gray level scales of the image (the heart area in Fig. 5.25

(b) has high intensity values compared to the background, so the "ring" artifacts

on the background cannot be seen).

20 40 60 80 100 120

(a)

20 40 60 80 100 120

(C)

20 40 60 80 100 120

(b)

20 40 60 80 100 120

(d)

Figure 5.25. Emission images using high quality smoothed transmission scan data

(a) Filtered back-projection method (b) EMML algorithm with corrected data (c)

LS reprojection method (d) ML method