Citation
Maximum likelihood estimation of detector efficiencies in positron emission tomography

Material Information

Title:
Maximum likelihood estimation of detector efficiencies in positron emission tomography
Creator:
Lee, Wen-Hsiung
Publication Date:
Language:
English
Physical Description:
vii, 100 leaves : ; 29 cm.

Subjects

Subjects / Keywords:
Data models ( jstor )
Estimation methods ( jstor )
Fisher information ( jstor )
Geometric angles ( jstor )
Image reconstruction ( jstor )
Pets ( jstor )
Photons ( jstor )
Positrons ( jstor )
Random variables ( jstor )
Simulations ( jstor )
Dissertations, Academic -- Electrical and Computer Engineering -- UF ( lcsh )
Electrical and Computer Engineering thesis, Ph.D ( lcsh )
Image processing -- Digital techniques ( lcsh )
Imaging systems in medicine -- Mathematical models ( lcsh )
Tomography, Emission ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph.D.)--University of Florida, 2001.
Bibliography:
Includes bibliographical references (leaves 96-99).
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by Wen-Hsiung Lee.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
028221204 ( ALEPH )
49240595 ( OCLC )

Downloads

This item has the following downloads:


Full Text










MAXIMUM LIKELIHOOD ESTIMATION OF DETECTOR EFFICIENCIES IN
POSITRON EMISSION TOMOGRAPHY


a


By


WEN-HSIUNG LEE


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


2001


























Copyright 2001



by



Wen-Hsiung Lee














ACKNOWLEDGMENTS


First I would like to thank my advisor and mentor, Professor John M. M. Anderson, for his endless effort and support on my research. His encouragement and suggestions have helped me through all the difficulties. More importantly, he taught me the methodology of solving problems, and he showed his diligence and patience in research. In short, he set a good example to me, which I believe will benefit me for all my life.

Sincere thanks go to my committee members, Dr. Edmonson, Dr. Li, and Dr. Pardalos, for their guidance and patience. Sincere thanks also go to Dr. Vatow for his constructive comments and invaluable help.

I also feel thankful toward Hong and Kerkil, my colleagues, who made me enjoy a pleasant team-working environment.

I owe a large measure of gratitude to my mother and my sisters for their unwavering support to help me go through tough times. Without them, I would not have had the strength to go on with my study at UF.


iii














TABLE OF CONTENTS


ACKNOWLEDGMENTS . . . . . . . . . . . .

ABSTRACT...... ..................

CHAPTERS

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


1.1 Overview of PET ......... 1.2 Basic Physics of PET ........ 1.3 Errors in PET . . . . . . . . . . .
1.3.1 Non-Collinearity . . . . .
1.3.2 Positron Range . . . . . .
1.3.3 Attenuation . . . . . . . .
1.3.4 Accidental Coincidences .
1.3.5 Scatter . . . . . . . . . . .
1.3.6 Detector Deadtime . . . .
1.3.7 Detector Inefficiency . . . 1.4 Mathematical Model for PET . .
1.4.1 Poisson Model for Emissio 1.4.2 Modified PET Model for E
1.5 Existing Methods for Estimating 1.6 Outline of Research . . . . . . . .


. . . . . . . . . . . . . . . . . . . . 2
. . . . . . . . . . . . . . . . . . . . 4
. . . . . . . . . . . . . . . . . . . . 5
. . . . . . . . . . . . . . . . . . . . 5
. . . . . . . . . . . . . . . . . . . . 5
. . . . . . . . . . . . . . . . . . . . 6
. . . . . . . . . . . . . . . . . . . . 6
. . . . . . . . . . . . . . . . . . . . 6
. . . . . . . . . . . . . . . . . . . . 7
. . . . . . . . . . . . . . . . . . . . 7
. . . . . . . . . . . . . . . . . . . . 8
ni Scan Data .. .. . ... . . ... 8
mission Scan Data ... . .. . .. 9
Detector Efficiencies .. .. . . ... 10 . . . . . . . . . . . . . . . . . . . . 13


2 ML ESTIMATION ALGORITHMS USING A POISSON DATA MODEL.
2.1 Poisson Model I ....................................
2.1.1 Estimation Problem Associated with Poisson Model I .....
2.1.2 Estimation Algorithms associated with Poisson Model I . . . .
2.2 Poisson M odel II . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Estimation Problem Associated with Poisson Model II . . . .
2.2.2 Estimation Algorithms Associated with Poisson Model II . . .

3 ML ESTIMATION ALGORITHMS USING A REFINED DATA MODEL.
3.1 Poisson Model III ...................................
3.1.1 Estimation Problem Associated with Poisson Model III . ...
3.1.2 Estimation Algorithms Associated with Poisson Model III . .


iv


iii vi


16 16
21 21 30 33
34

38 38
44 44









4 SIMULATION . . . . . . . . . . . . . . . . . . . . . .
4.1 Synthetic Data . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Comparison of Detector Efficiency Estimates . . . .
4.1.2 Comparison of Reconstructed Emission Images . . .
4.2 R eal D ata . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . .

6 FUTURE WORK . . . . . . . . . . . . . . . . . . . .

APPENDICES

A CONCAVITY OF THE LOG-LIKELIHOOD FUNCTIONS


A.1 Poisson Model I


A.2 Poisson Model II . . . . . . . .
A.3 Poisson Model III . . . . . . . .

B CRAMER-RAO LOWER BOUNDS
B.1 Poisson Model I . . . . . . . . .
B.2 Poisson Model II . . . . . . . .
B.3 Poisson Model III . . . . . . . .

C DERIVATION OF THE E STEP OF
C.1 Poisson Model I . . . . . . . . .
C.2 Poisson Model II . . . . . . . .
C.3 Poisson Model III . . . . . . . .

REFERENCES... ...........

BIOGRAPHICAL SKETCH . . . . .


THE EM ALGORITHMS . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .


v


48 49 50
54 56


. . . . . 75


77


78 78 81 82

84 84 86 88

89 89 92
94

96

100


. . . . .


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


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














Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

MAXIMUM LIKELIHOOD ESTIMATION OF DETECTOR EFFICIENCIES IN POSITRON EMISSION TOMOGRAPHY By

Wen-Hsiung Lee

December 2001
Chairman: John M. M. Anderson
Major Department: Electrical and Computer Engineering

Positron emission tomography (PET) is a medical imaging modality that enables physicians and researchers to study biochemical process within the human body and small animals. In PET, a subject is administered a radiopharmaceutical that is absorbed preferentially by the region-of-interest (e.g., brain). Due to the radiopharmaceutical and certain interactions at the cellular level, photon pairs are emitted in all directions and the number of photon pairs emitted from a region is proportional to the metabolic rate of the region. PET scanners are designed to detect and count the photon pairs. From this photon count data, images are reconstructed whereby the value of a pixel is proportional to the metabolic rate of the region associated with the pixel. The reconstructed images provide valuable information to help physicians detect abnormalities, such as tumors, which have extremely high metabolic rates.

There are many sources of error in PET that must be addressed in order to obtain accurate images. The focus of this research is the problem of detector inefficiency. Due to inherent detector non-uniformity, detector efficiencies are less than 100%. To obtain more precise emission images, the photon count data must be corrected to account for the effect of detector inefficiency.


vi








We introduce a maximum-likelihood method for estimating detector efficiencies in PET. First, we develop three Poisson models for blank scan data obtained from rotating rod sources. Then, for each model, we estimate the detector efficiencies using expectation-maximization algorithms, where the maximization step is solved using two optimization algorithms. As desired, the resulting algorithms have the property that the log-likelihood function is non-decreasing as the iteration number increases. For each data model, one of the proposed algorithms guarantees that the efficiency estimates always lie in the interval [0, 1]. Although the second algorithm for each data model does not have this property, in all of the experiments performed it produced efficiency estimates that were between zero and one. Simulation studies using synthetic data demonstrate that, based on various comparison criteria, the proposed estimation algorithms outperform two alternative approaches. Additionally, simulation studies using real data show that the third Poisson model leads to much better reconstructed emission images than the first two models.


vii














CHAPTER 1
INTRODUCTION


Medical imaging modalities are widely used to obtain anatomical and biochemical information about the human body. For example, X-ray computed tomography and magnetic resonance imaging are used to acquire detailed images of the anatomical structure of a region of interest (ROI) (e.g., brain, heart, thorax). In many situations, however, it is more important to obtain functional information than anatomical information [1]. For instance, the knowledge of blood flow rates, oxygen utilization rates, and glucose utilization rates of brain or heart might be imperative for diagnosing certain diseases.

Biochemical and functional information can be obtained via a class of imaging techniques known as nuclear medicine imaging. Among the techniques within nuclear medicine imaging are planar imaging [2, 3], single-photon emission computed tomography (SPECT) [4-6], and positron emission tomography (PET) [7].

Planar imaging suffers from relatively low availability of radiopharmaceuticals, poor sensitivity of cameras, and, most importantly, the so-called projection problem. The projection problem refers to situations where the ROI might be obscured due to activity in front of or behind the ROI. Moreover, photons originating from the ROI might be attenuated due to overlying tissue.

Although SPECT overcomes the projection problem and provides more sensitivity than planar imaging, it still has insufficient sensitivity and quantitative accuracy for certain applications. Furthermore, SPECT uses radiopharmaceuticals that incorporate long half-life radionuclides, such as 121 and 201TI. Because of the long half-lives


1






2


(ranging from hours to several days) of the radionuclides, patients are malignantly exposed to radioactivity for a long period of time.

PET avoids those afore-mentioned shortcomings and provides better images in the sense of efficiency, sensitivity and accuracy. In addition, positron-emitting isotopes of carbon, oxygen, nitrogen and fluorine can be readily incorporated into many useful radiopharmaceuticals because these compounds are naturally used by the body [8]. However, because radionuclides incorporated into radiopharmaceuticals for PET have short half-lives (2 minutes, 10 minutes, 20 minutes, and 110 minutes for 1C, 15o, 13N and 18F, respectively), on-site cyclotrons are needed to produce radionuclides. The high cost and high complexity of cyclotrons have prevented PET from dominating the field of nuclear medical imaging. Even so, PET still has a wide range of use in clinical and research applications due to its advantages over other modalities.

There are a number of error sources in PET, such as attenuation, accidental coincidences, and detector inefficiency. In this dissertation, we address the problem of detector inefficiency and propose algorithms, based on the maximum likelihood estimation principle, that estimate detector efficiencies.


1.1 Overview of PET

A PET study begins with the injection or inhalation of a radiopharmaceutical by a subject. A radiopharmaceutical is a substance that is labeled with a positronemitting isotope. The radiopharmaceutical is transported to and absorbed by the ROI. When the isotope decays, it emits positrons, which travel a short distance (on the order of tenth of millimeter) before annihilating with electrons. The areas with a high concentration of the radiopharmaceutical emit more positrons than the areas of low concentration. The annihilation of a positron with a nearby electron produces two high-energy (511 keV) photons that move away from each other in nearly opposite directions. The position where the annihilation takes place is unknown and






3


the photons' line of flight is omni-directional and completely random. If these two photons, after traveling through the subject, are detected by a pair of detectors on a detector ring surrounding the ROI within a small timing window r (~ iOns), a coincidence is recorded by that detector pair or tube (a detector pair is also referred to as a tube.) The observed emission data are the collection of the counts (i.e., total number of coincidences) recorded by each detector tube. An illustration of the scanning process is provided in Fig. 1.1, which shows a cross-section of the scanner.


Photon Pair




Tube

ROI _,Coincidence Dete 1or -g
Window (,c)



Positron
Annihilatio Location
Position


Figure 1.1: Annihilation of a positron with a nearby electron produces a pair of photons that propagate in nearly opposite directions.


The number of photons emitted from a region is proportional to the concentration of radiopharmaceutical in the region. Moreover, the concentration of the radiopharmaceutical in the region is proportional to the metabolic rate of the region. Using the emission data, reconstruction algorithms create images, called emission images, where the value of a pixel is proportional to the concentration of the radiopharma-






4


ceutical in the region associated with the pixel. Thus, emission images can be used to assess biochemical processes.


1.2 Basic Physics of PET

The cylindrical rings of detectors surrounding the ROI typically have a diameter of 80-100 cm and an axial extent of 10-20 cm. The detectors are shielded from outside the field of view (or ROI) by relatively thick, lead end-shields. Most scanners can switch between two operating modes. In the first mode, known as slice-collimated mode or 2-dimensional (2D) mode, axial collimation is provided by thin annular rings of tungsten called septa, and only direct plane sinograms involving coincidences between detectors on the same ring are recorded. In practice, the situation may be more complicated in that, in order to improve statistics, direct and cross plane sinograms can be combined in different ways [9]. The other mode is fully three-dimensional (3D) mode where septa are retracted and cross plane sinograms of coincidences between detectors on different rings are also recorded [7].

Detectors are critical components of a PET scanner [7]. Each detector ring usually contains several hundred detector modules, where each module consists of an array of crystals coupled to several photomultiplier tubes. When a photon interacts with a crystal, electrons are moved from the valence band to the conduction band. Then, the electrons return to the valence band at impurities in the crystal. This process causes numerous light photons to be emitted. A photomultiplier tube collects light photons and converts them into electrical signals, which are decoded so as to determine which detector module was actually hit by a photon. The material bismuth-germanate is widely used as the crystal because of its high attenuation density, high light output (2500 light photons per 511 keV photon), fast rise time, and short decay time (300 ns). Crystals with high attenuation densities are desirable because the probability of a photon interaction is very high. Moreover, a high light output improves positioning






5


accuracy, while a fast rise time is helpful for accurate timing, and a short decay time enables the scanner to handle high counting rates.


1.3 Errors in PET

There are two major categories of physical effects in PET: (1) effects due to the positron annihilation process and propagation of photons through the body, such as positron range, attenuation, Compton scatter, and accidental coincidences, and (2) effects that result from the interactions of photons with detectors, such as detector efficiency and dead-time [7].


1.3.1 Non-Collinearity

As mentioned previously, an emitted positron interacts with an electron as it travels through the body. Each interaction requires certain energy. When the momentum of the positron's energy is nearly zero, annihilation occurs and two annihilation photons are produced each with an energy of 511 keV. The photon pair flies off along nearly collinear paths, where the degree of non-collinearity depends on the energy left in the positron [7]. The energy of the positron after annihilation varies widely among isotopes. The divergence from collinearity is on the order of a degree or less. So it is usually ignored.


1.3.2 Positron Range

The distance the positron travels before annihilation is termed positron range. Typically, the positron range is about 3-5 Mm [7]. Since the resolution of reconstructed PET images is over 5 mm for most commercial scanners, the error introduced by positron range is not a serious problem. However, for high resolution experimental scanners, this becomes the fundamental limitation of PET.






6


1.3.3 Attenuation

Two interactions that may cause a photon pair to go undetected are photoelectric absorption and Compton scattering. For 511 keV photons, the incidence of photoelectric absorption is very low and can be ignored in most cases [7]. A Compton scatter happens when a photon interacts with an outer shell electron. In this case, the photon is deflected from its path. Consequently, the counts of the tube corresponding to the original path of a photon pair does not increase when one or both of the photons undergo Compton scattering. The total effect of photoelectric absorption and Compton scatter is called attenuation and undetected annihilations are termed scattered events.


1.3.4 Accidental Coincidences

If two photons resulting from different annihilations are detected within the coincidence timing window, they will be incorrectly interpreted as a coincidence event. Such events are called accidental coincidences or randoms. as true coincidences or trues.

The rate of accidental coincidences is related to the concentration of the radionuclide in the ROI. A common correction strategy is to subtract an estimate of the randoms by using an estimate of the single events or delayed coincidence circuitry [10].


1.3.5 Scatter

Suppose one or both photons of a photon pair undergo Compton scattering. If the photons remain in the detector plane and are detected, then the annihilation event is recorded by the wrong tube. In other words, since the photons' path is not collinear, the event is mis-positioned and additive noise results. However, due to the energy loss caused by the scattering process, scattered photons can be distinguished from






7


unscattered photons to a limited degree by setting a threshold such that photons with sufficiently low energy, which correspond to scattered photons. are ignored.


1.3.6 Detector Deadtime

The time required to process a single photon event limits the counting rate of a PET scanner [11]. Event processing begins with the rising edge of the pulse for the first detector involved. The pulse is integrated for some time interval, after which position calculations and energy discriminations are performed. The detectors are "dead" to new photon events during this time. Since the number of randoms increase as the square of the activity in the field of view [7], the number of detected counts may decrease as the flux of incident photons increases. This effect is accounted for in two ways. First, it becomes a dominant limitation on the injected dose. Second, a scale factor is used to correct data for the deadtime.


1.3.7 Detector Inefficiency

The efficiency of a detector pair is the probability that an incident photon pair is detected by the detector pair. Due to inherent detector non-uniformities, detector efficiencies are not equal and less than 100% [12,13]. The efficiency of current PET detectors depends on a number of geometric and non-geometric factors. The nongeometric factors include fluctuations in photomultiplier tube gains and variations of crystal characteristics [7]. Geometric factors include the angles subtended by detectors and the angles of incidence of photons. The angle subtended by a detector pair depends only on the perpendicular distance from the origin, and the angle subtended by the detector pair affects the angles of incidence of photons. And the length of intersection of a detector and the path of a photon depends on the angle of incidence of the photon. If the intersection length is larger, the probability that the photon is detected by the detector is higher.






8


1.4 Mathematical Model for PET


Since the data in PET follow from a counting process, it is assumed that the emission process is a spatial Poisson process. In 1982, Shepp and Vardi [14] proposed a Poisson model for PET that has been widely accepted as a valid mathematical model for the PET process. However, their model was based on the assumption that the data are error-free except for Poisson noise. In the following, the basic model for PET emission data is first introduced. Because this dissertation focuses on how to correct (or normalize) emission data for variations in detector efficiency, the basic data model is then modified to take into account the effect of detector inefficiency.


1.4.1 Poisson Model for Emission Scan Data

The ROI is divided into a grid of J voxels and the emission intensity is assumed to be constant over each voxel. Assuming the data are free of errors, the emission datum recorded by detector pair (k, 1) is a realization of a random variable Ykl:


YkI ~ Poisson([PX]kl), k = 1, 2, ... , I, I E fank,


where

[PX]kl~ Pkl,jXj,11
j=1

and where I is the number of detectors on a ring and the set fank is the set of detectors that are paired with detector k (see Fig. 1.2). In addition, Pkl,j is the probability that a photon pair emitted from voxel j is incident on detector pair (k, 1), xj is the mean number of photon pairs emitted from voxel j, and x = [xi, x2,..., xJ]T is the emission intensity vector.






9


Figure 1.2: Definition of fank.

1.4.2 Modified PET Model for Emission Scan Data

Again, the efficiency of a detector pair is the probability that an incident photon pair is detected by the detector pair. Let Ckl be the efficiency of detector pair (k, 1). Having a detector pair that is not 100% efficient is mathematically equivalent to thinning [15] the incident photons. To see what effect the thinning has on a Poisson random variable, we make use of the following well-known Corollary. Corollary Let Z be a Poisson random variable with mean 2. Suppose a random variable Y is created by thinning Z with probability p. Then, Y is Poisson with mean pZ.


Proof:


P(Y =y) = P(Y= y I Z = z)P(Z = z)


= (')py(1 - p)Z-y e-z(


= Z ! y ! ( . p ) Z . ( )z
Z=Y y!(z - y) Z!






10


_- 2P* (1 P)z-Y(Z)z y. ZY (z -y)!

Using the change of variable z' = z - y, the probability P(Y = y) becomes


P(Y Y) - ePY fI~ (l-P)Z'(Z)'+Y
P( )=e - E 1
y z'=0
= (Pz)y ((1 -- P))/ z'=O
y'- ___==e- 2 (P2) _,4P)
y!
e-pz(p2)y 12
(1.2

Thus, Y is a Poisson random variable with mean p2.

Q.E.D.


According to the above Corollary, the emission data, after thinning with probability Ekl, become

YkI ~ Poisson(Ek [PX]kl), k = 1, 2,. . ., I, I E fank. (1.3) Given the observed emission scan data y, which is a vector with elements {ykl}, the problem is to reconstruct the emission intensity vector x. Shepp and Vardi [14], the developers of the Poisson model for emission scan data, assumed that the detectors were 100% efficient. However, detector efficiencies must be estimated before emission images can be reconstructed, because, in reality, they are not ideal.


1.5 Existing Methods for Estimating Detector Efficiencies


To estimate a set of detector efficiencies, a scan is taken using either rotating rods or a uniform cylinder that contains positron emitting nuclides. The scan, which is referred to as a blank scan (no patient in the scanner), lasts about 2-4 hours and is usually performed once a month (see Fig. 1.3).






11


Figure 1.3: Blank scan.

Numerous methods [9,16-24] have been developed for estimating detector efficiencies. Unfortunately, existing methods for estimating the detector efficiencies are largely ad hoc and do not take into account the statistical properties of the blank scan data. For instance, an early technique discussed in [9] was to estimate the efficiency of a detector pair by computing the ratio of the number of photon pairs measured by that detector pair and the average number of photon counts per tube. It has been noted that the accuracy of this heuristic technique is limited by the statistical accuracy of the blank scan [9]. Another technique for estimating detector efficiencies is the well known fan-sum method by Hoffman et al. [16]. Hoffman et al. assumed that bkl, the number of photon pairs recorded by detector tube (k, 1), is equal to


bkl Ek--: inki, k = 1, 2, .. . , I, I E fank, (1.4)






12


where Ck is the efficiency of the k"h detector and nki is the unknown number of photon pairs transmitted along detector tube (k, 1). Hoffman et al. also assumed that



E finkl constant. (1.5)
lEfan,

Their reason for making the assumption in (1.5) is because the summation is over a large number of detectors and the variations would tend to average out. Under these two assumptions, the proposed estimate for the efficiency of the kth detector


ek E bkI (1.6)
lEfan,

is approximately equal to

ek - Ek . constant. (1.7)

The approximation in (1.5) is poor when there are large fluctuations in the detector efficiencies [17] or the number of terms in the set fank depends on k. In such instances, the fan-sum algorithm may not provide accurate estimates for the detector efficiencies. Also, the fan-sum algorithm does not account for the Poisson nature of the blank scan data.

Several methods [9,16-19] have been developed based on the assumption of noisefree (i.e., no Poisson noise) blank scan data. For example, an iterative method by Ferreira et al. [17] is based on the assumptions that the data are noise free and every detector views the same activity:



E nla = A, k 1,2,...,I, E fank. (1.8)
tEfank






13

Substituting (1.4) into (1.8) gives


k =1, 2,. ,I, I E fank. (1.9)
IEfank El

To get an approximate solution for the above system of equations, the following fixed-point iteration is used:


m+1 bkI
Ek =-E , k = 1, 2,., I, 1 E fank, (1.10)
lEfank

where c is the mth estimate of El. The constant factor 1/A is neglected because the efficiency estimates are normalized so that their mean is one. Like the fan-sum algorithm, Ferreira's algorithm assumes noise-free data and relies on a key assumption (see (1.8)) that is questionable when the number of terms in the set fank depends on k.

Another algorithm based on the assumption of noise-free data is the algorithm by Defrise et al. [9]. The algorithm by Defrise et al. parallels the fan-sum method with the exception that geometric means rather than arithmetic means of the data are used to estimate the mean number of unobserved number of transmitted photon pairs. Defrise et al. claimed that their algorithm produces unbiased efficiency estimates in the sense that they converge to the real efficiencies in the ideal case where the blank scan data are noise-free.


1.6 Outline of Research


All of the estimation methods we have discussed are sub-optimal in the sense that they assume noise free data and do not provide ML estimates. To capture the stochastic nature of PET data, we propose three Poisson models for blank scan data. At first we ignore the problem of detector penetration [25-28] when developing the






14


first two Poisson models for blank scan data. However, later we incorporate the effect of detector penetration and refine our original Poisson models for blank scan data.

Motivated by the optimality properties of ML estimators [29], we maximize the log-likelihood function associated with each Poisson model by using the expectationmaximization (EM) algorithm [30, 31]. To apply the EM algorithm, we choose the blank scan data and unobserved number of transmitted photon pairs for each detector tube to be the complete data. Two possibilities have been explored:

" Over-parameterized approach, where the unknowns are the efficiencies of the

detector pairs, {Ekl}, and certain emission means.

" Sufficiently-parameterized approach, where unknowns are the efficiencies of individual detectors, {Ek}, and certain emission means.

Because the over-parameterized approach does not take into consideration the assumption that EkI = EkEh, k = 1, 2, ... , I, 1 E fank, the number of unknowns is much higher than that for the sufficiently-parameterized approach. Closed form expressions for the E step and M step of the EM algorithm are derived for our first model using the over-parameterized formulation. However, simulations show that the over-parameterized approach has poor performance. Thus, the results for the over-parameterized approach will not be presented in Chapter 4, which contains our simulation results.

For each of our three proposed Poisson models for blank scan data, we use the sufficiently-parameterized formulation and develop two algorithms for estimating detector efficiencies. For each data model, a closed form expression for the E step of the EM algorithm is derived. Unfortunately, no closed form solutions to the maximization problems in the M step are available. For each data model, we propose solving the maximization problems in the M step iteratively by using a fixed-point iteration or the method of coordinate descent. All the resulting (EM) algorithms have






15


the property that the log-likelihood function increases monotonically with increasing iteration. It is desirable that a detector efficiency estimator provides estimates that lie in the interval [0, 1]. The coordinate descent based algorithms we propose for each data model guarantee that the efficiency estimates always lie in the interval [0, 1]. The fixed-point iteration based algorithms for each data model do not have this property. However, in all the experiments we performed the efficiency estimates did lie in the interval [0, 1]. In terms of certain objective criteria, simulation studies using synthetic data demonstrate that the proposed estimation algorithms outperform the fan-sum and Ferreira's algorithms. Additionally, our experiments show that emission images based on the proposed detector efficiency estimates have smaller mean squared error than emission images obtained using the detector efficiency estimates generated by the fan-sum and Ferreira's algorithms. Simulation studies using real data show that the estimation algorithms based on the third Poisson model provide much better reconstructed emission images.

This dissertation is organized as follows. In Chapter 2, we introduce two Poisson data models for blank scan data and develop estimation algorithms that provide ML estimates of the detector efficiencies. Then, in Chapter 3, we refine our data model to incorporate the effect of detector penetration and present the corresponding ML estimation algorithms. In Chapter 4, we discuss the results of our simulation studies. In Chapter 5, conclusions about our research are given. Finally, we suggest some topics for future work in Chapter 6.














CHAPTER 2
ML ESTIMATION ALGORITHMS USING A POISSON DATA MODEL


In this chapter, we first present two Poisson models for blank scan data obtained from rotating rod sources. Then, we present ML estimation algorithms based on these two Poisson models. The first model is based on the fact that the mean number of photon pairs transmitted along a detector pair only depends on its perpendicular distance from the origin, which is the center of field-of-view. The second model, which is a modification of the first model, takes into consideration that the rods are nearly uniform. Note, the second model has fewer unknowns than the first model.

In many scanners, the space between blocks of detectors is filled up with lead. Throughout our research, though, we make the simplifying assumption that there is no lead between blocks of detectors. Also, we assume the detector rings are perfectly circular.


2.1 Poisson Model I

Before discussing the first Poisson model, we introduce some notation. Let a denote the angle between the lines that connect the origin to detectors k and 1, and p(k, 1) denote the perpendicular distance from the origin to the line defined by the mid-points of detectors k and 1 (see Fig. 2.1). The detector indices k and I are defined such that k < 1, and the angle a is defined counterclockwise starting from detector k and ending at detector 1. Clearly, p(k, 1) = r cos(2), where a = (1 - k) 2I
and r is the radius of the detector ring. Therefore, for detector pairs (ki, 11) and (k2, 12), p(ki, 1I) = p(k2, 12) if 11 - k, = 12 - k2. In addition, p(ki, 1i) and p(k2, 12) will


16






17


k



Figure 2.1: Definitions of p(k, 1) and a.

be equal if (11 - k1) + (12 - k2) equals the number of detectors on one ring. To see this fact, consider Fig. 2.2, where detector pairs (1,114) and (50,321) of the SiemensCTI ECAT EXACT 921 scanner [26] are depicted (note: this scanner has I = 384 detectors per ring). Because a, = 1(114 - 1) and a2 = 2'4(321 - 50), their sum


Figure 2.2: An example to demonstrate how p(k, 1) is related to the value of 1 - k. equals a1+a2 = r(113+271) = 27r. Clearly a2+ = 27r as well. Therefore, a, = /, which implies p(l, 114) = p(50, 321).

For the Siemens-CTI ECAT EXACT 921 scanner, Table 2.1 illustrates the detector pairs of the first projection and the perpendicular distance of each pair. From the


271 50






(D321


114 11 3



"7-1 a
"--\ -----------






18


Table 2.1: Perpendicular distances of detector pairs for the first projection of the Siemens-CTI ECAT EXACT 921 scanner.


1st projection
member pair (k, 1) 1 - k p(k, l)(cm)
1st (152,264) 112 dl = 25.1114
2nd (152,265) 113 d2 = 24.8428
3rd (151,265) 114 d3 = 24.5726

80th (113,304) 191 d80 = 0.3375
81st (112,304) 192 d81 = 0
82nd (112,305) 193 d8o = 0.3375

159th (73,343) 270 d3 = 24.5726
160th (73,344) 271 d2 = 24.8428



discussion above, it can be seen that the value of I - k determines the perpendicular distance p(k, 1). Taking this property into account, it can be seen that there are 81 possible perpendicular distances from the origin.

To see how the number of transmitted photon pairs is related to p(k, 1), consider a scanner with 16 detectors, where the dashed circle stands for the path of the rods (see Fig. 2.3). It is easy to show that the arc of tube 2 is longer than the arc of tube 1. Thus, because the rods rotate at a constant speed, the rods will stay in tube 2 for a longer time than in tube 1. Observe that p(k, 1) of tube 2 is larger than p(k, 1) of tube 1. More generally, it can be proved using a simple geometric argument that the dwell time is proportional to the perpendicular distance from the center to the tube. And, because the activity viewed by a detector tube is proportional to the dwell time of the rods in that tube, it is reasonable to claim that the mean number of the transmitted photon pairs depends only on the perpendicular distance from the origin to the tube. For a blank scan, let BkI denote the observed number of photon pairs detected by detector pair (k, 1), and Nk, denote the unobserved number of transmitted photon pairs along detector pair (k, 1). Given the discussion above, we assume that






19


2











Figure 2.3: An example to show how the number of transmitted photon pairs depends on p(k, 1).

Bkj and Nke are Poisson with means CkjAp(k,) and A,(k,I), respectively:


Nk, ~ Poisson(Ap(k,l)), k = 1, 2,... I, 1 E fank, (2.1)

Bk1 ~ Poisson(EklAp(k,,)), k = 1, 2, . . , I, I C fank, (2.2)


where Ap(k,l) depends only on p(k, 1). We refer to the data model described by equations (2.1) and (2.2) as Poisson Model I. To see the suitability of Poisson Model I, we compare real blank scan data and synthetic data generated according to Poisson Model I in Figs. 2.4 and 2.5. It can be seen from the figures that the synthetic data bears a close resemblance to the real data.

Note, a rod passes through the same detector tube twice in each cycle. This fact is implicitly accounted for through the parameters {Ap(k,1)}.

From Poisson Model I, it follows that


E[BkL] = k1E[Nki].


(2.3)










20


600 550 500


450 400 350 300


250 200


150 100 so-


20 40 60 s0 100 120 140 160



Figure 2.4: Overlapped 192 projections of real blank scan data.


550500


400 360 300


250 200 10 100


20 40 60 80 100 120 140 10



Figure 2.5: Overlapped 192 projections of synthetic blank scan data generated according to Poisson Model I.





Hoffman's assumption in (1.4) follows by replacing the expectations E[Bki] and E[Nk] in (2.3) with their ML estimates bkl and n1k, respectively. The statement in (2.3) is more sound probabilistically than the assumption in (1.4) because detector efficiencies should be interpreted as probabilities. Moreover, all we can say about BkL and Nkz is that their means obey the relationship in (2.3).






21


2.1.1 Estimation Problem Associated with Poisson Model I

For convenience, let B, N, E, and A be vectors with elements {Bk1}, {NOl}, {E}, and {AP(k,l)}, respectively. Also, let b be a vector with elements {bki}, where bkl is an observation of Bkl. To estimate E and A, we use the ML method:


(P1) (i,i) = arg max 01(E, A) subject to 0 < E < 1 and A > 0, (2.4)


where 01(E, A) is the log-likelihood function

I
1 (E, A) A log P(B = b; E, A) = log ] e-kAp(kl) (EklA(k,l))2.5) k=1 Efank bk!

and 0 and 1 are I x 1 vectors of zeros and ones, respectively. In Appendix A, we show that Hessian of the log-likelihood function is not negative definite. Thus, the loglikelihood function is not concave. Given this fact, it is expected that the performance of iterative algorithms used to solve (P1) will depend on the initialization procedure.


Cramer-Rao Lower Bound

To compute the Cramer-Rao lower bound (CRB) for the estimators of E and A in

(P1), we need to first compute the Fisher information matrix. We prove in Appendix B that the Fisher information matrix is non-invertible for the probability density function associated with (P1). Therefore, the CRB for the estimators of e and A cannot be determined.


2.1.2 Estimation Algorithms associated with Poisson Model I

If n, an observation of N, was available, then the parameters E and A could be estimated using

E b, k = 1, 2,...,I, E fank, (2.6)
nak






22


E' nki
Ad '(k,1):p(k,1)=d} d = di, d d (27)
number of pairs (k, 1) such that p(k, 1) = d ' 2, - - , dD,

where {dj} are the possible perpendicular distances and D is the number of possible perpendicular distances. Given this observation, our choice for the complete data is X = [B, N]. It should be noted that an observation of N is not available. Consequently, as part of the EM algorithm, an estimate of an observation of N will be generated.

The E step of the EM algorithm requires the computation of the conditional expectation of the log-likelihood function of the complete data given the incomplete data and current estimates of E and A. Let U(e, A; E', A") denote this conditional expectation, where E' and A' are the estimates of the detector efficiencies and mean number of transmitted photon pairs, respectively, at the mt" iteration. Also, let Lc(b, n; E, A) denote the log-likelihood function of the complete data:


L,(b, n; E, A) = log P(B = b, N = n; E, A). (2.8)


In Appendix C, we show that


U(E, A; E", Am) A E[LC(B, N; E, A) I B = b; E', Am]

= Y Z [bkl log(E6) + (n' - bkl) log(1 - EWE)] k=1 iEfan k
dD
+ [-Ad + n' log Ad] + C1, (2.9)
d=di {(k,1):p(k,1)=d}

where C, is a constant that includes all of the terms that are independent of E and A, and n' is the conditional expectation of NkI given the incomplete data and the estimates E' and Am:


nkl = E[Nk, I Bk = bkl; E", A"], k = 1 E fank (2.10)






23


We also show in Appendix C that


n' = (1 - e'e') Ak,l) + bkl, k = 1, 2,. , I, I E fank. (2.11)


The quantity n' is an estimate of the number of photon pairs transmitted along detector pair (k, l). From (2.11), it can be seen that n' is estimated by adding an estimate of the number of photon pairs that are transmitted but not detected, (1 - eg )A,),1), to the number of photon pairs that are transmitted and detected, bkl.

We recognize from (2.9) that U(E, A; E', A') is decoupled in the sense that it can be written as the sum of functions fi and g that depend only on A and E, respectively. More specifically,


U(E, A;E",A") = f"(A) g(") (E) + C1, (2.12)


where
dD
f"m)(A) A Z Z [-Ad + n- log Ad], (2.13)
d=di {(k,1):p(k,l)=d}
and

g(")(E) Z [bka log(ekt) + (n' - bk) log(1 - EkE)]. (2.14)
k=1 lEfan k
From (2.12), it follows that the problems of maximizing U(E, A; EmA') with respect to E and A are decoupled.

Summarizing our results, the steps of the EM algorithm are:


Step I. (Initialization)

Set m=O.

Choose c0 and A0 s.t. 0 < E0 < 1 and A > 0.

Step II. (E step)






24


k1 = (I - e'el) + kb1, k = ,2 . ,lE fankStep III. (M step)

(PA) A+ arg max ff)(A).
A\>o
(Pe em1 Aarg m1ax g ().
Step IV.

If stopping criterion is not satisfied, set m = m + 1 and go to Step I.

Otherwise stop.

Solving (PA), which is equivalent to maximizing U(E, A; c'A') with respect to A, is straightforward and leads to a closed-form solution. Momentarily ignoring the constraints

Ad>O, d=di,d2,...,dD, (2.15)

we take the partial derivative of flm)(A) with respect to each Ad and set it to zero:


affm)(A) _1
= { : (-1 +n m -)= 0, d d d2,...,dD- (2.16)
0 d {(k,1):p(k,1)=d} kAd

Solving (2.16) yields the desired update

Z
Snkl
Am+1 {(k,1):p(k,1)=d}d=did d (.)
number of pairs (k, 1) such that p(k, 1) = d ' 2, ... , dD (2.17)


From (2.11), it can be seen that the constraints in (2.15) are satisfied if cm' E[O, 1] and A' > 0 for all m, k and d.

Unfortunately, a closed-form solution is not available for the maximizer of the function U(E, A; EmA'm) with respect to E (see (PE)). We iteratively solve this maximization problem using either the coordinate descent (CD) method [32, pages 283-287] or a fixed-point iteration. See Abatzoglou et al. [33] or Luo et al. [34] for information on the convergence properties of CD method.






25


EM-CD Algorithm

To solve (PE), we first use the CD method where the function g(m) is minimized iteratively on a coordinate by coordinate basis. More specifically, the desired minimizer, Em+l, is given by solving the following maximization problems iteratively for k = 1, 2,. . ., I:


(k1i+ arg mnax g C),m1,i 2g1,+ k+,-1 k+1 ,1i) (2.18)
O<(<1


where 6g+1,i+1 is the estimate of Em+1 at the (i + i)h iteration. A reasonable choice for the initial estimate of Em+1 is Em+l,O = Em. The iteration is terminated when successive estimates are virtually the same. In practice, the algorithm is terminated when
I m+1,i+1 m+1,i
( -E _)2 < oe, (2.19)
k=1 ek
where 6, is a small positive number. The iterate that results when the stopping criterion is satisfied is considered the solution to (PE).

To see how the maximization problems in (2.18) can be solved, let


01(E) = g(m)(zi, z, .. , zI1, z1+1, . . , zI), 1 = 1,2,..., I, (2.20)


where z is a vector whose elements, z1, z2,..., zj, all lie in the interval [0, 1]. It is straightforward to see that


_2_____ -bkl Z2
[S - (-( - bk) 1 2] (2.21)
092 kEfan, (1 Zk)


is less than zero, for 1 = 1, 2, ... , I, if Em E[0, 1] and A' > 0 for all m, k and d. Hence, 41 is a strictly concave function if E' E[0, 1] and A' > 0 for all m, k and d. Because of this concavity property, the problems in (2.18) can be readily solved with arbitrary






26


accuracy using 1-D optimization algorithms (e.g., dichotomous search method, golden section method, and Fibonacci search [32, pages 268-275]). The constraint that every Ek lies in the interval [0, 1] can be easily enforced by using the interval [0, 1] as the initial interval of uncertainty for the chosen 1-D minimization method. Hence, the algorithm is guaranteed to provide detector efficiency estimates that lie in the interval [0, 1] when co E[0, 1] and A' > 0 for all k and d. We refer to the EM algorithm that results from using the CD method to solve (PE) as the EM-CD algorithm. To solve the maximization problem in (2.18), we use the dichotomous search method [32].


EM-FP Algorithm

The second way we solve the problem (PE) is to use a fixed-point iteration. Although it guarantees nonnegative detector efficiency estimates, the fixed-point iteration is not theoretically guaranteed to have the property that the estimates all lie in the interval [0, 1]. Nevertheless, in all of our simulations, the detector efficiency estimates obtained using the fixed-point iteration did have this property. An advantage of the fixed-point iteration is that it converges faster in practice than the coordinate descent algorithm.

Ignoring the constraint in (PE), our approach is to use a fixed-point iteration to solve the system of equations


(')= 0, k = 1, 2,. . .., I. (2.22)
O9fk

Taking the partial derivative of g(") with respect to Ek gives


,9() b+ (n' - bkl) " 1. (2.23)
196k 4Ej k1 EkEl
lEfank






27


Thus, the resulting system of equations are


[-(n'kj-bk) ] = 0, k=1,2,...,I. (2.24)
lEfank Ek 1

Obviously, a closed form solution for Em+l is not possible. To find an approximate solution, we write (2.24) as


Sbkl(- + ) = , k = 1, 2, ..., I. (2.25)
lefank Ck 1 - fk(I Efank EkQ Simplifying the left-hand side of (2.25) yields


1k bk l -_ _ _ _ , k = 1 ,2 ....I1. (2 .26 )
1cfan, lEfank

The following fixed-point iteration is then used to find an approximate solution to the system of equations in (2.26):


Mr+l,i+i - 1cankb k I
+1,i+ le an , k = ,2, .. ., I. (2.27)

lEfank 1 I

As with the EM-CD algorithm, the initial estimate of Em+1 is Em+1,O = Em. Further, the solution to (PE) is the iterate obtained when the stopping criterion in (2.19) is satisfied.

If nm in (2.11) is positive for all (k, 1), then every term in (2.27) is nonnegative. Thus, the detector efficiency estimates will be nonnegative provided 0 < Em < 1 and A' > 0 for all m. We refer to the EM algorithm that results from using the discussed fixed-point iteration as the EM-FP algorithm.






28


Over-Parameterized Approach

We now assume the unknowns are the set of parameters {6ki}. This case is referred to as over-parameterized case because the number of unknowns in this case exceeds the number of data points.

Let e be the vector with elements {Ekl}. The conditional expectation of the loglikelihood function of the complete data given the incomplete data and the current estimates of e and A is


U,(e, A; e', Am) = > [bk log(Eki) + (n' - bkl) log(1 - EkI)
k=1 g6 fan k
dD
+ [Ad + n' log Ad] + Co, (2.28)
d=di {(k,i):p(k,1)=d}

where C, includes all the terms that are independent of e and A. The function U(e, A; e', A"m) is decoupled:


U.(e, A; e', Am) = g(m") (e) + fl") (A) + C', (2.29)


where

g")(e) = j [b log(ek) + (nm - bkl) log(1 -- El)] (2.30)
k=1 iefan k

and
dD
() [-Ad + nm log Ad]. (2.31)
d=dl {(k,1):p(k,1)=d}

Because the function U0(e, A; em, Am) is decoupled to two functions g(m)(e) and f(')(A), the problem of maximizing U0(e, A; em, Atm) with respect to e and A is decoupled.

The E step for the over-parameterized formulation is the same as the E step for the sufficiently-parameterized formulation except that EkE, is replaced by Ekl. In other words, for k = 1, 2,..., I, 1 Efank, the conditional expectation of N, given the






29


incomplete data and the estimates em and A' is


A
nki = E[Nkl | Bk1 = bkl; em, Em] = (1 - E4)A'k,) + bkl. (2.32)


However, the M step differs. New estimates for e and A are obtained by maximizing U0(e, A; em, A) with respect to e and A under the constraints


0 < Ckl <; 1, k = 1, 2, . .., I1, 1 E fank, (2.33)


and

Ad > 0, d = di,d2 ..., dD. (2.34)

Ignoring the constraints (2.33) and (2.34) momentarily, we take the derivatives of g(m) (e) and f(m) (A) with respect to every Eki and Ad, respectively, and set them to zero.


ag(m)(e) 1 _
= bkl - + (nm - bkl) = 0, k = 1, 2,. .., I, l E fank. (2.35)
aeki Eki ki -Ekl

=f~m)) - Z (-1+n -)=0, d=di,d2,...,dD. (2.36)
OAd {(k,1):p(k,1)=d} d


Solving (2.35) and (2.36), we have the following iteration.


m+1 -bkl
fk1 =- , k =1,2, ..., I, I E fank. (2.37)

nknl


number of pairs (k, 1) such that p(k, 1) = d , d = dD, d2, . . , dD- (2-38) For the constraints (2.33) and (2.34) to be satisfied, E01 must lie in the interval [0, 1] for k = 1, 2, ... , I, I E fank and A' must be nonnegative for d = di, d2,..., dD. As mentioned earlier, the results using this approach are poor. Consequently, we do not include this approach in our simulation studies.






30

2.2 Poisson Model I

As we mentioned earlier, the CRB for the estimators of E and A in (P1) cannot be determined because the Fisher information matrix for the PDF in (P1) is noninvertible. One of the possible approaches to changing the non-invertibility is to reduce the number of unknowns in the estimation problem [35]. In this section, we present a second Poisson model that incorporates a priori information about the blank scan data so that the number of unknowns in the estimation problem is reduced.


k


, r os8,rgsin0)

0
02
01










Figure 2.6: The exact way to calculate Ap(kl), given A(rrod cos 0, rrod sin 0).


For detector pair (k, 1) (see Fig. 2.6), the mean number of photon pairs that are transmitted along the detector pair, Ap(k,l), can be expressed as [37]


Ap(k,1) J A(x, y)P(k, 11 x, y)dxdy, (2.39)
rod path

where A(x, y) is the number of photon pairs emitted from the rod at point (x, y) and P(k, 1 1 x, y) is the probability that a photon pair emitted from the rod at point






31


(x, y) is incident on detector pair (k, 1). Because the rod path is circular, it is more convenient to use polar coordinate than Cartesian coordinate. From the scanner geometry, Ap(k,l) can be expressed as

/02
Ap(k,) = 2 ]1 A(rrod cos 9, rrod sin 9)P(k, lIrrod cos 9, rrod sin O)rroddO, (2.40)


where A(rrod cos 9, rrod sin 9) is the mean number of photon pairs that are emitted from the rod at the angle 9, rrod is the radius of the path of the rotating rod, and P(k, lrrod cos 9, rrod sin 9) is the probability that a photon pair emitted from the rod at the angle 9 is incident upon the detector pair (k, 1). The factor of two accounts for the fact that the rod passes through the same detector tube twice in one cycle. We assume that the rod rotates at a constant angular velocity and emits positron uniformly. Under these assumptions, it follows that A(rrod cos 9, rrod sin 9) is constant over 9. Let A A 2A(rrod cos 9, rrod sin O)rrod. Then, the integral in (2.40) may be approximated by the following sum (see Fig. 2.7):


Ap(k,I) E AP(nA9)A9 Ac(k, 1), (2.41)


where

c(k,l) PIA90, (2.42)
n

and P,, P(nA9) is the probability that a photon pair emitted from the rod at the nth position is incident upon the detector pair (k,l). This probability may be calculated using the angle-of-vies (AOV) method [14,36]. Specifically, the probability that a photon pair emitted from a point (x, y) is incident on detector pair (k, 1) is determined by the angles a,, a2, a3, and a4 (see Fig. 2.8). The AOV is the maximum






32


angle for which a photon pair would be incident on a detector pair and is given by


A
AOV = min{ai, a2, 7r - a3, 7r - a4}. (2.43)


Then, the probability that a photon pair emitted from the point (x, y) is incident on detector pair (k, 1) is equal to AOV/7r [36]. Note, the parameter c(k, 1) depends



k



P2










--- -- ---





Figure 2.7: A proposed way to calculate c(k, 1) to approximate Ap(k,l).


only on p(k, 1). The parameters {c(k, 1)} for the Siemens-CTI ECAT EXACT 921 scanner are plotted in Fig. 2.9 as a function of (1 - k).

Equipped with the above information, we assume that Nk, and Bkj are Poisson with means Ac(k, 1) and ekAc(k, 1), respectively:


NkI - Poisson(Ac(k, 1)), k = 1, 2,.. ., I, E fank, (2.44) BkI - Poisson(EklAc(k, 1)), k = 1, 2,..., I, I E fank. (2.45)






33


k









41









Figure 2.8: The four angles to determine the angle-of-view associated with the point (x, y) and detector pair (k, 1).

We refer to the data model described in equations (2.44) and (2.45) as Poisson Model II.


2.2.1 Estimation Problem Associated with Poisson Model II

As discussed in the previous section, the parameters {c(k, I) } are determined from the geometry of the scanner. The unknowns Adi, Ad2,. . . , Ad for Poisson Model I have been reduced to only one unknown, A, for Poisson Model II.

The estimation problem now becomes


(P2) (i, A) = arg max 02(C, A) subject to 0 < E < 1 and A > 0, (2.46)






34

c(kI)
0.08 k


0.0750.07


0.0650.06



120 130 140 150 160 170 180 190 I-k
Figure 2.9: The parameters {c(k, l)} for the Siemens-CTI ECAT EXACT 921 scanner calculated according to (2.42). where 02(E, A) is the log-likelihood function:


#2(E, A) = log P(B = b; E, A) = log 1i fi e-kL Ac(kI) (EkLAc(k, 2.b47 k=1 Efan bkl!

We show in Appendix A that the function #2(E, A) is not concave.

Unfortunately, although the number of unknowns is reduced in (P2), the Fisher information matrix for the probability density function associated with (P2) is still non-invertible (see Appendix B).


2.2.2 Estimation Algorithms Associated with Poisson Model II

We show in Appendix C that the conditional expectation of the log-likelihood function of the complete data X = [B, NJ given the incomplete data and the current






35


estimates E' and A' is


U'(E, A; E', A')


A


E[L'(B, N; E, A) I B = b; E', A']
I


= S S[bkI log(Ekfl) + (n' - bkI) log(1 - EkE)]
k=1 iEfan k

+ 5 5[-Ac(k, 1) + n- log Ac(k, 1)] + C2,
k=1 lfan k


where L'(B, N; E, A) is the log-likelihood function of the complete data:


L'(B, N; E, A) = P(B = b, N = n; E, A),


and C2 is independent of E and A. We also show in Appendix C that, for k 1, 2,..., I,l E fank,


nk = E[Nk, | BkI = bkI; Em, A'] =(1 - eE')A'c(k, 1) + bkl. (2.50)


The function U'(E, A; E', A') is given by


U'(e, A; E', A') = g(") (E) + f2-)(A) + C2,


I
9(- = (,5 [bkI log(ekEl) + (n' - bkI) log(1 - EkfI)]
k=1 lEfan k


f2") (A) A 5 [-Ac(k, l) + n' log Ac(k, 1)].
k=1 IEfan k


(2.48)


(2.49)


where


(2.51)


and


(2.52)


(2.53)






36


The steps of the EM algorithm are:


Step I. (Initialization)

Set m=O.

Choose c0 and A0 s.t. 0 < c0 < 1 and A0 > 0.

Step II. (E step)

n' = (1 - ee"')A'c(k, l) + bk, k 1,2,..., I, 1 E fan.

Step III. (M step)

(PA) A-+1 arg max f2")(A).
A>O
(PE) E--1 arg mnax g(")(,E).

Step IV.

If stopping criterion is not satisfied, set m = m + 1 and go to Step II.

Otherwise stop.


The problem (PE) is exactly the same optimization problem that arose when Poisson Model I was used. Therefore, the detector efficiencies can be estimated using (2.18) and (2.27). To find the solution to (PA), we ignore the constraint A > 0 momentarily. Taking the derivative of f ')(A) with respect to A and setting it to zero, we get

[-c(k, l) + ]= 0. (2.54)
k=1 LEfan k
Solving (2.54), the new estimate for A is



Am+' - k=1 lefan k (2.55)
Z c(k, l)
k=1 iEfan k






37

We refer to the estimation algorithm using (2.18) and (2.55) for Poisson Model II as the EM-CD2 algorithm. We refer to the estimation algorithm using (2.27) and (2.55) for Poisson Model II as the EM-FP2 algorithm.














CHAPTER 3
ML ESTIMATION ALGORITHMS USING A REFINED DATA MODEL


3.1 Poisson Model III

To get a more complete data model for the blank scan data, the effect of detector penetration must be incorporated. In Fig. 3.1, although the photons are incident


k k'


















Figure 3.1: Detector penetration.


upon detector pair (k', 1'), the photons may penetrate the detectors. Suppose the photon incident on detector k' penetrates detector k' and is detected by detector k, and the photon incident on detector 1' penetrates detector ' and is detected by detector 1. Then, a coincidence may be recorded erroneously by detector pair (k, 1). It should be noted that a coincidence may also be erroneously recorded by detector


38






39


pairs (k, ') or (k', 1) when detector penetration takes place, even though the photons are incident upon detector pair (k', '). Even if the photons stop within detectors k and 1, the photons might not be detected by detectors k and 1. For a photon to be detected by a detector, the photon must have enough energy after stopping to interact with the detector [25].

The probability that the photon penetrates detector k' (or survival rate) is e-Lk' [38, 39], where Lk' is the length of intersection of detector k' and the path of the photon (see Fig. 3.2) and p is the attenuation coefficient of the detector. Therefore,















Figure 3.2: Definitions of Lk and Lk'.


for the given path of the photon pair, the probability that the photon incident on detector k' is detected by detector k is equal to


Pke-Lk '( - e-Lky), (3.1)


where pk is the probability that a photon event is recorded given that a photon has been stopped by detector k, and (1 - e-LP) is the probability that the photon is stopped by detector k. The attenuation coefficient of detectors is assumed to be known. For BGO detectors, the linear attenuation coefficient is 0.96 cm-1 [25] and the intersection lengths can be calculated from the scanner geometry. Although we






40


focused on one photon of the photon pair, similar statements can be made about the other photon.

When the effect of detector penetration is considered, the tube associated with a detector pair is defined differently (see Fig. 3.3). The mean number of photon pairs
























Figure 3.3: The tube associated with detectors k and I is defined by the dashed lines when detector penetration is taken into account. Note that, to clearly depict the angles, detectors k and I are magnified.


that are detected by detector pair (k, 1) is given by


E[BkI]

= f A(x,y)PkpiP(k,lix,y)dxdy (3.2)
rod path

2 f4 A(rrod cos 0,rrod sin 0)pkpIP(k, 1Irrod cos 0, rrod sinO)rroddO,

where A(rrod cos 0, rrod sin 0) is the mean number of photon pairs that are emitted from the rod at angle 0, and P(k,1jrrodcos0,rrodsin0) is the probability that a






41


photon pair emitted from the rod at angle 0 is stopped within detector k and 1, respectively. Since many paths are possible for a photon pair that is incident on detector pair (k, 1), the probability P(k, lrrod cos 9, rrod sin 9) can be expressed as


P(k, lrrod cos 0, rrod sin )p)(
f o 3+AOV e-L,,(,)ju (1 - e-kI))e~l(Ta/~


where Lk (a) is the intersection length of the path of the photon traveling at angle a and detector k, and Lk'(a) is the intersection length of the path of the photon traveling at angle a and detector k' (see Figs. 3.4 and 3.5). The intersection lengths L,(7 + a) and L1, (r + a) are defined similarly. Note that a photon might have


k AOV












01 0

I I
Figure 3.4: Angles that are used to calculated P(k, lrrod cos 9, rrod sin 0) in (3.3).


penetrated several detectors before reaching detector k. In this case, (3.3) would have to be modified accordingly.






42


k AOV








I N





Figure 3.5: Angles that are used to calculated P(k, lrrod cos 0, rrod sin 9) in (3.3).


For convenience, we let A = 2A(rrod cos , Trod sin 0)rrod. In practice, the exact integral in (3.2) cannot be determined. Therefore, we approximate the integral by N9
APkP Z P(nAO)AO, (3.4)
n=1

where


P(nAO) = Zek(ta)(j - e-Lk (A))eLL,(1+iAQ)l e- , (3.5)


AO A (02 - 01)/N9, (3.6)

and
A
Aa = AOV/N. (3.7)

The quantity in (3.4) can be expressed as


APkpla(k, 1),


(3.8)






43


where


No
a(k, 1) E P(nAO)AO.
n=1


(3.9)


Observe that a(k, 1) depends only on p(k, 1). The parameters {a(k, l)} for the SiemensCTI ECAT EXACT 921 scanner are plotted in Fig. 3.6. The number of photon pairs


a(kI)


120 130 140 150
1-k


160 170 180 190


Figure 3.6: The parameters {a(k, l)} for the Siemens-CTI ECAT EXACT 921 scanner calculated according to (3.9).


that are detected by detector pair (k, 1) is therefore assumed to be an observation of the random variable Bk1, where


Bk1 - Poisson(APkpia(k, 1)), k = 1, 2, .. ., I, 1 E fank.


(3.10)


In order to utilize the same methodology for estimating unknown parameters of Poisson Model II, the means of random variables {Bkl} must be related to the means of some other Poisson random variables in the same manners that the means of {Bkl} are related to the means of {Nkl}. Let Nkl be the number of photon pairs that are


x 10-"


6





5


4


- -






44


stopped within detectors k and 1. Then, the random variables {Nkl} have means {Aa(k, l)}:

N kl ~ Poisson(Aa(k, 1)), k = 1, 2, ..., 1 E fank. (3.11) We refer to the data model described in equations (3.10) and (3.11) as Poisson Model III.


3.1.1 Estimation Problem Associated with Poisson Model III

Since the parameters {a(k, 1) } can be determined from the scanner geometry (see (3.5) and (3.9)), the estimation problem now becomes


(P3) (, A) = arg max 03(p, A) subject to 0 < p <; 1 and A > 0, (3.12)
P"\

where p is the vector with elements {Pk} and 03(p, A) is the log-likelihood function:


03(p, A) = log P(B = b; p, A) = log f1 11 e-kpIa,)(pkplk ()) k=1iEfan (

We show in Appendix A that the function 3(p, A) is not concave. In addition, we show in Appendix B that the Fisher information matrix for the probability density function associated with (P3) is non-invertible. Therefore, the Cramer-Rao lower bound for the estimators of p and A cannot be determined.


3.1.2 Estimation Algorithms Associated with Poisson Model III

Because Poisson Model III has the exact same functional form as Poisson Model II, the derivation of estimation algorithms parallels that for Poisson Model II.

We show in Appendix C that the conditional expectation of the log-likelihood function of the complete data X = [B, N] given the incomplete data and the current





45


estimates pm and A' is


U"(p, A; p', Am)


A


E[L' (B, N; p, A) I B = b; pm, Am)
I


[bkI log(pkpI) + (n' - bkl) log(1 - pkpI)]
k=1 icfan k

+ Z ~ [-Aa(k,1)+n-logAa(k,1)]+C3,
k=1 IEfan k


(3.14)


where L'(B, N; p, A) is the log-likelihood function of the complete data:


L''(B, N; p, A) = log P(B = b, N = n; p, A),


(3.15)


and where C3 is independent of p and A, We also show in Appendix C that, for k = 1, 2,..., I, l E fank,


nk = E[Nk| BkI = bkl; p', Am] = (1 - p'm, )Ama(k, 1) + bkl, The function U"(p, A; pm, At) is decoupled and given by


U"(p, A; p', Am) = g(m) (p) + f") (A) + C2,


(3.16)


(3.17)


(3.18)


dD f(") (A) A E
d=dl


(3.19)


) [-Aa(k, 1) + n' log Aa(k, 1)].
{ (k,l):p(k,l) =d}


where


and


g(m)(p) [ki log(pkpI) + (n' - bkI) log(1 - PkP)]
k=1 tEfan k






46


The steps of the EM algorithm are:


Step I. (Initialization)

Set m=.

Choose p0 and AO s.t. 0 < p0 < 1 and A0 > 0.

Step II. (E step)

nn = (1 - pip")Ama(k, l) + bk, k = 1, 2,..., I, 1 E fankStep III. (M step)

(PA) Am+' arg max f'm)(A).
A>0
(Pp) pm+l arg max g(m)(p).

Step IV.

If stopping criterion is not satisfied, set m = m + 1 and go to Step II.

Otherwise stop.

Repeating the steps in Section 2.2.2, the problem (Pp) is solved either by using a CD-based method


+,i+= arg a )m+1,i+1 m+1,i+1 m+1,i+1 m+1,i m+1,i

(3.20)
or a fixed-point iteration


+1,i+1 efan1pm+1, Pk E mmZi(.1

lEfan, I p where pj"'+Z+ is the estimate of pm+1 at the (i + 1)1h iteration for k 1, 2,..., I. The new estimate for A is given as



A"'1- k=1lEfank (3.22)
E a(k,1)
k=1 lefan k






47

We refer to the estimation algorithm using (3.20) and (3.22) for Poisson Model III as the EM-CD3 algorithm. We refer to the estimation algorithm using (3.21) and (3.22) for Poisson Model III as the EM-FP3 algorithm.














CHAPTER 4
SIMULATION


In this chapter, we assess the performance of the proposed algorithms using synthetic PET data and data from a thorax phantom by Data Spectrum Inc. For the purpose of making comparisons, we apply the fan-sum algorithm and Ferreira's algorithm to this data as well. Note although Ferreira's algorithm was developed for a cylindrical source, the underlying assumptions could also be made for rotating rod sources.

For all the iterative algorithms, the initial estimates for the detector efficiencies are d = 0.5 for all k. Given this choice, the initial estimate for A for the EM-CD and EM-FP algorithms is



Axo = {(k,1):p(k,1)=d} I d = di d d (4.1)
number of pairs (k, 1) s.t. p(k, 1) = d ' ' 2,... , d,

the initial estimate for A for the EM-CD2 and EM-FP2 algorithms is

b
0 _ k=1 gEfank k
- I (4.2) E E c(k, 1)
k=1 lEfan,

and the initial estimate for A for the EM-CD3 and EM-FP3 algorithms is

b

A0 = k=1 Efank a (k, l)
k=1 lEfank


48






49


The steps of the EM-CD and EM-FP algorithms are repeated until the following stopping criterion is met.

I m+1 m dD Am+' m
(ek 'Ek ) 2 Z d d )2<6(4
k=1 Ek d=d1 A

where 6 is a small positive number. In our simulations, the EM-FP algorithm converged much faster than the EM-CD algorithm. Consequently, we used 6, = 6 = 10-1 for the EM-FP algorithm and 6, = 6 = 10-' for the EM-CD algorithm. As for both the EM-CD2 and EM-FP2 algorithms, the stopping criterion is I m+1_ Am+'_ -m
k k 2 2 < 6, (4.5)
k=1 Ek Am

where 6 is a small positive number. In our experiments, 6, = 6 = 10- for the EM-FP2 algorithm and 6, = 6 = 10-2.5 for the EM-CD2 algorithm. By contrast, Ferreira's algorithm did not converge even after 300,000 iterations when the value 6, = 10- was used. Therefore, in our simulations Ferreira's algorithm was executed for 250 iterations as specified by Ferreira et al. [17].


4.1 Synthetic Data


Using Poisson Model II, synthetic data are created to simulate blank scan data obtained from the Siemens-CTI ECAT EXACT 921 scanner. The scanner has 384 detectors per ring, and the ring diameter is 84.5 cm. The region to be reconstructed is a 43.9x43.9 cm2 square consisting of (128)2 voxels. Each plane consists of 192 projections and 160 members per projection. More information of this scanner can be found in Wienhard et al. [26]. As discussed in Chapter 2, this scanner has 81 possible perpendicular distances. The synthetic data consist of observations of Poisson random






50


variables {BkI}:


Bk ~Poisson(EkqEjAc(k, 1)), k = 1, 2,. . ., 384,1 E fank, (4.6)


where

A = 69000, (4.7)

and {c(k, I)} as discussed in Chapter 2. The value of A in (4.7) is chosen so that the synthetic blank scan data corresponds to 2-hour scan. Uniform, piece-wise uniform, and nonuniform detector efficiencies are used in our experiments presented.

As customary [17], the detector efficiency estimates are normalized so that for each estimation algorithm the mean of the detector efficiency estimates is equal to one.

Because the EM-CD3 and EM-FP3 algorithms are based on Poisson Model III, which incorporates the effect of detector penetration, the EM-CD3 and EM-FP3 algorithms are not included. The reason is because synthetic data do not account for the effect of detector penetration. On the other hand, if the synthetic data did incorporate the effect of detector penetration, it would be unfair to compare the EMCD3 and EM-FP3 algorithms against estimation algorithms that do not take detector penetration into consideration.


4.1.1 Comparison of Detector Efficiency Estimates

For our first experiment, the detector efficiencies were uniform:


EkO= .8, k =1, 2,..., 384. (4.8)


The true efficiencies are shown in Fig. 4.1(a), and the estimates provided by the fan-sum, Ferreira's, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are plotted






51


in Figs. 4.1(b)-4.1(g), respectively. From the estimation results, we conclude for all four of the proposed algorithms that the detector efficiency estimates have smaller variance. Note that the fan-sum estimates vary over a relatively wider range and the estimates provided by Ferreira's algorithm show a pronounced sinusoidal pattern.

For our second experiment, the detector efficiencies were piece-wise uniform:



Ek 0.,k=1 ,.. 9 49
0.4, k = 193,194, 392

The true efficiencies are shown in Fig. 4.2(a), and the estimates provided by the fan-sum, Ferreira's, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are plotted in Figs. 4.2(b)-4.2(g), respectively. The fan-sum algorithm provides very poor estimates, while the estimates provided by Ferreira's algorithm again exhibits a sinusoidal pattern. From the experiments, the following conclusions are drawn about the proposed algorithms. The efficiency estimates for the first half of the detectors show little deviation from a constant value. This is also true for the second half of detectors. In addition, the ratio of the efficiency estimates for the first half of detectors to the efficiency estimates for the second half of detectors is very close to the true ratio

0.8/0.4.

To assess performance of estimation algorithms, we compare the mean-square error (MSE) of the efficiency estimates of 50 Monte Carlo (MC) runs. The error criterion MSE is defined as
I 1 50
MSE = - Z( k(i) - Ek), (4.10)
k=1 i=1

where ek(i) is the estimate for Ek in the i'6 MC run. The results are tabulated in Tables 4.1 and 4.2, respectively. Although each of four of the proposed algorithms provides detector efficiency estimates with smaller MSE than the fan-sum and Fer-






52


reira's algorithms, the difference between all algorithms is small. Further comparison is made based on a different criterion in the following. Table 4.1: MSE of the detector efficiency estimates using synthetic data and uniform detector efficiencies in the first experiment.


Algorithm Fan-sum Ferreira's EM-CD EM-FP EM-CD2 EM-FP2
MSE 15.3721 15.4538 15.3623 15.3623 15.3630 15.3630




Table 4.2: MSE of the detector efficiency estimates using synthetic data and piecewise uniform detector efficiencies in the second experiment.


Algorithm Fan-sum Ferreira's EM-CD EM-FP EM-CD2 EM-FP2
MSE 69.3189 69.5227 68.0414 68.2702 68.2334 68.1167



Ideally, the detector efficiency estimates would be such that the ratio Ek/Ek is approximately equal to a constant for all k. Therefore, as a merit of figure, we use the variance of the ratio of the estimates and true values:


VR = sample variance of {/Ei, 2/E2, . . . , e/f}. (4.11)


The idea is that the smaller the VR the better. To further compare the algorithms, 50 MC runs were performed for each of the first and second experiments and the VR was computed for each run. The plots of the VR versus MC runs are shown in Figs. 4.3 and 4.4 for the first and second experiments, respectively. From the plots, it can be clearly seen that all four of the proposed algorithms have comparable performance and each of them have much smaller VR than the fan-sum algorithm and Ferreira's algorithm.






53


For our third experiment, we used nonuniform detector efficiencies that were generated as follows:


Ck{ if 0 1 k = 1)2,..., 384) (4.12)
1 if (a> 1

where

= 0.5 + v/0.008A(0, 1), k = 1, 2,...,384. (4.13)

Obviously, only positive values are used for {M}. The parameters in (4.13) are chosen so that with high probability the detector efficiencies lie in the interval [0.3, 0.7], which is a typical range for detector efficiencies.

The ratio of detector efficiency estimates to true detector efficiencies is plotted in Figs. 4.5(a)-4.5(f), where the detector efficiency estimates are provided by the fansum, Ferreira's, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms, respectively. Based on Fig. 4.5, we confer that all four of the proposed algorithms provide much better detector efficiency estimates than the fan-sum and Ferreira's algorithms. Tabulated in Table 4.3 is the MSE of the detector efficiency estimates generated by using each of the estimation algorithms. All four of the proposed algorithms provide detecTable 4.3: MSE of the detector efficiency estimates using synthetic data and nonuniform detector efficiencies in the third experiment.


Algorithm Fan-sum Ferreira's EM-CD EM-FP EM-CD2 EM-FP2
MSE 100.8640 100.7207 100.2472 100.2704 100.2215 100.2070



tor efficiency estimates with smaller MSE than the fan-sum and Ferreira's algorithms. We also compute the VR in each of 50 MC runs. The results are plotted in Fig. 4.6. It can be seen that all four of the proposed algorithms have comparable performance, and each of them is clearly better than the fan-sum algorithm and Ferreira's algo-






54


rithm in the sense that they provide estimates with smaller VR in every MC run than the fan-sum and Ferreira's algorithms.

Another comparison criterion adopted by other researchers is the relative error [9,17]:
- Ek L ~ ,.. (4.14)
Ek

The maximum relative errors are plotted in Fig. 4.7, and the average relative errors are plotted separately in Fig. 4.8. Based on the maximum relative error, all four of the proposed algorithms have comparable performance and each of them is better than the fan-sum and Ferreira's algorithms. In terms of mean relative error, the EM-FP algorithm has the best performance, and both the EM-CD and EM-CD2 algorithms are better than the fan-sum and Ferreira's algorithms in every MC run. To remove the bias factor from the comparisons, we normalize the detector efficiency estimates so that their mean is equal to the mean of the true detector efficiencies and then compute the relative errors. Note that in practice this cannot be done because the true detector efficiencies are not known. The maximum and average relative errors associated with the normalized detector efficiency estimates are plotted in Figs. 4.9 and 4.10. The plots confirm that all four of the proposed algorithms have comparable performance, and each of them outperform the fan-sum and Ferreira's algorithms. Specifically, the proposed algorithms provide estimates with smaller maximum and average absolute errors than the fan-sum algorithm and Ferreira's algorithm in every MC run.


4.1.2 Comparison of Reconstructed Emission Images

The efficiency estimates were then used to reconstruct a synthetic emission image. Assuming only errors due to detector inefficiency, the emission data are assumed to be a realization of the set of Poisson random variables {Yk,} given in (1.3). For






55


convenience, the data model for {Ykj} is expressed again. The emission data, {YkI}, are modeled as observations of the independent random variables {Ykl}:


Ykj~Poisson(ekel[PX]sk), k = 1, 2,...,384,1 E fank,


where
[Px]kiA Pii (4.15)
j=1
and xj is the mean number of photon pairs emitted from voxel j and J is the number of voxels. In addition, Pklj is the probability that a photon pair emitted from voxel j is incident on detector pair (k, 1) and x = [x1, x2, ... xJ]' is referred to as the emission intensity vector. Given the observed emission data y, which is a vector with elements {yA}, the problem is to estimate the emission intensity vector x. Given the synthetic emission intensity x in Fig. 4.11, we use the EMML algorithm [1, 14,40,41] to estimate x from y:
, n~j) (n+1)
S c EL [pci(f)~ (4.16)
k,1
where

Pkck,j = EklEPkIj. (4.17)

The efficiency estimates were used in (4.17) to reconstruct the emission image.

In Figs. 4.12(a)-4.12(d) are the plots of the reconstructed emission images at the 15th iteration of the EMML algorithm using true detector efficiencies, no normalization, fan-sum algorithm, and Ferreira's algorithm, respectively. (The EMML algorithm virtually converges after 15 iterations.) And in Figs. 4.13(a)-4.13(d) are the plots of the reconstructed emission images at the 15th iteration of the EMML algorithm using the EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms, respectively. Because all the reconstructed images look very similar to each other, 50 MC runs were used to compare the MSE of the reconstructed images. The criterion MSE is






56


defined for each iteration of the EMML algorithm as


MSE = 1 (i) - Xj)2, (4.18)


where Ji(i) is the estimate for xj at the ith MC run. The curves of the MSE of reconstructed images versus iteration number of the EMML algorithm are shown in Fig. 4.14. Because the mean of true detector efficiencies is different from the mean of estimates of each estimation algorithm, bias becomes an issue and the MSE corresponding to true detector efficiencies is not included in the plot. While all

Table 4.4: MSE of the reconstructed emission images from synthetic data.


Iteration of EMML 16 17 18 19 20
Fan-sum 786.1761 785.5547 788.4473 794.3706 802.9343
Ferreira's 786.1472 785.5140 788.3954 794.3079 802.8612
EM-CD 786.1328 785.4880 788.3587 794.2610 802.8044
EM-FP 786.1363 785.4823 788.3440 794.2375 802.7723
EM-CD2 786.2022 785.3626 788.0479 793.7720 802.1426
EM-FP2 786.1835 785.2883 787.9220 793.5967 801.9190


the estimation algorithms improve the reconstructed emission images over using no normalization, the curves for all the estimation algorithms are close together. The values of MSE are tabulated in Table 4.4 for iterations 16-20 of the EMML algorithm for clarity.


4.2 Real Data


Real PET data were obtained from Dr. John Votaw of the PET Laboratory at the Emory University School of Medicine in Atlanta, GA. The PET scanner is the Siemens-CTI ECAT EXACT 921 model. The collected data consist of 47 planes, where each plane consists of 192 projections and 160 members per projection. The






57


scanner has 384 detectors per ring and the region to be reconstructed is a 43.9x43.9 cm2 square consisting of (128)2 voxels.

The fan-sum method, Ferreira's method, and the proposed algorithms were used to estimate the detector efficiencies for plane 10. The duration of the blank scan was 30 minutes. Since the true detector efficiencies are unknown, we have no objective way to compare the detector efficiency estimates. The efficiency estimates were used in conjunction with the EMML algorithm to reconstruct the emission image for plane 10. The duration of the emission scan was 7 minutes.

In Figs. 4.15(a)-(d) and 4.16 (a)-(b), the reconstructed emission images at the 20th iteration of the EMML algorithm using the fan-sum, Ferreira's, EM-CD, EM-FP, EM-CD2, EM-FP2 algorithms, respectively, are shown. The reconstructed emission images are indistinguishable. As the EM-CD3 and EM-FP3 algorithms are based on Poisson Model III which incorporates the effect of detector penetration, when applying the EM-CD3 and EM-FP3 algorithms to reconstruct emission images, the probability matrix P in (1.3) must be defined differently. Let Pkl,3 stand for the probability that a photon pair is stopped by detectors k and 1 given that the photon pair is emitted from voxel j. The methods by Shepp et al. [14] and Carroll [36] can be extended to compute these probabilities [42]. The reconstructed emission images using the EM-CD3, EM-FP3 algorithms and no normalization (i.e., {k = 1}) are plotted in Figs 4.17(a)-(c). The reconstructed emission images for the EM-CD3 and EM-FP3 algorithms show details (e.g. heart, lungs, and spine) much more clearly than the reconstructed emission images for the other estimation algorithms. Nevertheless, it is difficult to objectively determine whether using the EM-CD3 and EM-FP3 algorithms lead to better reconstructed emission images than using no normalization.

To date, there is no accepted methodology for assessing the performance of detector efficiency estimation algorithms when real data are used. Dr. Votaw suggested the following experiment. A cylindrical phantom with uniform activity could be scanned






58


for an extremely long time ( > 48 hours). Then, the emission images corresponding to each detector efficiency estimation algorithm would be reconstructed. Finally, the average variance of the pixel values for each image would be computed. The detector efficiency estimation algorithm that gave rise to the smallest variance would be considered to be the best algorithm.

The shortcoming of the experiment described above is that only a single realization of the blank scan data is used. To account for the statistical properties of the detector efficiency estimation algorithms, multiple realizations of the blank scan data would be needed. Unfortunately, it is difficult to obtain multiple realizations because of the extremely high scanning costs.























(a)


100 200 300


(e)














100 200 300


1.04


(b)


1.02 [


0.98


0.961


1.04


1.02


(c)


1.04


(d)


1.021


1


0.98


S0.96' - 0.96'
100 200 300 100 200 300 100 200 300


()


0.98 [


0.96


100 200 300


1.04


(g)


100 200 300


Figure 4.1: True efficiencies of synthetic data in the first experiment and the estimates:
(a) true efficiencies; the others are the detector efficiency estimates generated by using
(b) fan-sum algorithm; (c) Ferreira's algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.


59


1.04 1.02


1


0.98


1.5


0.5



0


1.04 1.02


1


0.98 0.96


1.02 [


1


0.98[


0.96'


1


1


1


O v


J."19AL.L11LIJAh
7FTTTM


00%40


A*"























(a)














100 200 300


(e)














100 200 300


(b)
1.6


1.4 1.2

1

0.8 0.6

0.4


1.6

1.4 1.2


1


0.8 0.6

0.4


, , 100 200 300


(f)














100 200 300


1.6 1.4 1.2


0.8 0.6 0.4


1.6

1.4 1.2


1[


0.8 0.6

0.4


(c)














100 200 300


(g)














100 200 300


1.6 1.4 1.2


1


0.8 0.6 0.4


(d)














100 200 300


Figure 4.2: True efficiencies of synthetic data in the second experiment and the estimates: (a) true efficiencies; the others are the detector efficiency estimates generated by using (b) fan-sum algorithm; (c) Ferreira's algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.


60


1


0.8 0.6


0.4 0.21


1


1.6

1.4 1.2-


0.8

0.6

0.4-


1






61


X 10-4


-- - (a) Fan-sum algorithm x (b) Ferreira's algorithm
(c) EM-CD
0 (d) EM-FP
o (e) EM-CD2
(f) EM-FP2

x x X x
xx x XX X x x x XXXXX X x x
x x x) X X


A s~ .~ A ~q A A A A - A A A A A A A A - A - A A A - A A A


5 10 15 20 25
MC run numbe


30 35 40 45 50


Figure 4.3: VR of the detector efficiency estimates in the first experiment versus MC run number using synthetic data, uniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.


5


2



1


5RI19S S95 ,S9e@009M


4b
*x


"






































-4 k


C XXX XX X XX X XX X X X X XXX XXXX XXXX XXX X XXXXX XXX X X X X X X XX X$


)0 00 00000000OOOO0O000 OO0 0 O0000000000000OO0


5 10 15 20 25
MC run number


30 35 40 45 50


Figure 4.4: VR of the detector efficiency estimates in the second experiment versus MC run number using synthetic data, piece-wise uniform detector efficiencies, and:
(a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.


62


2



0



-2


- - - (a) Fan-sum algorithm
x (b) Ferreira's algorithm
(c) EM-CD
o (d) EM-FP
o (e) EM-CD2
A (f) EM-FP2


0


-6



-8



-10


'IV























(a)














100 200 300


(d)






- 2 100 200 300


(b)
2.1 '


2.05


2


1.95 1.9


2.1 2.05


2


1.95 1.9


100 200 300


(e)














100 200 300


2.1


2.05


2


(c)


1.95[


1.9


2.1


100 200 300


(f)


2.05[


2


1


.951


1.9


100 200 300


Figure 4.5: The ratio of the detector efficiency estimates and true detector efficiencies, where the detector efficiency estimates are generated using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.


63


2.1 2.05


2


1.95 1.9


2.1 2.05


2


1.95 1.9







64


x


I I I .. _ -


X X
x


X
X X X X X
X


XX X XX X X XX X X X X XX X X XX X XX XXX X X
X


0. - - -- - - - - -- - .. - -


'9


S41MRSAAg~gAAAAAAARN~n n E q Elk*


5 10 15 20 25
MC run number


30 35 40 45 50


Figure 4.6: VR of the detector efficiency estimates in the third experiment versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.


- - - (a) Fan-sum algorithm
x (b) Ferreira's algorithm
(c) EM-CD
o (d) EM-FP
o (e) EM-CD2
(f) EM-FP2


I1.


1:






































x x X x x XX X
- X ~ XX x X XX


A
0


0


11 A A 13
0 13 A 01 6130 13 6 [03& 13 13

0 0 JR A
131, IRI! R 1 0
AD 00 !11, ;k 0 A
0
0 0 IR 00 10 00 1 O0 0 I&S, 'I
00 0 00 00 0
0 0 0 0


5 10 15 20 25
MC run numbei


30 35 40 45 50


- -- - (a) Fan-sum algorithm
(b) Ferreira's algorithm
(c) EM-CD
0 (d) EM-FP
o (e) EM-CD2
A (f) EM-FP2


X

X x x
XXXX X X X X X
XX XX XXX X X X X X X


65


Figure 4.7: Maximum relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.


1.11 1.1


III I I I I I I


1.09 -


1.08 1.07 1.06


2




E
E
CO


1.05 -


/
-~ .-, '.


1.04 1.03


1.02






66














1.0161

- - (a) Fan-sum algorithm x (b) Ferreira's algorithm
1.016- (c) EM-CD
0 (d) EM-FP
3 (e) EM-CD2
1.0159- A (f) EM-FP2


o 1.0158A A '








0 00 o 000 00 0 00
1.0157 - O




00 OO 00 0 OO O OO 9
<1.0156 x x x
00~ 0





5 10 15 20 25 30 35 40 45 50
MC run number

Figure 4.8: Average relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.
























I I I . .


X X XX X X X X X X X X X
-x x V


X

X XXX X x x X x X XSx X X XX XA


A
0


5 10 15


0


Q A A 0
IoA 0 I K 3


20 25 30 35 40 45
MC run number


Figure 4.9: Maximum relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm. Note that the detector efficiency estimates are normalized so that their mean equals the mean of true efficiencies.


67


0.05


0.045


0.04 .


0.035


- - - (a) Fan-sum algorithm
x (b) Ferreira's algorithm
(c) EM-CD
o (d) EM-FP
o (e) EM-CD2
A (f) EM-FP2


6

4)
0.03
a)
E = 0,025


0.02 L


0.015-


0.01


50


_























0.02


0.018


0.016 [


I I I I I I


0.014 -


0.012 0.01


XX X x x x x x x x X x x x x x X x x x x x xx x x x




-


0.008 -


0.006 0.004-


0.002Oo 0 0 0 00 0 00 00 0 01


I I


5 10 15 20


25 30 35 40 45
MC run number


50


Figure 4.10: Average relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm. Note that the detector efficiency estimates are normalized so that their mean equals the mean of true efficiencies.


68


- - - (a) Fan-sum algorithm x (b) Ferreira's algorithm
(c) EM-CD o (d) EM-FP
0 (e) EM-CD2
A (f) EM-FP2


Z)


a)






69


Figure 4.11: Synthetic thorax phantom.


















(b)


4 he


(d)


4w


'4


Figure 4.12: Reconstructed images of the thorax phantom in Fig. 4.11 at the 15A iteration of the EMML algorithm using: (a) true efficiencies; (b) no normalization;
(c) fan-sum algorithm; and (d) Ferreira's algorithm.


70


(a)


t 4,


(c)


Di






71


(a)


(b)


4'


I


(c)


(d)


I


Figure 4.13: Reconstructed images of the thorax phantom in Fig. 4.11 at the 15th iteration of the EMML algorithm using: (a) EM-CD algorithm; (b) EM-FP algorithm;
(c) EM-CD2 algorithm; and (d) EM-FP2 algorithm.






72


12000

0 (a) No normalization
- - - (b) Fan-sum algorithm x (c) Ferreira's algorithm
10000- 0 (d) EM-CD
o (e) EM-FP
0 (f) EM-CD2
(g) EM-FP2
8000


w
U) 6000



4000



2000



0
0 2 4 6 8 10 12 14 16 18 20
Iteration number

Figure 4.14: MSE of the reconstructed emission images versus the iteration number of the EMML algorithm using: (a) no normalization; (b) fan-sum algorithm;
(c) Ferreira's algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.






73


(a)


(b)


4 ~


(c)


4


(d)


;4~ I


K ~


Figure 4.15: Reconstructed emission images of plane 10 of a thorax phantom at the 20th iteration of EMML algorithm using: (a) the fan-sum algorithm; (b) Ferreira's algorithm; (c) EM-CD algorithm; and (d) EM-FP algorithm.


(a)


(b)


;~L;.~ fA~
~


Figure 4.16: Reconstructed emission images of plane 10 of a thorax phantom at the 20th iteration of EMML algorithm using: (a) EM-CD2 algorithm; and (b) EM-FP2 algorithm.


- 103L,






74


(a)


(b)


i~3)


(c)


AI


Figure 4.17: Reconstructed emission images of plane 10 of a thorax phantom at the 20th iteration of EMML algorithm using: (a) EM-CD3 algorithm; (b) EM-FP3 algorithm; and (c) no normalization (i.e., {J& = 1}).


Ilk
41














CHAPTER 5
CONCLUSIONS


Algorithms for computing ML estimates of detector efficiencies in PET have been proposed. The algorithms follow from the EM algorithm and three Poisson models that were developed for the PET blank scan data. The refined Poisson model (i.e., Poisson Model III) directly accounts for the Poisson nature of blank scan data, certain geometric properties of scanners, and detector penetration. The close resemblance between real data and the synthetic data generated according to the proposed models provides some justification for the models.

For each Poisson model, a closed form expression for the E step of the EM algorithm was derived. The corresponding M step was solved iteratively using the coordinate descent method and a fixed-point iteration. The detector efficiency estimation algorithms that were based on the coordinate descent method have the desirable property that the estimates were guaranteed to lie in the interval [0, 1]. The fixed-point iteration based algorithms do not have this property. However, they converge about 20 times faster than the coordinate descent based algorithms. All of the algorithms have the property that the corresponding log-likelihood function increase monotonically with increasing iterations.

Experiments using synthetic and real data were used to assess the proposed algorithms. For the synthetic data case, the following objective criteria were used: (1) variance of the ratio of the efficiency estimates and true efficiencies, (2) maximum relative error, and (3) average relative error. In terms of these criteria, the proposed algorithms greatly outperformed the fan-sum and Ferreira's algorithms. However, the reconstructed emission images for all of the estimation algorithms were visually iden-


75






76


tical. In addition, the mean-square error (MSE) values associated with the emission images were also similar, but the proposed algorithms did produce the smallest MSE.

There is no accepted approach for comparing the performance of detector efficiency estimation algorithms when real data are used. We expect the EM-CD3 and EM-FP3 algorithms to have a large impact because Poisson Model III takes the effect of detector penetration into consideration. The reconstructed emission images for the fan-sum, Ferreira's, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are indistinguishable. When applying the EM-CD3 and EM-FP3 algorithms to reconstruct emission images, a different probability matrix P is used. The reconstructed emission images for the EM-CD3 and EM-FP3 algorithms show details much more clearly than the reconstructed emission images for the other estimation algorithms. While it is difficult to objectively determine whether using the EM-CD3 and EM-FP3 algorithms lead to better reconstructed emission images than using no normalization, we conclude that the effect of detector penetration is a significant factor when reconstructing emission images.














CHAPTER 6
FUTURE WORK


While the proposed algorithms converge in practice, the convergence properties of the proposed algorithms are not well understood. Theoretic analysis of the proposed algorithm will be pursued to determine their convergence properties.

As we stated in previous chapters, the performance of the proposed algorithms depends on the initialization step because the objective functions are not concave. It is a common practice to use a constant value as the initial detector efficiency estimate for iterative algorithms (e.g., Ferreira's algorithm [17]). From our experiments, we found that the average values of the detector efficiency estimates obtained from the proposed algorithms are very close to the constant value that was used to initialize the algorithms. However, if the constant value is scaled, detector efficiency estimates are not correspondingly scaled. Therefore, the effect of the initialization step must be investigated.


77














APPENDIX A
CONCAVITY OF THE LOG-LIKELIHOOD FUNCTIONS


A.1 Poisson Model I The log-likelihood function for the blank scan data is given by


( log P(B = b; ,A) = E [-EkEAp(k,) + bkl log(EkEAp(k,1))] k iefank

= [--EkElAd + bkl log(EkcjAd)]. (A.1)
d {(k,1):p(k,1)=d}

To determine whether a function is concave or not, we check if its Hessian is negative definite. To determine if a hermitian matrix is negative definite, we use the following theorem [43].


Theorem 1 The following condition is necessary and sufficient for an n x n hermitian H to be negative definite. The principal minors consisting of the the determinants of the k x k matrices in the top left-hand corner of H (k = 1, ... , n) are all negative:


Hi1 H12 H13
H11 < 0, Hil H12 <0, H21 H22 H23 <0,--H21 H22
H31 H32 H33

The (i, j) element of Hessian for the log-likelihood function is


Hi 201 (B = b; 0)(A2 H =(A.2)


78






79


where 0 = [r, A]. The first partial derivatives of 01(E, A) are


/1 (E, A)
O k


bki
- [EI (k,l) + -] k 1,2,... I,


and


D1(E, A)
Dd


bki
Z: [ f f +--], d=di,d2,...,d.
{(k,I):p(k,I)=d} d


And, the second partial derivatives are


if k= ifl1E fank


otherwise


: [q6(p(k, 1) - d)],
Iefank


E [- W]
{(kL):p(kI)=d) d
0


a2,p1(OE,A) _aAd8fka2A) - aBdaAm-


From (A.3), we know that


1
Hil = 2 E bil < 0,
i lEfan,


H12 = H21 = 0,


and


1
H22 = 2 E b2 < 0.
E2 IEfan2


Therefore,


Hn1 H21


H12 H22


- Efank


0


(A.3)


if d=m

otherwise


= Hn H22 - H12H21 > 0.






80

So we conclude that H is not negative definite and that the objective function of

(P1), 11(E, A), is not concave.

For the over-parameterized approach, the same conclusion can be reached by a similar derivation. The log-likelihood function for the blank scan data is given by


log P(B =b; e, A) = log 1I e-'fkIA(k,) (EklAp(kl))bkI
k LEfank bkl!

= log 1i f e'EklAd (klAd)bkl
d {(k,L):p(k,I)=d} bkI! = I Z [-EkIAp(k,L) + bkI log(EkIAp(k,)) + log(bkl!)] k lEfank

=E E [-EIkAd + bkI log(fkAd) + log(bkl!)]. (A.4)
d {(k,I):p(k,I)=d}


The first partial derivatives of log-likelihood function are


a log P(B = b; e,A) bkl
Ap(k,) , k = 1,2,...I, l E fank,
OEki EkI


(A.5)


and


0 log P(B = b; e,A)
OAd


(A.6)


z [- kI + -], d = di, d2, ..., dD
{(k,I):p(k,I)=d} d


Further, the second partial derivatives are


02 log P(B = b; e, A)
DEkI O6mn

02 log P(B = b; e, A)
o9Adg=k


92 log P(B =b; e,)aAdaAm,


-+ if (k, 1) (m, n) {0 if (k,l1) $(m, n)

k- if p(k,1)d 0 if p(k, 1) d ' [-b I if d = m {(k,I):p(k,i)=d} d
0 if d # m


(A.7) (A.8)


(A.9)






81


From (A.7)-(A.9), it can be seen that


H11 < 0, H12 = H21 = 0, H22 < 0.


Therefore, Hi1 H21


H12
= H11H22 H22


- H12H21 > 0.


So we conclude that the log-likelihood function is not concave.


A.2 Poisson Model II

The log-likelihood function for the blank scan data is given by


02(E, A)


log P(B = b; e, A)


E E [-yEkAc(k, 1) + bkl log(EkqEAc(k, 1))]. k of#an a

The first partial derivatives Of 0k2(c, A) are


902(e, A)
= E


be
Z [-EiAc(k,1) + -], k = 1,2, ..., I, icfan, E


and


02 (4E,A) IbkI
[A + A], k 1, 2,..., I.
k=1 iEfan,


(A.10)






82


And, the second partial derivatives are


I-


,2 2(E,) _. a~aek


- E if k = l
LEfank f
-Ac(k, 1) if I E fank


0


otherwise


Z Elc(k, 1), lEfank


92 02 (E,A)= I
OA2 _1 b
k=1 lEfank \


From (A.11), we know that


1
H1= - (
E1 LEfan1


bil < 0,


H12 = H21 = 0,


and


1
H22 = 2 E b2 < 0.
62 lEfan2


Therefore,


Hi1 H21


H12 H22


= H11H22 - H12H21 > 0.


So we conclude that H is not negative definite and that the objective function of

(P2), 02(E, A), is not concave.



A.3 Poisson Model III The objective function of the problem (P3) is


log P(B = b; p, A)


(A.11)


03 (p, A)






83


= >1 [-pkpAa(k, 1) + bkl log(pkpAa(k, 1))]. k iefank


(A.12)


Because Poisson Model III has the exact same functional form as Poisson Model II, the second derivatives of 03(p, A) are


- N if k = I
lEfank Pk
-Aa(k, 1) if 1 E fank

0 otherwise

- pia(k,1), lEfank
I g k k=1 Efank


(A.13)


From (A.13), we know that


1
Hil = - 1 bit < 0,
i lEfan,


H12 = H21 = 0,


and


1
H22 = 2 E b2t < 0.
P2 IEfan2


Therefore, Hil H12 = HlIH22 - H12H21 > 0. H21 H22 So we conclude that H is not negative definite and that the objective function of

(P3), 43(p, A), is not concave.


-9Pk 9P1


82 3(P,A)
aOPkp


8202(p,A)
i9A2














APPENDIX B
CRAMER-RAO LOWER BOUNDS To compute the Cramer-Rao lower bound, the Fisher information matrix must be computed. In this appendix, we derive the Fisher information matrix for the sufficiently-parameterized approach associated with each Poisson data model.


B.1 Poisson Model I For Poisson Model I, the probability density function (PDF) is


Ii - (kEjApck,j) ) bk
P(B =b; e, A) = J 1 e-k(k) bkl!
k=1 wt fan k

From (A. 3), we know that the second partial derivatives of log P (B = b; E, A) are


- Z if k =1 021ogP(B=b;4E,A) - IEfank '
Jfgf -A(k,l) if 1 E fank
0 otherwise

02lo P b-A) - E [c J(p(k, 1) - d)], lEfank
E [--] if d = m
821og P(B=b;E, - {(k,1):p(k,1)=d} d
a)daAm { 0 otherwise


(B.1)


84






85


The expected value of bkl is EkcA(k,L), so we can calculate the expected values of the second derivatives as follows.


E [ o g P(B=b;E,A)
t9Ckafl Ib=B]


I


a2og P(B=b;E,A) a adaEk b=B] -


E[9021gP(B=b;EA)
eAdOAm


- E A(kI) if k=l
IEfank
-Ap(k,l) if 1 E fank

0 otherwise

E [c,6(p(k, 1)-d)],
IEfan,


-]if d= m {( k,):p(k,I)=d} d 0b=Be
0 otherwise


Consider a simplified example, where there are only 4 detectors on the ring (Fig. B.1). Detector pairs (1,3) and (2,4) have the same perpendicular distance from the




2 1









3 4



Figure B.1: A simplified example to derive the Fisher information matrix.


origin and therefore have the same A, as the mean of the transmitted photon count. Detector pairs (1,4), (2,1), (3,2), and (4,3) have the same perpendicular distance from the origin and therefore have the same A2 as the mean of the transmitted photon count. The data {b13, b24, b14, b21, b32, b43} are collected. The problem is to estimate


(B.2)






86


{I1, 62, E3, E4, A,, A2}. Then the Fisher information matrix can be found to be


(6, A) -


T1 A2 A1

A2 C3 E2 + C4


A2

T2 A2 A1

C4 El + E3


A1 A2

T3 A2

E1

E2 + E4


A2 A, A2

T4 E2 El + 63


E3 f4 61 C2

T5

0


E2 E4 C1 + E3 E2 + E4 El + E3

0
624


where


1
T, = - (2 + 014 + 042),
El


E2 1
T3 = --(E1A1 + 62A2 + 4A2),
E3

T4 = -(iA2 + f2A, + 32),
64

T5 = 1 (Ei3 + 62E4),



T6= (ClE4 + C2C1 + E3(2 + f4E3).


and


We proved that the determinant of I(E, A) is equal to zero. That is, I(E, A) is noninvertible. The same conclusion can be extended to all higher-dimensional cases.



B.2 Poisson Model II For Poisson Model II, the PDF is


P(B = b; E, A) = 11 1J e -kfIAc(kL) (CkCAc(k, k=1 Efank bkl!






87


From (A.11), we know the second partial derivatives of log P(B = b; E, A) are


if k = 1


-
a2 iogP(B=bE,A) _ IEfank '
99fk(9f-Ac(k,l1) if

0 o
a21og P(Bb; EA) =-c(k, 1),
lEfan
.21og P(B=b;tA) I=b
2 E2 *
k=1 lefank


thErfank therwise


(B.3)


The expected values of the second partial derivatives of log P(B = b; E, A) are


eJAc(k,i)

E [ 0210g (B=b;E,,\) Ifn,
l OE, Ib=B] = -Ac(kl)

0
E [ alog P(B=b;E,) Ib=B = - Z cjc(k, l), le fan ,
E l aog P(B=b;CE,) I k)
[ A2 Ib=B] = - L L c
k=1 Efan,


if k = 1 if I E fank otherwise


(B.4)


Consider the example in Fig. B.1 again. There are two possible values for {c(k, l)}: ci for detector pairs (1,3) and (2,4), because these two detector pairs have the same perpendicular distance from the origin; c2 for detector pairs (1,4), (2,1), (3,2),and (4,3), because these four detector pairs have the same perpendicular distance from the origin. Therefore, the Fisher information matrix is


I(E, A) =


Si Ac2

Ac, Ac2 S2


Ac2 S3

Ac2

Ac, S4


Ac, Ac2 55

Ac2 S6


Ac2 Ac2 Ac2 S7

S8


S2

S4 S6 S8 59






88


where

S = -(E3c1 + E2C2 + E4c2), El

S2 = 1(E3C1+ E2c2 + E4c2),
A
S3 = -(Eic2 + E3C2 + E4c1), E2

S4 = (EiC2 + E3c2 + E4c1),



A
S5 = -(Eici + (2c2 + Eac2), E3




S7 = (61C2 + e2c1 + E3c2),
E4

S8 = 1 ElC2 + E2Ci + E3C2),

and

S9 = EEiCc + E2E4c1 + EIE4C2 + E2fic2 + f3E2c2 + E4E3c2).

It can be readily proved that, regardless of the values ci and c2, the fifth row of the matrix I(E, A) can be expressed as a linear combination of the first four rows. Therefore, the Fisher information matrix I(E, A) has a determinant equal to zero and is non-invertible. The same conclusion can be reached for the general case.



B.3 Poisson Model III Because Poisson Model III has the exact same functional form as Poisson Model II, we have the same conclusion for Poisson Model III regarding the Fisher information matrix. In other words, the Fisher information matrix is non-invertible.














APPENDIX C
DERIVATION OF THE E STEP OF THE EM ALGORITHMS


In this appendix, we prove that the E step is given by (2.9). Said another way, we determine the conditional expectation of the complete data, X = [B, N], given the incomplete data, B, and current estimates of E and A.


C.1 Poisson Model I


Using Bayes' Rule [44], the log-likelihood function of the complete data is given


by


L,(b, n; E, A)


- logP(B = b,N = n;c,A)

= logP(B = b I N = n; E, A)P(N = n; E,A)
I
= log) JJ (b)(lkEi)bkl(1 - EkEL)nk~-bk
k=1 IEfank
dD nki
+ log e- 1
d=d1 {(k,1):p(k,1)=d} nki= S S [log(nk|) + bkl log(Ekl) + (nkt - bkl)log(1 - EW)1
k=1 lEfank
dD
+ 5 E [-Ad + nki log(Ad) + log(nki!)]. (C.1)
d=di {(k,1):p(k,1)=d}


89






90


Let U(E, A; r', Am) denote the conditional expectation of the log-likelihood function of the complete data given the incomplete data and estimates E' and A". Then,


U(E, A; E', Am)

i E[L,(B, N;,E, A) I B = b; E, Am],
I Z E(NkI )IB
= E E (Ef an B = b; Em, A"] + log(EkE)E[Bkl I B = b; Em, A']
k=1 IEfan, (C.2)
+ log(1 - EkQI)E[Nki - BkI I B = b; Em, Am])
dD
+ E E (-Ad + log(Ad)E [N I B = b; e', Am]
d=dl {(k,):p(k,i)=d}
+ E[log(Nki!) I B = b; E", Am]).


To maximize U(E, A; E', Am) with respect to E, we take into consideration only the terms log(EkfI)E[BI | B = b; E', A'] and log(1 - EkE)E[N, - BkI I B = b; E', Am] as the remaining terms are independent of E. Similarly, we only need to consider the terms -Ad and log(Ad)E[N, I B = b; em, Am] in order to maximize U(E, A; E', Am) with respect to A. Let n' A E[Nk, I B = b; em, Am]. Because E[BkI | B = b; E"', A"'] = bk1,


U(,E, A; Em, AtM)

= i Z [bka log(fkEl) + (n' - bkL) log(1 - ekEl)] (C.3) k=1 tsfan(
dD
+ Z Z [-Ad + nm log Ad]+ C1, d=dl {(k,1):p(k,I)=d}

where C1 is independent of E and A. To finish the E step, an expression for nQ must be obtained.






91


To find n', we need to derive the conditional probability of Nk, given the incomplete data.


The second equality therefore


P(Nk = njk B = b; E, A)

P(Nkl = nkI BkI = bkI; E, AN)
_ P(Bkl=bkl INkL=nkI)P(N l=nkl)
P(Bkl=bkl) (C.4)
- kkj' )(Ekc)bkl (1-EkEl)"kl-bkle p(kI) nkl
ek l) p(k,L) (CkeLAp(k,I))bkl
e bkl!
-(1-fkl)Ap(k,) [(1 EkEL)Ap(k,l) Iki-bkl (nkI--bkl)!
follows from Bayes' Rule. The conditional expectation of NkL is


E[Nk, I BkI = bkI; E, A] E nkIP(NkI = nkI | BkI = bkl; C, A)
l!bkl
= E (nki - bkI + bkI)P(Nk= nkI I BkL = bkI; E, A)
nkI !bkj
= E (nkI - bkl)P(Nkl -- ni BkI = bkI; E, A) (C.5)
nklk bkl
+ E bkIP(NI = nkI I Bk1 = bkI; C, A)
nkL bkl
E n je-(1-EkfI)Ap(k,) [(1--kf)Ap(k,)]|"k nk, > >0 nk
+ bkL Z: e-( 1-kE)Apk,1 [(1-kE)p(k ,)]nhk, n'k;>OI
where n'l = nkI - bki. The first summation is equivalent to the mean of a Poisson random variable with parameter (1 - EkI)Ap(k,1), and the second summation is the sum of the probability over all possible values of n'. Therefore,


E[NkI I BkI = bkI; E, A\] = (1 - fkEI)Ap(k,) + bkl.


(C.6)


Given the previous estimates Em and A',


n F = E[Nk, I BkI = bkI; Em, Am] = (1 - EE7)A + bk,


(C.7)






92


for k = 1, 2,..., I, 1 E fank.


C.2 Poisson Model II With c(k, 1) known for k = 1, 2, ... , I, I E fank, the log-likelihood function of the complete data X = [B, N] is given by


L'(b, n; E, A)


= log P(B = b, N = n; e, A)

= logP(B= bI N = n;E, A)P(N = n;e, A)
I
= log k (")(Eke Cbk1(1-ekel)kI
k=1 IEfank

+ log y e-Ac(k,1) (Ac(k, 1))nki
k=1 I fank fkL!

S)[log((| + bkI log(Ekl) + (nkI - bkl)log(1 - Eke)]
k=1 LEfank

+ E [-Ac(k, 1) + nki log(Ac(k, 1)) + log(nkl!)]. (C.8)
k=1 IEfank


Let U'(E, A; E", A"-) denote the conditional expectation of the log-likelihood function of the complete data given the incomplete data and estimates E' and Am. Then,


U'(E, A; r'm, Am)

A E[L'(B, N; E, A) B b; E', A'],

= ean [() B =b; E, A'] + log(Ekeq)E[Bkl I B = b; E', A']
k=1 icank
+ log(1 - Ek()E[Nk, - BkI | B = b; E", A"-])
I
+ L E (-Ac(k, 1) + log(Ac(k, 1))E[Nk, I B = b; Em, Am]
k=1 LEfank
+ E[log(Nk!) I B = b; Em, A'm]).


(C.9)


To maximize U'(E, A; E', Am) with respect to E, we take into consideration only the terms log(EkqE)E[Bkj | B = b; Em, Am] and log(1 - EkEI)E[Nk, - BkI I B = b; Em, Am]






93


as the remaining terms are independent of c. Similarly, we only need to consider the terms -Ac(k, 1) and log(Ac(k, 1))E[Nk, I B = b; ', A'] in order to maximize U'(E, A; E', A') with respect to A. Let nm = E[Nk, I B = b; Em, A'm]. Because E[Bl I B = b; Em, A'] = bk,


U'(E, A; Em, Am) = [bkI L log(eE) + (nm - bkl) log(1 - EkEL)] (C.1O)
k=1 g~an kI C.10
I
[-Ac(k, 1) + nm log Ac(k, 1)+ C2, k=1 lEfank


where C2 is independent of E and A. To finish the E step, an expression for nm must be obtained.

To find n', we need to derive the conditional probability of Nkl given the incomplete data.
P(N, = ngk I B = b; f, A) = P(N, = nkL Bkj = bki; E, A)
_ P(BkL=bklINkl=nkI)P(Nkl=nka) P(Bki=bka) (C-11)

e-ekel'Wck,') (e'-f"\c(k,'))bk
- ~bkI,~E~E/b nkl

-(1--kEL)Ac(k,) [(1-eke)Ac(kL)|"kI-bk1
(nkj-bkj)!




Full Text

PAGE 1

MAXIMUM LIKELIHOOD ESTIMATION OF DETECTOR EFFICIENCIES IN POSITRON EMISSION TOMOGRAPHY i i . ) ! | By AQr > | WEN-HSIUNG LEE 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 2001

PAGE 2

Copyright 2001 by Wen-Hsiung Lee

PAGE 3

ACKNOWLEDGMENTS First I would like to thank my advisor and mentor, Professor John M. M. Anderson, for his endless effort and support on my research. His encouragement and suggestions have helped me through all the difficulties. More importantly, he taught me the methodology of solving problems, and he showed his diligence and patience in research. In short, he set a good example to me, which I believe will benefit me for all my life. Sincere thanks go to my committee members, Dr. Edmonson, Dr. Li, and Dr. Pardalos, for their guidance and patience. Sincere thanks also go to Dr. Vatow for his constructive comments and invaluable help. I also feel thankful toward Hong and Kerkil, my colleagues, who made me enjoy a pleasant teamworking environment. I owe a large measure of gratitude to my mother and my sisters for their unwavering support to help me go through tough times. Without them, I would not have had the strength to go on with my study at UF. in

PAGE 4

TABLE OF CONTENTS ACKNOWLEDGMENTS iii ABSTRACT vi CHAPTERS 1 INTRODUCTION 1 1.1 Overview of PET 2 1.2 Basic Physics of PET 4 1.3 Errors in PET 5 1.3.1 Non-Collinearity 5 1.3.2 Positron Range 5 1.3.3 Attenuation 6 1.3.4 Accidental Coincidences 6 1.3.5 Scatter 6 1.3.6 Detector Deadtime 7 1.3.7 Detector Inefficiency 7 1.4 Mathematical Model for PET 8 1.4.1 Poisson Model for Emission Scan Data 8 1.4.2 Modified PET Model for Emission Scan Data 9 1.5 Existing Methods for Estimating Detector Efficiencies 10 1.6 Outline of Research 13 2 ML ESTIMATION ALGORITHMS USING A POISSON DATA MODEL . 16 2.1 Poisson Model I 16 2.1.1 Estimation Problem Associated with Poisson Model I 21 2.1.2 Estimation Algorithms associated with Poisson Model I . . . . 21 2.2 Poisson Model II 30 2.2.1 Estimation Problem Associated with Poisson Model II ... . 33 2.2.2 Estimation Algorithms Associated with Poisson Model II . . . 34 3 ML ESTIMATION ALGORITHMS USING A REFINED DATA MODEL . 38 3.1 Poisson Model III 38 3.1.1 Estimation Problem Associated with Poisson Model III ... . 44 3.1.2 Estimation Algorithms Associated with Poisson Model III . . 44 IV

PAGE 5

4 SIMULATION 48 4.1 Synthetic Data 49 4.1.1 Comparison of Detector Efficiency Estimates 50 4.1.2 Comparison of Reconstructed Emission Images 54 4.2 Real Data 56 5 CONCLUSIONS 75 6 FUTURE WORK 77 APPENDICES A CONCAVITY OF THE LOG-LIKELIHOOD FUNCTIONS 78 A.l Poisson Model I 78 A. 2 Poisson Model II 81 A. 3 Poisson Model III 82 B CRAMER-RAO LOWER BOUNDS 84 B. l Poisson Model I 84 B.2 Poisson Model II 86 B. 3 Poisson Model III 88 C DERIVATION OF THE E STEP OF THE EM ALGORITHMS .... 89 C. l Poisson Model I 89 C.2 Poisson Model II 92 C.3 Poisson Model III 94 REFERENCES 96 BIOGRAPHICAL SKETCH 100 v

PAGE 6

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MAXIMUM LIKELIHOOD ESTIMATION OF DETECTOR EFFICIENCIES IN POSITRON EMISSION TOMOGRAPHY By Wen-Hsiung Lee December 2001 Chairman: John M. M. Anderson Major Department: Electrical and Computer Engineering Positron emission tomography (PET) is a medical imaging modality that enables physicians and researchers to study biochemical process within the human body and small animals. In PET, a subject is administered a radiopharmaceutical that is absorbed preferentially by the region-of-interest (e.g., brain). Due to the radiopharmaceutical and certain interactions at the cellular level, photon pairs are emitted in all directions and the number of photon pairs emitted from a region is proportional to the metabolic rate of the region. PET scanners are designed to detect and count the photon pairs. From this photon count data, images are reconstructed whereby the value of a pixel is proportional to the metabolic rate of the region associated with the pixel. The reconstructed images provide valuable information to help physicians detect abnormalities, such as tumors, which have extremely high metabolic rates. There are many sources of error in PET that must be addressed in order to obtain accurate images. The focus of this research is the problem of detector inefficiency. Due to inherent detector non-uniformity, detector efficiencies are less than 100%. To obtain more precise emission images, the photon count data must be corrected to account for the effect of detector inefficiency. vi

PAGE 7

We introduce a maximum-likelihood method for estimating detector efficiencies in PET. First, we develop three Poisson models for blank scan data obtained from rotating rod sources. Then, for each model, we estimate the detector efficiencies using expectation-maximization algorithms, where the maximization step is solved using two optimization algorithms. As desired, the resulting algorithms have the property that the log-likelihood function is non-decreasing as the iteration number increases. For each data model, one of the proposed algorithms guarantees that the efficiency estimates always lie in the interval [0, 1]. Although the second algorithm for each data model does not have this property, in all of the experiments performed it produced efficiency estimates that were between zero and one. Simulation studies using synthetic data demonstrate that, based on various comparison criteria, the proposed estimation algorithms outperform two alternative approaches. Additionally, simulation studies using real data show that the third Poisson model leads to much better reconstructed emission images than the first two models. Vll

PAGE 8

CHAPTER 1 INTRODUCTION Medical imaging modalities are widely used to obtain anatomical and biochemical information about the human body. For example, X-ray computed tomography and magnetic resonance imaging are used to acquire detailed images of the anatomical structure of a region of interest (ROI) (e.g., brain, heart, thorax). In many situations, however, it is more important to obtain functional information than anatomical information [1]. For instance, the knowledge of blood flow rates, oxygen utilization rates, and glucose utilization rates of brain or heart might be imperative for diagnosing certain diseases. Biochemical and functional information can be obtained via a class of imaging techniques known as nuclear medicine imaging. Among the techniques within nuclear medicine imaging are planar imaging [2,3], single-photon emission computed tomography (SPECT) [4-6], and positron emission tomography (PET) [7]. Planar imaging suffers from relatively low availability of radiopharmaceuticals, poor sensitivity of cameras, and, most importantly, the so-called projection problem. The projection problem refers to situations where the ROI might be obscured due to activity in front of or behind the ROI. Moreover, photons originating from the ROI might be attenuated due to overlying tissue. Although SPECT overcomes the projection problem and provides more sensitivity than planar imaging, it still has insufficient sensitivity and quantitative accuracy for certain applications. Furthermore, SPECT uses radiopharmaceuticals that incorporate long half-life radionuclides, such as 123 I and 201 TI. Because of the long half-lives 1

PAGE 9

2 (ranging from hours to several days) of the radionuclides, patients are malignantly exposed to radioactivity for a long period of time. PET avoids those afore-mentioned shortcomings and provides better images in the sense of efficiency, sensitivity and accuracy. In addition, positron-emitting isotopes of carbon, oxygen, nitrogen and fluorine can be readily incorporated into many useful radiopharmaceuticals because these compounds are naturally used by the body [8]. However, because radionuclides incorporated into radiopharmaceuticals for PET have short half-lives (2 minutes, 10 minutes, 20 minutes, and 110 minutes for 11 C, ls O, 13 N and 18 F, respectively), on-site cyclotrons are needed to produce radionuclides. The high cost and high complexity of cyclotrons have prevented PET from dominating the field of nuclear medical imaging. Even so, PET still has a wide range of use in clinical and research applications due to its advantages over other modalities. There are a number of error sources in PET, such as attenuation, accidental coincidences, and detector inefficiency. In this dissertation, we address the problem of detector inefficiency and propose algorithms, based on the maximum likelihood estimation principle, that estimate detector efficiencies. 1.1 Overview of PET A PET study begins with the injection or inhalation of a radiopharmaceutical by a subject. A radiopharmaceutical is a substance that is labeled with a positronemitting isotope. The radiopharmaceutical is transported to and absorbed by the ROI. When the isotope decays, it emits positrons, which travel a short distance (on the order of tenth of millimeter) before annihilating with electrons. The areas with a high concentration of the radiopharmaceutical emit more positrons than the areas of low concentration. The annihilation of a positron with a nearby electron produces two high-energy (511 keV) photons that move away from each other in nearly opposite directions. The position where the annihilation takes place is unknown and

PAGE 10

3 the photonsÂ’ line of flight is omni-directional and completely random. If these two photons, after traveling through the subject, are detected by a pair of detectors on a detector ring surrounding the ROI within a small timing window r (~ 10ns), a coincidence is recorded by that detector pair or tube (a detector pair is also referred to as a tube.) The observed emission data are the collection of the counts (i.e., total number of coincidences) recorded by each detector tube. An illustration of the scanning process is provided in Fig. 1.1, which shows a cross-section of the scanner. Figure 1.1: Annihilation of a positron with a nearby electron produces a pair of photons that propagate in nearly opposite directions. The number of photons emitted from a region is proportional to the concentration of radiopharmaceutical in the region. Moreover, the concentration of the radiopharmaceutical in the region is proportional to the metabolic rate of the region. Using the emission data, reconstruction algorithms create images, called emission images , where the value of a pixel is proportional to the concentration of the radiopharma-

PAGE 11

4 ceutical in the region associated with the pixel. Thus, emission images can be used to assess biochemical processes. 1.2 Basic Physics of PET The cylindrical rings of detectors surrounding the ROI typically have a diameter of 80-100 cm and an axial extent of 10-20 cm. The detectors are shielded from outside the field of view (or ROI) by relatively thick, lead end-shields. Most scanners can switch between two operating modes. In the first mode, known as slice-collimated mode or 2-dimensional (2D) mode, axial collimation is provided by thin annular rings of tungsten called septa, and only direct plane sinograms involving coincidences between detectors on the same ring are recorded. In practice, the situation may be more complicated in that, in order to improve statistics, direct and cross plane sinograms can be combined in different ways [9]. The other mode is fully three-dimensional (3D) mode where septa are retracted and cross plane sinograms of coincidences between detectors on different rings are also recorded [7]. Detectors are critical components of a PET scanner [7]. Each detector ring usually contains several hundred detector modules, where each module consists of an array of crystals coupled to several photomultiplier tubes. When a photon interacts with a crystal, electrons are moved from the valence band to the conduction band. Then, the electrons return to the valence band at impurities in the crystal. This process causes numerous light photons to be emitted. A photomultiplier tube collects light photons and converts them into electrical signals, which are decoded so as to determine which detector module was actually hit by a photon. The material bismuth-germanate is widely used as the crystal because of its high attenuation density, high light output (2500 light photons per 511 keV photon), fast rise time, and short decay time (300 ns). Crystals with high attenuation densities are desirable because the probability of a photon interaction is very high. Moreover, a high light output improves positioning

PAGE 12

5 accuracy, while a fast rise time is helpful for accurate timing, and a short decay time enables the scanner to handle high counting rates. 1.3 Errors in PET There are two major categories of physical effects in PET: (1) effects due to the positron annihilation process and propagation of photons through the body, such as positron range, attenuation, Compton scatter, and accidental coincidences, and (2) effects that result from the interactions of photons with detectors, such as detector efficiency and dead-time [7]. 1.3.1 Non-Collinearity As mentioned previously, an emitted positron interacts with an electron as it travels through the body. Each interaction requires certain energy. When the momentum of the positronÂ’s energy is nearly zero, annihilation occurs and two annihilation photons are produced each with an energy of 511 keV. The photon pair flies off along nearly collinear paths, where the degree of non-collinearity depends on the energy left in the positron [7]. The energy of the positron after annihilation varies widely among isotopes. The divergence from collinearity is on the order of a degree or less. So it is usually ignored. 1.3.2 Positron Range The distance the positron travels before annihilation is termed positron range. Typically, the positron range is about 3-5 pm [7]. Since the resolution of reconstructed PET images is over 5 mm for most commercial scanners, the error introduced by positron range is not a serious problem. However, for high resolution experimental scanners, this becomes the fundamental limitation of PET.

PAGE 13

6 1.3.3 Attenuation Two interactions that may cause a photon pair to go undetected are photoelectric absorption and Compton scattering. For 511 keV photons, the incidence of photoelectric absorption is very low and can be ignored in most cases [7]. A Compton scatter happens when a photon interacts with an outer shell electron. In this case, the photon is deflected from its path. Consequently, the counts of the tube corresponding to the original path of a photon pair does not increase when one or both of the photons undergo Compton scattering. The total effect of photoelectric absorption and Compton scatter is called attenuation and undetected annihilations are termed scattered events. 1.3.4 Accidental Coincidences If two photons resulting from different annihilations are detected within the coincidence timing window, they will be incorrectly interpreted as a coincidence event. Such events are called accidental coincidences or randoms, as true coincidences or trues. The rate of accidental coincidences is related to the concentration of the radionuclide in the ROI. A common correction strategy is to subtract an estimate of the randoms by using an estimate of the single events or delayed coincidence circuitry [10]. 1.3.5 Scatter Suppose one or both photons of a photon pair undergo Compton scattering. If the photons remain in the detector plane and are detected, then the annihilation event is recorded by the wrong tube. In other words, since the photonsÂ’ path is not collinear, the event is mis-positioned and additive noise results. However, due to the energy loss caused by the scattering process, scattered photons can be distinguished from

PAGE 14

7 unscattered photons to a limited degree by setting a threshold such that photons with sufficiently low energy, which correspond to scattered photons, are ignored. 1.3.6 Detector Deadtime The time required to process a single photon event limits the counting rate of a PET scanner [11]. Event processing begins with the rising edge of the pulse for the first detector involved. The pulse is integrated for some time interval, after which position calculations and energy discriminations are performed. The detectors are “dead” to new photon events during this time. Since the number of randoms increase as the square of the activity in the field of view [7], the number of detected counts may decrease as the flux of incident photons increases. This effect is accounted for in two ways. First, it becomes a dominant limitation on the injected dose. Second, a scale factor is used to correct data for the deadtime. 1.3.7 Detector Inefficiency The efficiency of a detector pair is the probability that an incident photon pair is detected by the detector pair. Due to inherent detector non-uniformities, detector efficiencies are not equal and less than 100% [12, 13]. The efficiency of current PET detectors depends on a number of geometric and non-geometric factors. The nongeometric factors include fluctuations in photomultiplier tube gains and variations of crystal characteristics [7]. Geometric factors include the angles subtended by detectors and the angles of incidence of photons. The angle subtended by a detector pair depends only on the perpendicular distance from the origin, and the angle subtended by the detector pair affects the angles of incidence of photons. And the length of intersection of a detector and the path of a photon depends on the angle of incidence of the photon. If the intersection length is larger, the probability that the photon is detected by the detector is higher.

PAGE 15

8 1.4 Mathematical Model for PET Since the data in PET follow from a counting process, it is assumed that the emission process is a spatial Poisson process. In 1982, Shepp and Vardi [14] proposed a Poisson model for PET that has been widely accepted as a valid mathematical model for the PET process. However, their model was based on the assumption that the data are error-free except for Poisson noise. In the following, the basic model for PET emission data is first introduced. Because this dissertation focuses on how to correct (or normalize) emission data for variations in detector efficiency, the basic data model is then modified to take into account the effect of detector inefficiency. 1.4.1 Poisson Model for Emission Scan Data The ROI is divided into a grid of J voxels and the emission intensity is assumed to be constant over each voxel. Assuming the data are free of errors, the emission datum recorded by detector pair ( k,l ) is a realization of a random variable Y k i'. Yki ~ Poisson ([P*]fci), A: = 1,2, ...,/, I e fan*, where J [f*]H=EVj. (1.1) 3=1 and where I is the number of detectors on a ring and the set fan* is the set of detectors that are paired with detector k (see Fig. 1.2). In addition, P klJ is the probability that a photon pair emitted from voxel j is incident on detector pair (k,l), x 3 is the mean number of photon pairs emitted from voxel j , and x = [x\, x 2 > . . . , x j] T is the emission intensity vector.

PAGE 16

9 Figure 1.2: Definition of fan*;. 1.4.2 Modified PET Model for Emission Scan Data Again, the efficiency of a detector pair is the probability that an incident photon pair is detected by the detector pair. Let be the efficiency of detector pair (k,l). Having a detector pair that is not 100% efficient is mathematically equivalent to thinning [15] the incident photons. To see what effect the thinning has on a Poisson random variable, we make use of the following well-known Corollary. Corollary Let Z be a Poisson random variable with mean Z. Suppose a random variable Y is created by thinning Z with probability p. Then, Y is Poisson with mean pZ. Proof: OO P(Y = y) = J E p ( Y = y \Z = z)P(Z = z) z=y z=y

PAGE 17

10 a -2P u ^ (i-p) z y (zy = e y' 2 / ! ^ (* ~ y)\ Using the change of variable z' = z-y, the probability P(Y = y) becomes P{Y = V) = e ~zt ^ (1 -py\zy'+y y ! zho = c -z(pZ) y ^ ((i -P )zy z ’= 0 V'_ c -z (P-^) j/ c (i-p)z 2/! ( 1 . 2 ) Thus, y is a Poisson random variable with mean pZ. Q.E.D. According to the above Corollary, the emission data, after thinning with probability €m, become Y k i ~ Poisson(e kl [Px] kl ), k = 1 , 2 ,...,/, I e fan fc . (1.3) Given the observed emission scan data y, which is a vector with elements {yu}, the problem is to reconstruct the emission intensity vector x. Shepp and Vardi [14], the developers of the Poisson model for emission scan data, assumed that the detectors were 100% efficient. However, detector efficiencies must be estimated before emission images can be reconstructed, because, in reality, they are not ideal. 1.5 Existing Methods for Estimating Detector Efficiencies To estimate a set of detector efficiencies, a scan is taken using either rotating rods or a uniform cylinder that contains positron emitting nuclides. The scan, which is referred to as a blank scan (no patient in the scanner), lasts about 2-4 hours and is usually performed once a month (see Fig. 1.3).

PAGE 18

11 Figure 1.3: Blank scan. Numerous methods [9, 16-24] have been developed for estimating detector efficiencies. Unfortunately, existing methods for estimating the detector efficiencies are largely ad hoc and do not take into account the statistical properties of the blank scan data. For instance, an early technique discussed in [9] was to estimate the efficiency of a detector pair by computing the ratio of the number of photon pairs measured by that detector pair and the average number of photon counts per tube. It has been noted that the accuracy of this heuristic technique is limited by the statistical accuracy of the blank scan [9]. Another technique for estimating detector efficiencies is the well known fan-sum method by Hoffman et al. [16]. Hoffman et al. assumed that b ki , the number of photon pairs recorded by detector tube (k, /), is equal to bid = fcejnjfef, k = 1, 2, . . . , I, l e fan*, (1.4)

PAGE 19

12 where e k is the efficiency of the k th detector and n k i is the unknown number of photon pairs transmitted along detector tube (k, l). Hoffman et al. also assumed that 53 ti^ki ~ constant. (1.5) iefan fc Their reason for making the assumption in (1.5) is because the summation is over a large number of detectors and the variations would tend to average out. Under these two assumptions, the proposed estimate for the efficiency of the k th detector £k — 51 ( 1 6 ) /efan*. is approximately equal to h « e k • constant. (1.7) The approximation in (1.5) is poor when there are large fluctuations in the detector efficiencies [17] or the number of terms in the set fan* depends on k. In such instances, the fan-sum algorithm may not provide accurate estimates for the detector efficiencies. Also, the fan-sum algorithm does not account for the Poisson nature of the blank scan data. Several methods [9,16-19] have been developed based on the assumption of noisefree (i.e., no Poisson noise) blank scan data. For example, an iterative method by Ferreira et al. [17] is based on the assumptions that the data are noise free and every detector views the same activity: 53 n k i = A, k = 1,2,...,/, Z 6 fan*, iefan*. ( 1 . 8 )

PAGE 20

13 Substituting (1.4) into (1.8) gives efc = -r k = 1,2,...,/, I e fan fc . (1.9) iefan* € ' To get an approximate solution for the above system of equations, the following fixed-point iteration is used: e? +1 = E % fc= 1,2,...,/, iefan t , (1.10) Jefan*. * where e™ is the m th estimate of e*. The constant factor 1 /A is neglected because the efficiency estimates are normalized so that their mean is one. Like the fan-sum algorithm, Ferreira’s algorithm assumes noise-free data and relies on a key assumption (see (1.8)) that is questionable when the number of terms in the set fan* depends on k. Another algorithm based on the assumption of noise-free data is the algorithm by Defrise et al. [9]. The algorithm by Defrise et al. parallels the fan-sum method with the exception that geometric means rather than arithmetic means of the data are used to estimate the mean number of unobserved number of transmitted photon pairs. Defrise et al. claimed that their algorithm produces unbiased efficiency estimates in the sense that they converge to the real efficiencies in the ideal case where the blank scan data are noise-free. 1.6 Outline of Research All of the estimation methods we have discussed are sub-optimal in the sense that they assume noise free data and do not provide ML estimates. To capture the stochastic nature of PET data, we propose three Poisson models for blank scan data. At first we ignore the problem of detector penetration [25-28] when developing the

PAGE 21

14 first two Poisson models for blank scan data. However, later we incorporate the effect of detector penetration and refine our original Poisson models for blank scan data. Motivated by the optimality properties of ML estimators [29], we maximize the log-likelihood function associated with each Poisson model by using the expectationmaximization (EM) algorithm [30,31]. To apply the EM algorithm, we choose the blank scan data and unobserved number of transmitted photon pairs for each detector tube to be the complete data. Two possibilities have been explored: • Over-parameterized approach, where the unknowns are the efficiencies of the detector pairs, {e*,;}, and certain emission means. • Sufficiently-parameterized approach, where unknowns are the efficiencies of individual detectors, {ejt}, and certain emission means. Because the over-parameterized approach does not take into consideration the assumption that eki = ekCi, k = 1,2 ,...,/, I e fan*, the number of unknowns is much higher than that for the sufficiently-parameterized approach. Closed form expressions for the E step and M step of the EM algorithm are derived for our first model using the over-parameterized formulation. However, simulations show that the over-parameterized approach has poor performance. Thus, the results for the over-parameterized approach will not be presented in Chapter 4, which contains our simulation results. For each of our three proposed Poisson models for blank scan data, we use the sufficiently-parameterized formulation and develop two algorithms for estimating detector efficiencies. For each data model, a closed form expression for the E step of the EM algorithm is derived. Unfortunately, no closed form solutions to the maximization problems in the M step are available. For each data model, we propose solving the maximization problems in the M step iteratively by using a fixed-point iteration or the method of coordinate descent. All the resulting (EM) algorithms have

PAGE 22

15 the property that the log-likelihood function increases monotonically with increasing iteration. It is desirable that a detector efficiency estimator provides estimates that lie in the interval [0, 1]. The coordinate descent based algorithms we propose for each data model guarantee that the efficiency estimates always lie in the interval [0, 1]. The fixed-point iteration based algorithms for each data model do not have this property. However, in all the experiments we performed the efficiency estimates did lie in the interval [0, 1]. In terms of certain objective criteria, simulation studies using synthetic data demonstrate that the proposed estimation algorithms outperform the fan-sum and FerreiraÂ’s algorithms. Additionally, our experiments show that emission images based on the proposed detector efficiency estimates have smaller mean squared error than emission images obtained using the detector efficiency estimates generated by the fan-sum and FerreiraÂ’s algorithms. Simulation studies using real data show that the estimation algorithms based on the third Poisson model provide much better reconstructed emission images. This dissertation is organized as follows. In Chapter 2, we introduce two Poisson data models for blank scan data and develop estimation algorithms that provide ML estimates of the detector efficiencies. Then, in Chapter 3, we refine our data model to incorporate the effect of detector penetration and present the corresponding ML estimation algorithms. In Chapter 4, we discuss the results of our simulation studies. In Chapter 5, conclusions about our research are given. Finally, we suggest some topics for future work in Chapter 6.

PAGE 23

CHAPTER 2 ML ESTIMATION ALGORITHMS USING A POISSON DATA MODEL In this chapter, we first present two Poisson models for blank scan data obtained from rotating rod sources. Then, we present ML estimation algorithms based on these two Poisson models. The first model is based on the fact that the mean number of photon pairs transmitted along a detector pair only depends on its perpendicular distance from the origin, which is the center of field-of-view. The second model, which is a modification of the first model, takes into consideration that the rods are nearly uniform. Note, the second model has fewer unknowns than the first model. In many scanners, the space between blocks of detectors is filled up with lead. Throughout our research, though, we make the simplifying assumption that there is no lead between blocks of detectors. Also, we assume the detector rings are perfectly circular. 2.1 Poisson Model I Before discussing the first Poisson model, we introduce some notation. Let a denote the angle between the lines that connect the origin to detectors k and /, and p(k, l) denote the perpendicular distance from the origin to the line defined by the mid-points of detectors k and l (see Fig. 2.1). The detector indices k and l are defined such that k < /, and the angle a is defined counterclockwise starting from detector k and ending at detector l. Clearly, p(k,l) = rcos(f), where a = k) and r is the radius of the detector ring. Therefore, for detector pairs (ki,li) and (fa, l 2 ), p(fa, h) — p(k 2 , l 2 ) if h~ fa = l 2 fa. In addition, p(ki, h) and p(k 2 , l 2 ) will 16

PAGE 24

17 Figure 2.1: Definitions of p(k, l) and a. be equal if (/ x ki) + (Z 2 & 2 ) equals the number of detectors on one ring. To see this fact, consider Fig. 2.2, where detector pairs (1,114) and (50,321) of the SiemensCTI ECAT EXACT 921 scanner [26] are depicted (note: this scanner has / = 384 detectors per ring). Because on = g(114 1) and a 2 = §^(321 50), their sum Figure 2.2: An example to demonstrate how p(k, l ) is related to the value of l k. equals ot\-\-a 2 — ^(113 + 271) = 27t. Clearly Q! 2 +/d = 2n as well. Therefore, a.\ = /3, which implies p(l, 114) = p(50, 321). For the Siemens-CTI ECAT EXACT 921 scanner, Table 2.1 illustrates the detector pairs of the first projection and the perpendicular distance of each pair. From the

PAGE 25

18 Table 2.1: Perpendicular distances of detector pairs for the first projection of the Siemens-CTI ECAT EXACT 921 scanner. 1st projection member pair ( k , l) l-k p(k, l)(cm) 1st (152,264) 112 d x = 25.1114 2nd (152,265) 113 d 2 = 24.8428 3rd (151,265) 114 d 3 = 24.5726 80th (113,304) 191 d 80 = 0.3375 81st (112,304) 192 dgi = 0 82nd (112,305) 193 d 8 o = 0.3375 159th (73,343) 270 d 3 = 24.5726 160th (73,344) 271 d 2 = 24.8428 discussion above, it can be seen that the value of l — k determines the perpendicular distance p(k,l). Taking this property into account, it can be seen that there are 81 possible perpendicular distances from the origin. To see how the number of transmitted photon pairs is related to p(k,l), consider a scanner with 16 detectors, where the dashed circle stands for the path of the rods (see Fig. 2.3). It is easy to show that the arc of tube 2 is longer than the arc of tube 1. Thus, because the rods rotate at a constant speed, the rods will stay in tube 2 for a longer time than in tube 1. Observe that p(k , l ) of tube 2 is larger than p(k , /) of tube 1. More generally, it can be proved using a simple geometric argument that the dwell time is proportional to the perpendicular distance from the center to the tube. And, because the activity viewed by a detector tube is proportional to the dwell time of the rods in that tube, it is reasonable to claim that the mean number of the transmitted photon pairs depends only on the perpendicular distance from the origin to the tube. For a blank scan, let B^i denote the observed number of photon pairs detected by detector pair ( k , /), and N^i denote the unobserved number of transmitted photon pairs along detector pair (k,l). Given the discussion above, we assume that

PAGE 26

19 Figure 2.3: An example to show how the number of transmitted photon pairs depends on p(k, l). B kl and N ki are Poisson with means e k i\{k,i) and \ p ( k ,i), respectively: N ki ~ Poisson(Ap(i i /)), A: = 1,2, I e fan*, (2.1) B kl ~ Poisson {e kl \ p ( ktl) ), k = 1, 2, . . . , I, l e fan*, (2.2) where \ p ( k ,i) depends only on p(k,l). We refer to the data model described by equations (2.1) and (2.2) as Poisson Model I. To see the suitability of Poisson Model I, we compare real blank scan data and synthetic data generated according to Poisson Model I in Figs. 2.4 and 2.5. It can be seen from the figures that the synthetic data bears a close resemblance to the real data. Note, a rod passes through the same detector tube twice in each cycle. This fact is implicitly accounted for through the parameters {A p (* > q}. From Poisson Model I, it follows that E[B kt ] = e kl E[Nki). (2.3)

PAGE 27

20 50 L 1 ' 1 1 1 1 1 20 40 60 80 100 120 140 160 Figure 2.4: Overlapped 192 projections of real blank scan data. Figure 2.5: Overlapped 192 projections of synthetic blank scan data generated according to Poisson Model I. HoffmanÂ’s assumption in (1.4) follows by replacing the expectations E[B k i ] and E[N k i ] in (2.3) with their ML estimates b k i and n k i, respectively. The statement in (2.3) is more sound probabilistically than the assumption in (1.4) because detector efficiencies should be interpreted as probabilities. Moreover, all we can say about B kt and N ki is that their means obey the relationship in (2.3).

PAGE 28

21 2.1.1 Estimation Problem Associated with Poisson Model I For convenience, let B, N, e, and A be vectors with elements {N k i}, {e fc }, and {A p ( fc) /)}, respectively. Also, let b be a vector with elements {&*.*}, where b k i is an observation of B k i . To estimate e and A, we use the ML method: (PI) (e,A) = argmax 0, (2.4) e,A where (f> i(e, A) is the log-likelihood function *1 (e, A) = log P(B = b, e, A) = log II g ^kl^p(kyl) ( M ^ j k=l /Gfan fc bkh and 0 and 1 are /xl vectors of zeros and ones, respectively. In Appendix A, we show that Hessian of the log-likelihood function is not negative definite. Thus, the loglikelihood function is not concave. Given this fact, it is expected that the performance of iterative algorithms used to solve (Pi) will depend on the initialization procedure. Cramer-Rao Lower Bound To compute the Cramer-Rao lower bound (CRB) for the estimators of e and A in (PI), we need to first compute the Fisher information matrix. We prove in Appendix B that the Fisher information matrix is non-invertible for the probability density function associated with (Pi). Therefore, the CRB for the estimators of e and A cannot be determined. 2.1.2 Estimation Algorithms associated with Poisson Model I If n, an observation of N , was available, then the parameters e and A could be estimated using hi = — , k = 1 , 2 ,...,/, / € fan fc , n ki (2.6)

PAGE 29

22 A d — X/ ; {(fc,/):p(fc,0=d} , d — d\, rf 2 , . . . , d#, (2.7) number of pairs (&, /) such that p(k, l) = d where {dj} are the possible perpendicular distances and D is the number of possible perpendicular distances. Given this observation, our choice for the complete data is X — [B, N], It should be noted that an observation of N is not available. Consequently, as part of the EM algorithm, an estimate of an observation of N will be generated. The E step of the EM algorithm requires the computation of the conditional expectation of the log-likelihood function of the complete data given the incomplete data and current estimates of e and A. Let U(e, A; e m , A m ) denote this conditional expectation, where e m and A m are the estimates of the detector efficiencies and mean number of transmitted photon pairs, respectively, at the m th iteration. Also, let L c (b , n; e, A) denote the log-likelihood function of the complete data: L c (b, n; e, A) = log P(B = 6 , IV = n; e, A). ( 2 . 8 ) In Appendix C, we show that U(e, A; e m , A m ) = E[L C (B, Ne, A) | B = 6; e m , A m ] / = Y, S lo S ( e *T<) + ( n kl b ki ) log(l €fcC|)] k=1 iefan k d D + £ ]£ [—Ad + log + Ci, (2.9) d=d\ {( k,l):p(k,l)=d } where C\ is a constant that includes all of the terms that are independent of e and A, and is the conditional expectation of Nki given the incomplete data and the estimates e m and A m : nZ = E[N ki | B kl = b kl ; e m , A m ], k = 1, 2, 6 fan*. (2.10)

PAGE 30

23 We also show in Appendix C that n ki — (1 ~ e ™ e T)^™{k,i) + hi, k = 1 , 2 , G fan*. (2.11) The quantity is an estimate of the number of photon pairs transmitted along detector pair ( k,l ). From (2.11), it can be seen that n £} is estimated by adding an estimate of the number of photon pairs that are transmitted but not detected, (1 — er e r)^( *:,*)> t° the number of photon pairs that are transmitted and detected, hiWe recognize from (2.9) that U(e, A; e m , A m ) is decoupled in the sense that it can be written as the sum of functions /i and g that depend only on A and e, respectively. More specifically, U(e, A; e m , A"*) = /<">( A) + S < m »(e) + C u (2.12) where fi A) = S + nj’J log A,(], d=dl {(k,l):p(k,l)=d} (2.13) and £ (m) (e) = 5Z 53 lo s( e i + ( n ki hi) log(l €*€|)]. k=1 iefan k (2.14) From (2.12), it follows that the problems of maximizing U(e, A;e m A m ) with respect to e and A are decoupled. Summarizing our results, the steps of the EM algorithm are: Step I. ( Initialization ) Set m—0. Choose e° and A 0 s.t. 0 < e° < 1 and A 0 > 0. Step II. ( E step )

PAGE 31

24 «« = (!tk*?)Km + 6 '«> ^ = 1)2,...,/, / e /an,. Step III. (M step ) ( p A) Am+1 = argmax/ 1 (m) (A). A>o (Pe) e m+1 = arg max g^ m \e). v ' & o 0, d — d\, ^2) • • • > djr), (2.15) we take the partial derivative of f[ m \ A) with respect to each Ad and set it to zero: dfW(\) 1 L = E (-l+ 0 for all m, k and d. Unfortunately, a closed-form solution is not available for the maximizer of the function U(e, A; e m A m ) with respect to e (see (Pe))We iteratively solve this maximization problem using either the coordinate descent (CD) method [32, pages 283-287] or a fixed-point iteration. See Abatzoglou et al. [33] or Luo et al. [34] for information on the convergence properties of CD method.

PAGE 32

25 EM-CD Algorithm To solve (Pe), we first use the CD method where the function g is minimized iteratively on a coordinate by coordinate basis. More specifically, the desired minimizer, e m+1 , is given by solving the following maximization problems iteratively for k = 1,2,...,/: C k = arg max q ^ ° 0 m+l,i+l m+l,t e fc-l ) e fc+l ) rn+l,i e I ), (2.18) where e™ +1,l+1 is the estimate of e™ +1 at the (i + l) th iteration. A reasonable choice for the initial estimate of e m+1 is e m+1 ’° = e m . The iteration is terminated when successive estimates are virtually the same. In practice, the algorithm is terminated when I m+l,i+l _ m+l,i £(m+M * ) 2 < t 2 -^) k = 1 e k where 5 t is a small positive number. The iterate that results when the stopping criterion is satisfied is considered the solution to (Pe)To see how the maximization problems in (2.18) can be solved, let M e ) = * 2 , • , zi-i, e, zi+ 1 , . . . , z/), / = 1,2,...,/, (2.20) where z is a vector whose elements, z\, z 2 , . . . , zj, all lie in the interval [0, 1]. It is straightforward to see that d 2 0 for all m, A: and d. Hence, i is a strictly concave function if ejj 1 e[0, 1] and \™ > 0 for all m, A: and d. Because of this concavity property, the problems in (2.18) can be readily solved with arbitrary

PAGE 33

26 accuracy using 1-D optimization algorithms (e.g., dichotomous search method, golden section method, and Fibonacci search [32, pages 268-275]). The constraint that every efc lies in the interval [0, 1] can be easily enforced by using the interval [0, 1] as the initial interval of uncertainty for the chosen 1-D minimization method. Hence, the algorithm is guaranteed to provide detector efficiency estimates that lie in the interval [0, 1] when e Q k £[0, 1] and A]] > 0 for all k and d. We refer to the EM algorithm that results from using the CD method to solve (P € ) as the EM-CD algorithm. To solve the maximization problem in (2.18), we use the dichotomous search method [32]. EM-FP Algorithm The second way we solve the problem (P e ) is to use a fixed-point iteration. Although it guarantees nonnegative detector efficiency estimates, the fixed-point iteration is not theoretically guaranteed to have the property that the estimates all lie in the interval [0, 1]. Nevertheless, in all of our simulations, the detector efficiency estimates obtained using the fixed-point iteration did have this property. An advantage of the fixed-point iteration is that it converges faster in practice than the coordinate descent algorithm. Ignoring the constraint in (Pe), our approach is to use a fixed-point iteration to solve the system of equations = 0, k = 1,2,...,/. (2.22) %7 (m) (c) dt k Taking the partial derivative of g ( m ' ) with respect to e k gives dg {m) (e) de k E [7 + (»B /€ fan, k — hi ) 1 ~ e fc e/ (2.23)

PAGE 34

27 Thus, the resulting system of equations are E [7 i + («S-W r =5 -]=0, * = 1,2,...,/. (2.24) K fan, e ‘ 1 e * e ' Obviously, a closed form solution for e m+1 is not possible. To find an approximate solution, we write (2.24) as E + E fefan* k k 1 idan k 1 CfcCi ^ = 1 , 2 ,...,/. (2.25) Simplifying the left-hand side of (2.25) yields £ £ T^7T> k = 1,2, . . . , i iefan*. iefan fc (2.26) The following fixed-point iteration is then used to find an approximate solution to the system of equations in (2.26): m+l,i+l £ fefan* r-HTft m+l,i L c i k — 1,2, ... ,1. ie fan* 1 * k m+J,i m+1,1 (2.27) As with the EM-CD algorithm, the initial estimate of e m+1 is € m+1, ° = e m . Further, the solution to (P e ) is the iterate obtained when the stopping criterion in (2.19) is satisfied. If n™ t in (2.11) is positive for all (k,l), then every term in (2.27) is nonnegative. Thus, the detector efficiency estimates will be nonnegative provided 0 < e m < 1 and A m > 0 for all m. We refer to the EM algorithm that results from using the discussed fixed-point iteration as the EM-FP algorithm.

PAGE 35

28 Over-Parameterized Approach We now assume the unknowns are the set of parameters {e fci }. This case is referred to as over-parameterized, case because the number of unknowns in this case exceeds the number of data points. Let e be the vector with elements {e kt }. The conditional expectation of the loglikelihood function of the complete data given the incomplete data and the current estimates of e and A is U 0 {e, A; e m , A m ) = £ ^ [b kl \og(e kl ) + (njj} b kt ) log(l e w )] k=1 iefan k + [~ x d + n ki log A d ] + Co , (2.28) d=d\ {( k,l):p(k,l)=d } where C Q includes all the terms that are independent of e and A. The function U 0 (e , A; e m , A m ) is decoupled: U 0 (e, A; e™, A m ) = g^\e) + /<"*>( A) + C\ (2.29) where «i m) (e) = E E log(e u ) + K b u ) log(l e„)] (2.30) fc=1 iefan k and /i m) (A) = E E [-Ai + nSlogAj). (2.31) d—di {( k,l)'p(k,l)=d } Because the function U Q (e, A; e m , A m ) is decoupled to two functions g( m \e) and fo m \ A), the problem of maximizing U 0 (e, A; e m , A m ) with respect to e and A is decoupled. The E step for the over-parameterized formulation is the same as the E step for the sufficiently-parameterized formulation except that e k ei is replaced by e k i. In other words, for k — 1, 2, / Gfan*,, the conditional expectation of N k i given the

PAGE 36

29 incomplete data and the estimates e m and A m is n kl & E[N kt | B u = b u ; e m , e m ] = (1 ejj) + b kl . (2.32) However, the M step differs. New estimates for e and A are obtained by maximizing U a (e, A; e m , A m ) with respect to e and A under the constraints 0 0, d = di, c?2 • • • ) dp. (2-34) Ignoring the constraints (2.33) and (2.34) momentarily, we take the derivatives of 5o m) ( e ) and /i m) (A) with respect to every e kt and X d , respectively, and set them to zero. 3o (m) (e) 1 -1 + — 0, * = 1,2,...,/, Ze fan*. (2.35) otki e k i 1 e kl -kl d/i m) ( A) dA d F, (— l + n^|— )— 0, d — di,di, . . . ,d D . {( k,iy. P (k,i)=d } (2.36) Solving (2.35) and (2.36), we have the following iteration. ^ +1 = -S, * = 1,2,...,/, Ze fan*. n kl (2.37) n {( k,iy. P (k,i)=d } fcZ ^ 771+1 _ number of pairs ( k , /) such that p(*, /) = d ——7, d — d\, d 2 , . . , dp. (2.38) For the constraints (2.33) and (2.34) to be satisfied, e kl must lie in the interval [0, 1] for * = 1, 2, . . . , I, l £ fan* and X° d must be nonnegative for d = d\, d 2 , . . . , d D . As mentioned earlier, the results using this approach are poor. Consequently, we do not include this approach in our simulation studies.

PAGE 37

30 2.2 Poisson Model II As we mentioned earlier, the CRB for the estimators of e and A in (PI) cannot be determined because the Fisher information matrix for the PDF in (PI) is noninvertible. One of the possible approaches to changing the non-invertibility is to reduce the number of unknowns in the estimation problem [35]. In this section, we present a second Poisson model that incorporates a priori information about the blank scan data so that the number of unknowns in the estimation problem is reduced. Figure 2.6: The exact way to calculate \ p (k,i), given A(r rod cos0,r rod sin0). For detector pair ( k,l ) (see Fig. 2.6), the mean number of photon pairs that are transmitted along the detector pair, X p ^,i), can be expressed as [37] / (2.39) rod path where A (x,y) is the number of photon pairs emitted from the rod at point (x, y) and P(k , l | x, y) is the probability that a photon pair emitted from the rod at point

PAGE 38

31 (x,y) is incident on detector pair ( k,l ). Because the rod path is circular, it is more convenient to use polar coordinate than Cartesian coordinate. From the scanner geometry, X p (k,i) can be expressed as r $: 2 V,0 = 2 J 0i ^( r rod cos 9 > r rod sin 9 ) P (*> *l r rod cos r rod sin d ) r iod d9 Â’ (2-40) where A(r r0( j cos 9, r roc j sin 9) is the mean number of photon pairs that are emitted from the rod at the angle 9, r ro( j is the radius of the path of the rotating rod, and P(k,l\r T0 ^ cos 9, r roc jsin0) is the probability that a photon pair emitted from the rod at the angle 9 is incident upon the detector pair (&,/). The factor of two accounts for the fact that the rod passes through the same detector tube twice in one cycle. We assume that the rod rotates at a constant angular velocity and emits positron uniformly. Under these assumptions, it follows that A(r ro{ jcos0, r^sinf?) is constant over 9. Let A = 2A(r rO( jcos0, r r0( j sin#)r r0( j. Then, the integral in (2.40) may be approximated by the following sum (see Fig. 2.7): A p (fc,() = \P(nA9)A9 = \c(k, /), (2.41) n where c(M) = E p nA0, (2.42) n and P n = P{nA9) is the probability that a photon pair emitted from the rod at the n th position is incident upon the detector pair (k,l). This probability may be calculated using the angle-of-vies (AOV) method [14,36]. Specifically, the probability that a photon pair emitted from a point (x, y) is incident on detector pair (k, l ) is determined by the angles a u a 2 , a 3 , and q 4 (see Fig. 2.8). The AOV is the maximum

PAGE 39

32 angle for which a photon pair would be incident on a detector pair and is given by Then, the probability that a photon pair emitted from the point ( x , y ) is incident on detector pair ( k,l ) is equal to AOV/7T [36]. Note, the parameter c(k,l) depends Figure 2.7: A proposed way to calculate c(k,l) to approximate A p(fc)/ ). only on p(k,l). The parameters {c(k,l)} for the Siemens-CTI ECAT EXACT 921 scanner are plotted in Fig. 2.9 as a function of (/ k). Equipped with the above information, we assume that N kt and B kt are Poisson with means Xc(k,l) and e k i\c(k,l), respectively: (2.43) N k i ~ Poisson(Ac(A;,/)), k = 1,2, e fan fc , B ki ~ Poisson(e fc iAc(M)), k = 1,2 , e fan*. (2.44) (2.45)

PAGE 40

33 Figure 2.8: The four angles to determine the angle-of-view associated with the point ( x,y ) and detector pair (k,l). We refer to the data model described in equations (2.44) and (2.45) as Poisson Model II. 2.2.1 Estimation Problem Associated with Poisson Model II As discussed in the previous section, the parameters {c(A;, /)} are determined from the geometry of the scanner. The unknowns X dl ,X d2 ,..., X do for Poisson Model I have been reduced to only one unknown, A, for Poisson Model II. The estimation problem now becomes (P2) (e, A) = argmax 2 (€, A) subject to 0 < e < 1 and A > 0, (2.46)

PAGE 41

34 c(k,l) Figure 2.9: The parameters {c(k, /)} for the Siemens-CTI ECAT EXACT 921 scanner calculated according to (2.42). where 0 2 (e> A) is the log-likelihood function: 2 (e, A) is not concave. Unfortunately, although the number of unknowns is reduced in (P2), the Fisher information matrix for the probability density function associated with (P2) is still non-invertible (see Appendix B). 2.2.2 Estimation Algorithms Associated with Poisson Model II We show in Appendix C that the conditional expectation of the log-likelihood function of the complete data X = [B,N] given the incomplete data and the current

PAGE 42

35 estimates e m and A m is U'{e, A; € m , A m ) = E[L' C (B, N\e, X) \ B = 6; e m , A m ] / = S H [5jt< log(e fc e,) + (nJ5 b kl ) log(l e fc e,)] k=1 Jefan k + H 5Z [-Ac(M)+n51ogAc(M)] + C 2 , (2-48) fc=1 iefan fc where L' C (B, iV; e, A) is the log-likelihood function of the complete data: L' C (B, N ; e, A) = P(B = 6, TV = n; e, A), (2.49) and C 2 is independent of e and A. We also show in Appendix C that, for k = 1,2,...,/,/ £ fan fe , nS = ^[Mh I B kl = b kl e m , A m ] = (1 e^)X m c(k, l ) + 6 W . (2.50) The function £/'(e, A; e m , A m ) is given by U'(e, A; e m , A”) = 9 <“>(e) + ft’W + C 2 , (2.51) where 5 M( e ) _ ^ [6 W log(e fc C|) + (n)J| 6 W ) log(l etc,)] (2.52) fc=1 Jefan * and / 2 (m) (A)=^ 2 [~MM) +n2}l°gAc(fc,0]. (2.53) fc=1 (efan it

PAGE 43

36 The steps of the EM algorithm are: Step I. ( Initialization ) Set m=0. Choose e° and X° s.t. 0 < e° < 1 and X° > 0. Step II. ( E step) nft = (1 “ e?e?)X m c(k, l ) + b kl , k = 1, 2, U (Pe) e m+1 = arg max g^ie). v ' °o 0 momentarily. Taking the derivative of X ) with respect to A and setting it to zero, we get I qr\ CH E E [-c(M) + -f ] = 0. (2.54) fc=1 iefan k Solving (2.54), the new estimate for A is A m+ 1 £ £ fe=1 iefan k E E c(k,l) k=1 /efan it (2.55)

PAGE 44

37 We refer to the estimation algorithm using (2.18) and (2.55) for Poisson Model II as the EM-CD2 algorithm. We refer to the estimation algorithm using (2.27) and (2.55) for Poisson Model II as the EM-FP2 algorithm.

PAGE 45

CHAPTER 3 ML ESTIMATION ALGORITHMS USING A REFINED DATA MODEL 3.1 Poisson Model III To get a more complete data model for the blank scan data, the effect of detector penetration must be incorporated. In Fig. 3.1, although the photons are incident upon detector pair ( k',l '), the photons may penetrate the detectors. Suppose the photon incident on detector k' penetrates detector k' and is detected by detector k, and the photon incident on detector l' penetrates detector l' and is detected by detector l. Then, a coincidence may be recorded erroneously by detector pair ( k,l ). It should be noted that a coincidence may also be erroneously recorded by detector 38

PAGE 46

39 pairs (k, l') or (k 1 , l) when detector penetration takes place, even though the photons are incident upon detector pair ( k',l '). Even if the photons stop within detectors k and l, the photons might not be detected by detectors k and l. For a photon to be detected by a detector, the photon must have enough energy after stopping to interact with the detector [25]. The probability that the photon penetrates detector k' (or survival rate) is [38,39], where L k > is the length of intersection of detector k! and the path of the photon (see Fig. 3.2) and p, is the attenuation coefficient of the detector. Therefore, for the given path of the photon pair, the probability that the photon incident on detector k! is detected by detector k is equal to p k e~ L ^{\ e~ L ^), (3.1) where p k is the probability that a photon event is recorded given that a photon has been stopped by detector k, and (1 e ~ Lktl ) is the probability that the photon is stopped by detector k. The attenuation coefficient of detectors is assumed to be known. For BGO detectors, the linear attenuation coefficient is 0.96 cm -1 [25] and the intersection lengths can be calculated from the scanner geometry. Although we

PAGE 47

40 focused on one photon of the photon pair, similar statements can be made about the other photon. When the effect of detector penetration is considered, the tube associated with a detector pair is defined differently (see Fig. 3.3). The mean number of photon pairs Figure 3.3: The tube associated with detectors k and l is defined by the dashed lines when detector penetration is taken into account. Note that, to clearly depict the angles, detectors k and l are magnified. that are detected by detector pair ( k , l) is given by E[B kl } = I \(x,y)p k piP(k,l\x,y)dxdy ( 3 . 2 ) rod path = 2 if A ( r rod cos 9 Â’ r rod sin 9 )pkPiP{k, l\r vod cos d, r rod sin 0)r vod d9 , where A(r rod cos0, r rod sin0) is the mean number of photon pairs that are emitted from the rod at angle 6, and P(k,l\r Tod cos 9, r vod sind) is the probability that a

PAGE 48

41 photon pair emitted from the rod at angle 6 is stopped within detector k and /, respectively. Since many paths are possible for a photon pair that is incident on detector pair (k,l), the probability P(k, l\r YQ( ^cosd, r ro( jsin0) can be expressed as p ( k > l \ r rod cos r rod sin e ) — e -L k i(a)n^ _ e -L fc (a)/j^ e -L ( ,(7r+a)/x^ _ g-L; (7r+a)/j) da where Lk(a) is the intersection length of the path of the photon traveling at angle a and detector k, and Lk>(ai) is the intersection length of the path of the photon traveling at angle a and detector k' (see Figs. 3.4 and 3.5). The intersection lengths Li(tt + a) and L;/(7 t 4a) are defined similarly. Note that a photon might have k AOV Figure 3.4: Angles that are used to calculated P(k,l\r TQ ^ cos 9, r rO( jsin0) in (3.3). penetrated several detectors before reaching detector k. In this case, (3.3) would have to be modified accordingly.

PAGE 49

42 k AOV Figure 3.5: Angles that are used to calculated P(k,l\r rod cos 0,r rod sin 9) in (3.3). For convenience, we let A = 2A(r rod cos0, r rod sin0)r rod . In practice, the exact integral in (3.2) cannot be determined. Therefore, we approximate the integral by N e A pkpi P(nA9)A6, n= 1 (3.4) where N a P{nA6) — '^2e~ Lk ' ( ' iAa \l — e ~ Lk ^ iAo ‘^)e~ Ll ^' n+iAa \l — c L ibr+»Aot)) ^ 5 ) i — I 7T i= 1 AO = (0 2 0, )/Avy, (3.6) and An = AOV/N a . (3.7) The quantity in (3.4) can be expressed as A p k pia(k,l), (3.8)

PAGE 50

43 where N e a(k,l ) = ^P(nA0)A0. 71=1 (3.9) Observe that a(k, l) depends only on p{k, l). The parameters {a(k, /)} for the SiemensCTI ECAT EXACT 921 scanner are plotted in Fig. 3.6. The number of photon pairs a(k.l) x 10Figure 3.6: The parameters {a(k, l)} for the Siemens-CTI ECAT EXACT 921 scanner calculated according to (3.9). that are detected by detector pair ( k , l ) is therefore assumed to be an observation of the random variable B kt , where Bid ~ Poisson(A p k p t a(k, /)), k = 1,2, . . . ,1, l e fan fc . (3.10) In order to utilize the same methodology for estimating unknown parameters of Poisson Model II, the means of random variables {Pit;} must be related to the means of some other Poisson random variables in the same manners that the means of {B kl } are related to the means of {N kt }. Let N kt be the number of photon pairs that are

PAGE 51

44 stopped within detectors k and l. Then, the random variables {-/V^} have means {A a(k, l )}: N ki ~ Poisson (Aa( Ac, /)), k = 1, 2, . . . , J, l e fan*,. (3.11) We refer to the data model described in equations (3.10) and (3.11) as Poisson Model III. 3.1.1 Estimation Problem Associated with Poisson Model III Since the parameters {a(k,l)} can be determined from the scanner geometry (see (3.5) and (3.9)), the estimation problem now becomes (P3) (p, A) = argmax0 3 (p, A) subject to 0 < p < 1 and A > 0, (3.12) where p is the vector with elements {p k } and
PAGE 52

45 estimates p m and A m is U"(p, A; p m , A m ) = E[L" C {B : N-p,X)\B = b, p m , X m ] / = E E i b ki MPfcP/) + hi) log(l PkPl )\ k ~ l ie fan k i + [~Aa(M) + <} log Aa(A;, /)] + C 3 , ( 3 . 14 ) * =1 iefan fc where L"(B, IV; p, A) is the log-likelihood function of the complete data: L"{B, N; p,X) = log P{B = b, N = n, p, A), (3.15) and where C 3 is independent of p and A, We also show in Appendix C that, for k = 1 , 2 , e fan fc , I = 6 W ; p m , A™] = (1 p-p-)A m a(A:, /) + b kl , (3.16) The function U"(p, A; p m , A m ) is decoupled and given by u"(p, A”*) = gW(p) + /<”»(A) + C 2 , (3.17) where 9 [m \p) = £ £ log(PfcPi) + (nJ3 M log(l p kPl )] (3.18) k=:1 iefan k and /s"”(A)= £ £ (-Ao(*,i) + nSlogAo(t,l)]. d=di {{k,l):p(k,l)=d} (3.19)

PAGE 53

46 The steps of the EM algorithm are: Step I. ( Initialization ) Set m=0. Choose p° and A 0 s.t. 0 < p° < 1 and A 0 > 0. Step II. ( E step) n ki = { 1 Pk ' 0 + bki, k = 1, 2, . . . , I, l e fan k . Step III. (M step) (P A ) A m+1 = argrnax/ 3 (m) (A). A>o ( p p) P m+1 = arg Q ma x ^ m )(p). Step IV. If stopping criterion is not satisfied, set m = m + 1 and go to Step II. Otherwise stop. Repeating the steps in Section 2.2.2, the problem (Pp) is solved either by using a CD-based method Pk = arg max 0
PAGE 54

47 We refer to the estimation algorithm using (3.20) and (3.22) for Poisson Model III as the EM-CD3 algorithm. We refer to the estimation algorithm using (3.21) and (3.22) for Poisson Model III as the EM-FP3 algorithm.

PAGE 55

CHAPTER 4 SIMULATION In this chapter, we assess the performance of the proposed algorithms using synthetic PET data and data from a thorax phantom by Data Spectrum Inc. For the purpose of making comparisons, we apply the fan-sum algorithm and Ferreira’s algorithm to this data as well. Note although Ferreira’s algorithm was developed for a cylindrical source, the underlying assumptions could also be made for rotating rod sources. For all the iterative algorithms, the initial estimates for the detector efficiencies are e° k = 0.5 for all k. Given this choice, the initial estimate for A for the EM-CD and EM-FP algorithms is A 0 A d — E {(k,iy. P (k,i)=d} number of pairs (k, l ) s.t. p(k, l) = d , d — di, d 2 , • • • , dp, (4.1) the initial estimate for A for the EM-CD2 and EM-FP2 algorithms is E E _ k=1 zefan^ k * E E c(k,l) k ~ 1 /efan fc (4.2) and the initial estimate for A for the EM-CD3 and EM-FP3 algorithms is A 0 = E E fc -Uefan t k k E E a(k,l) k ~ l zefan* (4.3) 48

PAGE 56

49 The steps of the EM-CD and EM-FP algorithms are repeated until the following stopping criterion is met. £( : fc=i ,m+ 1 — e k \2 dD \m+l \m + £e L ^P L ) 2 <<*, d=d i \m A d (4.4) where 8 is a small positive number. In our simulations, the EM-FP algorithm converged much faster than the EM-CD algorithm. Consequently, we used S e = 8 = 10 -11 for the EM-FP algorithm and 8 e = 8 = 10 -3 for the EM-CD algorithm. As for both the EM-CD2 and EM-FP2 algorithms, the stopping criterion is £( ! A:=l m+1 _ ,m \m+l _ \m e k \2 /T A \2 n / ' V A m f<6, (4.5) where £ is a small positive number. In our experiments, 6 e = 8 = 10 -3-5 for the EM-FP2 algorithm and 8 t = 5 = 10 -2 5 for the EM-CD2 algorithm. By contrast, Ferreira’s algorithm did not converge even after 300,000 iterations when the value 8 e = 10 -1 was used. Therefore, in our simulations Ferreira’s algorithm was executed for 250 iterations as specified by Ferreira et al. [17]. 4.1 Synthetic Data Using Poisson Model II, synthetic data are created to simulate blank scan data obtained from the Siemens-CTI EC AT EXACT 921 scanner. The scanner has 384 detectors per ring, and the ring diameter is 84.5 cm. The region to be reconstructed is a 43.9x43.9 cm 2 square consisting of (128) 2 voxels. Each plane consists of 192 projections and 160 members per projection. More information of this scanner can be found in Wienhard et al. [26]. As discussed in Chapter 2, this scanner has 81 possible perpendicular distances. The synthetic data consist of observations of Poisson random

PAGE 57

50 variables { B ki }: B ki ~ Poisson(e fc e ( Ac(A:, /)), k = 1,2,..., 384, l £ fan fc , (4.6) where A = 69000, (4.7) and {c(k,l)} as discussed in Chapter 2. The value of A in (4.7) is chosen so that the synthetic blank scan data corresponds to 2-hour scan. Uniform, piece-wise uniform, and nonuniform detector efficiencies are used in our experiments presented. As customary [17], the detector efficiency estimates are normalized so that for each estimation algorithm the mean of the detector efficiency estimates is equal to one. Because the EM-CD3 and EM-FP3 algorithms are based on Poisson Model III, which incorporates the effect of detector penetration, the EM-CD3 and EM-FP3 algorithms are not included. The reason is because synthetic data do not account for the effect of detector penetration. On the other hand, if the synthetic data did incorporate the effect of detector penetration, it would be unfair to compare the EMCD3 and EM-FP3 algorithms against estimation algorithms that do not take detector penetration into consideration. 4.1.1 Comparison of Detector Efficiency Estimates For our first experiment, the detector efficiencies were uniform: The true efficiencies are shown in Fig. 4.1(a), and the estimates provided by the fan-sum, Ferreira’s, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are plotted Cfc = 0.8, k = 1,2,..., 384. (4.8)

PAGE 58

51 in Figs. 4.1(b)-4.1(g), respectively. From the estimation results, we conclude for all four of the proposed algorithms that the detector efficiency estimates have smaller variance. Note that the fan-sum estimates vary over a relatively wider range and the estimates provided by Ferreira’s algorithm show a pronounced sinusoidal pattern. For our second experiment, the detector efficiencies were piece-wise uniform: £k 0.8, k = 1,2,..., 192 0.4, k = 193,194,392 (4.9) The true efficiencies are shown in Fig. 4.2(a), and the estimates provided by the fan-sum, Ferreira’s, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are plotted in Figs. 4.2(b)-4.2(g), respectively. The fan-sum algorithm provides very poor estimates, while the estimates provided by Ferreira’s algorithm again exhibits a sinusoidal pattern. From the experiments, the following conclusions are drawn about the proposed algorithms. The efficiency estimates for the first half of the detectors show little deviation from a constant value. This is also true for the second half of detectors. In addition, the ratio of the efficiency estimates for the first half of detectors to the efficiency estimates for the second half of detectors is very close to the true ratio 0 . 8 / 0 . 4 . To assess performance of estimation algorithms, we compare the mean-square error (MSE) of the efficiency estimates of 50 Monte Carlo (MC) runs. The error criterion MSE is defined as I i 50 MSE = £™Efe«-c*) 2 , (4.io) fc=i ou i=i where c k (i) is the estimate for e k in the i th MC run. The results are tabulated in Tables 4.1 and 4.2, respectively. Although each of four of the proposed algorithms provides detector efficiency estimates with smaller MSE than the fan-sum and Fer-

PAGE 59

52 reiraÂ’s algorithms, the difference between all algorithms is small. Further comparison is made based on a different criterion in the following. Table 4.1: MSE of the detector efficiency estimates using synthetic data and uniform detector efficiencies in the first experiment. Algorithm Fan-sum FerreiraÂ’s EM-CD EM-FP EM-CD2 EM-FP2 MSE 15.3721 15.4538 15.3623 15.3623 15.3630 15.3630 Table 4.2: MSE of the detector efficiency estimates using synthetic data and piecewise uniform detector efficiencies in the second experiment. Algorithm Fan-sum FerreiraÂ’s EM-CD EM-FP EM-CD2 EM-FP2 MSE 69.3189 69.5227 68.0414 68.2702 68.2334 68.1167 Ideally, the detector efficiency estimates would be such that the ratio i k /e k is approximately equal to a constant for all k. Therefore, as a merit of figure, we use the variance of the ratio of the estimates and true values: VR= sample variance of {ei/ei,e 2 /e 2 , ,e//e/}. (4.11) The idea is that the smaller the VR the better. To further compare the algorithms, 50 MC runs were performed for each of the first and second experiments and the VR was computed for each run. The plots of the VR versus MC runs are shown in Figs. 4.3 and 4.4 for the first and second experiments, respectively. From the plots, it can be clearly seen that all four of the proposed algorithms have comparable performance and each of them have much smaller VR than the fan-sum algorithm and FerreiraÂ’s algorithm.

PAGE 60

53 For our third experiment, we used nonuniform detector efficiencies that were generated as follows: & if 0 < £* < 1 i if e* > i k = 1,2,..., 384, (4.12) where & = 0.5 + V0.008Af(0, 1), k — 1,2, ...,384. (4.13) Obviously, only positive values are used for {£ fe }. The parameters in (4.13) are chosen so that with high probability the detector efficiencies lie in the interval [0.3, 0.7], which is a typical range for detector efficiencies. The ratio of detector efficiency estimates to true detector efficiencies is plotted in Figs. 4.5(a)-4.5(f), where the detector efficiency estimates are provided by the fansum, Ferreira’s, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms, respectively. Based on Fig. 4.5, we confer that all four of the proposed algorithms provide much better detector efficiency estimates than the fan-sum and Ferreira’s algorithms. Tabulated in Table 4.3 is the MSE of the detector efficiency estimates generated by using each of the estimation algorithms. All four of the proposed algorithms provide detecTable 4.3: MSE of the detector efficiency estimates using synthetic data and nonuniform detector efficiencies in the third experiment. Algorithm Fan-sum Ferreira’s EM-CD EM-FP EM-CD2 EM-FP2 MSE 100.8640 100.7207 100.2472 100.2704 100.2215 100.2070 tor efficiency estimates with smaller MSE than the fan-sum and Ferreira’s algorithms. We also compute the VR in each of 50 MC runs. The results are plotted in Fig. 4.6. It can be seen that all four of the proposed algorithms have comparable performance, and each of them is clearly better than the fan-sum algorithm and Ferreira’s algo-

PAGE 61

54 rithm in the sense that they provide estimates with smaller VR in every MC run than the fan-sum and FerreiraÂ’s algorithms. Another comparison criterion adopted by other researchers is the relative error [9,17]: The maximum relative errors are plotted in Fig. 4.7, and the average relative errors are plotted separately in Fig. 4.8. Based on the maximum relative error, all four of the proposed algorithms have comparable performance and each of them is better than the fan-sum and FerreiraÂ’s algorithms. In terms of mean relative error, the EM-FP algorithm has the best performance, and both the EM-CD and EM-CD2 algorithms are better than the fan-sum and FerreiraÂ’s algorithms in every MC run. To remove the bias factor from the comparisons, we normalize the detector efficiency estimates so that their mean is equal to the mean of the true detector efficiencies and then compute the relative errors. Note that in practice this cannot be done because the true detector efficiencies are not known. The maximum and average relative errors associated with the normalized detector efficiency estimates are plotted in Figs. 4.9 and 4.10. The plots confirm that all four of the proposed algorithms have comparable performance, and each of them outperform the fan-sum and FerreiraÂ’s algorithms. Specifically, the proposed algorithms provide estimates with smaller maximum and average absolute errors than the fan-sum algorithm and FerreiraÂ’s algorithm in every MC run. 4.1.2 Comparison of Reconstructed Emission Images The efficiency estimates were then used to reconstruct a synthetic emission image. Assuming only errors due to detector inefficiency, the emission data are assumed to be a realization of the set of Poisson random variables {!*/} given in (1.3). For

PAGE 62

55 convenience, the data model for {Yu} is expressed again. The emission data, {yu}, are modeled as observations of the independent random variables {}*/}: Y ki ~ PoissonfoejPx]*/), k = 1,2,..., 384, l e fan*, where l Px U = ± p ujXj, j 1 (4.15) and Xj is the mean number of photon pairs emitted from voxel j and J is the number of voxels. In addition, Pu,j is the probability that a photon pair emitted from voxel j is incident on detector pair ( k,l ) and x = [xi,x 2 ,...,xj] T is referred to as the emission intensity vector. Given the observed emission data y, which is a vector with elements {yu}, the problem is to estimate the emission intensity vector x. Given the synthetic emission intensity x in Fig. 4.11, we use the EMML algorithm [1,14,40,41] to estimate x from y: a(«+i) _ Xj ~( n + J ) ti £ P kl,j k,l E PgjVu [P c i ( “V (4.16) where P kl,j ~ e k^l P kl,j(4.17) The efficiency estimates were used in (4.17) to reconstruct the emission image. In Figs. 4.12(a)-4.12(d) are the plots of the reconstructed emission images at the 15 th iteration of the EMML algorithm using true detector efficiencies, no normalization, fan-sum algorithm, and Ferreira’s algorithm, respectively. (The EMML algorithm virtually converges after 15 iterations.) And in Figs. 4.13(a)-4. 13(d) are the plots of the reconstructed emission images at the 15
PAGE 63

56 defined for each iteration of the EMML algorithm as J i 50 MSE = £^E(*;»-*i) 2 , (4.18) j=l 0U i=l where xj(i) is the estimate for xj at the i th MC run. The curves of the MSE of reconstructed images versus iteration number of the EMML algorithm are shown in Fig. 4.14. Because the mean of true detector efficiencies is different from the mean of estimates of each estimation algorithm, bias becomes an issue and the MSE corresponding to true detector efficiencies is not included in the plot. While all Table 4.4: MSE of the reconstructed emission images from synthetic data. Iteration of EMML 16 17 18 19 20 Fan-sum 786.1761 785.5547 788.4473 794.3706 802.9343 Ferreira’s 786.1472 785.5140 788.3954 794.3079 802.8612 EM-CD 786.1328 785.4880 788.3587 794.2610 802.8044 EM-FP 786.1363 785.4823 788.3440 794.2375 802.7723 EM-CD2 786.2022 785.3626 788.0479 793.7720 802.1426 EM-FP2 786.1835 785.2883 787.9220 793.5967 801.9190 the estimation algorithms improve the reconstructed emission images over using no normalization, the curves for all the estimation algorithms are close together. The values of MSE are tabulated in Table 4.4 for iterations 16-20 of the EMML algorithm for clarity. 4.2 Real Data Real PET data were obtained from Dr. John Votaw of the PET Laboratory at the Emory University School of Medicine in Atlanta, GA. The PET scanner is the Siemens-CTI ECAT EXACT 921 model. The collected data consist of 47 planes, where each plane consists of 192 projections and 160 members per projection. The

PAGE 64

57 scanner has 384 detectors per ring and the region to be reconstructed is a 43.9x43.9 cm 2 square consisting of (128) 2 voxels. The fan-sum method, FerreiraÂ’s method, and the proposed algorithms were used to estimate the detector efficiencies for plane 10. The duration of the blank scan was 30 minutes. Since the true detector efficiencies are unknown, we have no objective way to compare the detector efficiency estimates. The efficiency estimates were used in conjunction with the EMML algorithm to reconstruct the emission image for plane 10. The duration of the emission scan was 7 minutes. In Figs. 4. 15(a)-(d) and 4.16 (a)-(b), the reconstructed emission images at the 20 th iteration of the EMML algorithm using the fan-sum, FerreiraÂ’s, EM-CD, EM-FP, EM-CD2, EM-FP2 algorithms, respectively, are shown. The reconstructed emission images are indistinguishable. As the EM-CD3 and EM-FP3 algorithms are based on Poisson Model III which incorporates the effect of detector penetration, when applying the EM-CD3 and EM-FP3 algorithms to reconstruct emission images, the probability matrix P in (1.3) must be defined differently. Let P k ij stand for the probability that a photon pair is stopped by detectors k and l given that the photon pair is emitted from voxel j. The methods by Shepp et al. [14] and Carroll [36] can be extended to compute these probabilities [42]. The reconstructed emission images using the EM-CD3, EM-FP3 algorithms and no normalization (i.e., {p k = 1}) are plotted in Figs 4.17(a)-(c). The reconstructed emission images for the EM-CD3 and EM-FP3 algorithms show details (e.g. heart, lungs, and spine) much more clearly than the reconstructed emission images for the other estimation algorithms. Nevertheless, it is difficult to objectively determine whether using the EM-CD3 and EM-FP3 algorithms lead to better reconstructed emission images than using no normalization. To date, there is no accepted methodology for assessing the performance of detector efficiency estimation algorithms when real data are used. Dr. Votaw suggested the following experiment. A cylindrical phantom with uniform activity could be scanned

PAGE 65

58 for an extremely long time ( > 48 hours). Then, the emission images corresponding to each detector efficiency estimation algorithm would be reconstructed. Finally, the average variance of the pixel values for each image would be computed. The detector efficiency estimation algorithm that gave rise to the smallest variance would be considered to be the best algorithm. The shortcoming of the experiment described above is that only a single realization of the blank scan data is used. To account for the statistical properties of the detector efficiency estimation algorithms, multiple realizations of the blank scan data would be needed. Unfortunately, it is difficult to obtain multiple realizations because of the extremely high scanning costs.

PAGE 66

59 (a) i c\a (b) (C) (d) 1 .U4 1.04 1.02 1.02 Jlllt 1.02 1 yiUL 1 A., 1 0.98 n qa 0.98 V 0.98 100 200 300 100 200 300 u.yo 100 200 300 0.96 100 200 300 (e) (f) (9) 1 .04 1.02 1.02 1 1 0.98 0.96 0.98 100 200 300 100 200 300 U.yo 100 200 300 Figure 4.1: True efficiencies of synthetic data in the first experiment and the estimates: (a) true efficiencies; the others are the detector efficiency estimates generated by using (b) fan-sum algorithm; (c) FerreiraÂ’s algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.

PAGE 67

60 1 0.8 0.6 0.4 0.2 ( a ) (b) (c) (d) 100 200 300 100 200 300 100 200 300 100 200 300 (e) (0 (g) 1.6 1.4 1.2 1 0.8 0.6 0.4 1.6 1.4 1.2 1 0.8 0.6 0.4 1.6 1.4 1.2 1 0.8 0.6 0.4 100 200 300 100 200 300 100 200 300 Figure 4.2: True efficiencies of synthetic data in the second experiment and the estimates: (a) true efficiencies; the others are the detector efficiency estimates generated by using (b) fan-sum algorithm; (c) FerreiraÂ’s algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.

PAGE 68

61 x 10 i r 5 — (a) Fan-sum algorithm X (b) Ferreira’s algorithm (c) EM-CD o (d) EM-FP (e) EM-CD2 A (f) EM-FP2 4|x x x x X X X X xxxx „ x „xx X X X x X x x X X * X v x X X"x x „" x ''xx x "xxx A "XX x x x x x X >3 2 1 10 15 20 25 30 MC run number 35 40 45 50 Figure 4.3: VR of the detector efficiency estimates in the first experiment versus MC run number using synthetic data, uniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 69

62 o 1 1 1 1 — (a) Fan-sum algorithm _ X (b) Ferreira’s algorithm (c) EM-CD o (d) EM-FP (e) EM-CD2 A (f) EM-FP2 -2 . cc -4 > D) O -6 xxxxxxxxxxxxxxxxxxxxxxxxxxxx x xxxxxxxxxxxxxxxx x xx -8 -12 1 L 10 15 20 25 30 MC run number 35 40 45 50 Figure 4.4: VR of the detector efficiency estimates in the second experiment versus MC run number using synthetic data, piece-wise uniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 70

63 (a) (b) (c) (d) 2.05 2.05 1.95 100 200 300 1.95 100 200 300 Figure 4.5: The ratio of the detector efficiency estimates and true detector efficiencies, where the detector efficiency estimates are generated using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) FerreiraÂ’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 71

64 MC run number Figure 4.6: VR of the detector efficiency estimates in the third experiment versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) FerreiraÂ’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 72

65 i.ii i.i 1.09 1.08 "i r ~i r — (a) Fan-sum algorithm X (b) Ferreira’s algorithm (c) EM-CD o (d) EM-FP (e) EM-CD2 A (f) EM-FP2 o k_ CD CD .a 1.07 j5 at x x x x X XX X X X X X * X X x x XX x X x X x X X X X X v X X x X X E D E a I 1.06 / \ r . /> /\ __ ' / " N / N 7 ' '\ /\ r \ '\ ^ . J. 1.05 / N ^ \ ^ " s “ 1 \ E __ A JflA o [1 a P O o 5? °i < V? ajjo 20 25 30 MC run number Figure 4.7: Maximum relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 73

66 1.0161 1.016 1.0159 o 1.0158 <5 1 r ”i r — (a) Fan-sum algorithm X (b) Ferreira’s algorithm (c) EM-CD o (d) EM-FP (e) EM-CD2 A (f) EM-FP2 h I \ I \ A / t \ \ , . , /\ , A 1 ' ' / ' 1 / V/N/'A ' / » / ' x' J A V '/ / '/ f\ K \ ' A ' ' 7 x \/ _ I A / ' / \ / * 20 25 30 MC run number Figure 4.8: Average relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) Ferreira’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm.

PAGE 74

67 Figure 4.9: Maximum relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) FerreiraÂ’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm. Note that the detector efficiency estimates are normalized so that their mean equals the mean of true efficiencies.

PAGE 75

68 Figure 4.10: Average relative error of the detector efficiency estimates versus MC run number using synthetic data, nonuniform detector efficiencies, and: (a) fan-sum algorithm; (b) FerreiraÂ’s algorithm; (c) EM-CD algorithm; (d) EM-FP algorithm; (e) EM-CD2 algorithm; and (f) EM-FP2 algorithm. Note that the detector efficiency estimates are normalized so that their mean equals the mean of true efficiencies.

PAGE 76

69 Figure 4.11: Synthetic thorax phantom.

PAGE 77

70 (a) (b) Figure 4.12: Reconstructed images of the thorax phantom in Fig. 4.11 at the 15
PAGE 78

71 (a) (b) Figure 4.13: Reconstructed images of the thorax phantom in Fig. 4.11 at the 15 th iteration of the EMML algorithm using: (a) EM-CD algorithm; (b) EM-FP algorithm; (c) EM-CD2 algorithm; and (d) EM-FP2 algorithm.

PAGE 79

72 Figure 4.14: MSE of the reconstructed emission images versus the iteration number of the EMML algorithm using: (a) no normalization; (b) fan-sum algorithm; (c) FerreiraÂ’s algorithm; (d) EM-CD algorithm; (e) EM-FP algorithm; (f) EM-CD2 algorithm; and (g) EM-FP2 algorithm.

PAGE 80

Figure 4.15: Reconstructed emission images of plane 10 of a thorax phantom at the 20 th iteration of EMML algorithm using: (a) the fan-sum algorithm; (b) Ferreira’s algorithm; (c) EM-CD algorithm; and (d) EM-FP algorithm. (a) (b) >* lAr § . . § -v' -’J. > t * «ti *L.-. v. * * ;V f% ;u : $:?<* K H-f' Figure 4.16: Reconstructed emission images of plane 10 of a thorax phantom at the 20 th iteration of EMML algorithm using: (a) EM-CD2 algorithm; and (b) EM-FP2 algorithm.

PAGE 81

74 (c) Figure 4.17: Reconstructed emission images of plane 10 of a thorax phantom at the 20
PAGE 82

CHAPTER 5 CONCLUSIONS Algorithms for computing ML estimates of detector efficiencies in PET have been proposed. The algorithms follow from the EM algorithm and three Poisson models that were developed for the PET blank scan data. The refined Poisson model (i.e., Poisson Model III) directly accounts for the Poisson nature of blank scan data, certain geometric properties of scanners, and detector penetration. The close resemblance between real data and the synthetic data generated according to the proposed models provides some justification for the models. For each Poisson model, a closed form expression for the E step of the EM algorithm was derived. The corresponding M step was solved iteratively using the coordinate descent method and a fixed-point iteration. The detector efficiency estimation algorithms that were based on the coordinate descent method have the desirable property that the estimates were guaranteed to lie in the interval [0, 1]. The fixed-point iteration based algorithms do not have this property. However, they converge about 20 times faster than the coordinate descent based algorithms. All of the algorithms have the property that the corresponding log-likelihood function increase monotonically with increasing iterations. Experiments using synthetic and real data were used to assess the proposed algorithms. For the synthetic data case, the following objective criteria were used: (1) variance of the ratio of the efficiency estimates and true efficiencies, (2) maximum relative error, and (3) average relative error. In terms of these criteria, the proposed algorithms greatly outperformed the fan-sum and FerreiraÂ’s algorithms. However, the reconstructed emission images for all of the estimation algorithms were visually iden75

PAGE 83

76 tical. In addition, the mean-square error (MSE) values associated with the emission images were also similar, but the proposed algorithms did produce the smallest MSE. There is no accepted approach for comparing the performance of detector efficiency estimation algorithms when real data are used. We expect the EM-CD3 and EM-FP3 algorithms to have a large impact because Poisson Model III takes the effect of detector penetration into consideration. The reconstructed emission images for the fan-sum, FerreiraÂ’s, EM-CD, EM-FP, EM-CD2, and EM-FP2 algorithms are indistinguishable. When applying the EM-CD3 and EM-FP3 algorithms to reconstruct emission images, a different probability matrix P is used. The reconstructed emission images for the EM-CD3 and EM-FP3 algorithms show details much more clearly than the reconstructed emission images for the other estimation algorithms. While it is difficult to objectively determine whether using the EM-CD3 and EM-FP3 algorithms lead to better reconstructed emission images than using no normalization, we conclude that the effect of detector penetration is a significant factor when reconstructing emission images.

PAGE 84

CHAPTER 6 FUTURE WORK While the proposed algorithms converge in practice, the convergence properties of the proposed algorithms are not well understood. Theoretic analysis of the proposed algorithm will be pursued to determine their convergence properties. As we stated in previous chapters, the performance of the proposed algorithms depends on the initialization step because the objective functions are not concave. It is a common practice to use a constant value as the initial detector efficiency estimate for iterative algorithms (e.g., FerreiraÂ’s algorithm [17]). From our experiments, we found that the average values of the detector efficiency estimates obtained from the proposed algorithms are very close to the constant value that was used to initialize the algorithms. However, if the constant value is scaled, detector efficiency estimates are not correspondingly scaled. Therefore, the effect of the initialization step must be investigated. 77

PAGE 85

APPENDIX A CONCAVITY OF THE LOG-LIKELIHOOD FUNCTIONS A.l Poisson Model I The log-likelihood function for the blank scan data is given by (e, A) = log P[B = b\ e, A) = H + bki\og(e k ei\ P (k,i))] k iefan* = -Ib k i\og(e k eiXd)]. (A.l) d {(k,l):p(k,l)=d} To determine whether a function is concave or not, we check if its Hessian is negative definite. To determine if a hermitian matrix is negative definite, we use the following theorem [43]. Theorem 1 The following condition is necessary and sufficient for annxn hermitian H to be negative definite. The principal minors consisting of the the determinants of the k x k matrices in the top left-hand corner of H (k = 1, . . . , n) are all negative: Hn < 0 , H n H\ 2 H n H n h 12 <0, h 21 h 22 H 2 3 H 2 i h 22 H 3l h 32 h 33 The ( i,j ) element of Hessian for the log-likelihood function is Hij = d 2 M B = b; 0) ddfiOj (A.2) 78

PAGE 86

79 where 6 = [e, A]. The first partial derivatives of 0i(e, A) are 90i(e, A) de k X [-<^(*,1) + — ], k = 1, 2, . . . , I, iefan* ejfc and 901 (e, A) _ ^ 2^ l — e k£i + t-J) d — d\, d 2 , . {( k,l):p(k,l)=d } A d d\ d .,d DAnd, the second partial derivatives are 9 2 i(e,A) de k dei a 2 i(e,A) dXddtk d 2 i(e,A) d\dd\ m * E ^ if* = Z iefan* k -A p (jfc,/) if / € fan* ’ 0 otherwise E M(p(M)-<0], ie fan* E {{k,l):p(k,l)=d} [-^f] if d = m 0 otherwise (A.3) From (A.3), we know that 1 #11 — -~2 X ^li < 61 lefani #12 — #21 = 0 , and Therefore, #22 — — 2 XI ^21 < 0 . «efan 2 #11 #12 #21 #22 — # 11#22 — # 12#21 > 0 .

PAGE 87

80 So we conclude that H is not negative definite and that the objective function of (PI), 0i (e, A), is not concave. For the over-parameterized, approach, the same conclusion can be reached by a similar derivation. The log-likelihood function for the blank scan data is given by logP(B = b;e, A) = logU n e~ e *' A r(M) ( £ *< Wo)*** k zefan*. ^ kl ' = logn n d {(k,l):p(k,l)=d] hi= XI H Wfci'W.O + hi log(e k i\ p (k,i)) + log(6 w !)] k zefan*. = XI Y2 [—tki^d + hi log(ejkjAd) + log(6 fc /!)]. (A. 4) d {( k,iy. P (k,i)=d } The first partial derivatives of log-likelihood function are dlogP(B = 6;e,A) , , hi , , „ T , _ , jr = -'W-0 + — , k = 1. 2, • • • , I, le fan*,, (A.5) oe kl tkl and <91ogP(B = 6;e, A) dX d {( k,l):p{k,l)=d } 1 e ki + t^]> d — di, d 2 , . . . , do*d Further, the second partial derivatives are (A.6) d 2 log P(B = b; e, A) dt-kld^-mri = < if (k, l) = (m, n) 0 if (k, l ) 7 ^ (m, n) (A.7) d 2 log P(B = 6; e, A) dX dde kt = < — 1 if p(k, l) = d 0 if p(k, l) ^ d < 5 (A.8) d 2 log P(B = 6; e, A) dX d dX m = < E [-^] {(k,l):p(k,l)=d] ^ 0 if d = m if d ^ m (A.9)

PAGE 88

81 From (A.7)-(A.9), it can be seen that # n < 0, #12 = #21 = 0 , #22 < 0 . Therefore, So we conclude that the #11 #12 — # n #22 — # 12#21 > 0 . #21 #22 log-likelihood function is not concave. A. 2 Poisson Model II The log-likelihood function for the blank scan data is given by &(e,A) = logP(jB = 6; e, A) = E E [-e fc €|Ac(fc, I) + & fc/ log(e*e,Ac(A:,Z))]k iefan* The first partial derivatives of 0 2 (e, A) are £ [-6,Ac(t,J) + -],* = 1.2 /, iefan fe and d(j> 2 (e, A) dX = E E + * = 1 , 2 , fc=1 Jefan* (A-10)

PAGE 89

82 And, the second partial derivatives are d 2 4 > 2(6, A) dt k dei E ^ if k = l le fan*. * — Xc(k,l ) if l e fan*, ’ 0 otherwise d 2 <(> 2(6, A) d\de k d 2 Me,\) d\ 2 ~ E £ic{k,l), le fan*. ~ E E A 2 fc=1 iefan* (A.11) From (A. 11), we know that #11 — —“2 13 ^li < 0) Cl /efani #12 = #21 = 0 , and # 22 ~2 13 b 2 i < o. 2 iefan 2 Therefore, #11 #12 #21 #22 — #11 #22 — # 12#21 > 0 . So we conclude that if is not negative definite and that the objective function of (P2), 2 (e, A), is not concave. A. 3 Poisson Model III The objective function of the problem (P3) is A h (p, A) logP(B = 6; p, A)

PAGE 90

83 = E E [-pkPi^a(k,l) + b k i\og(p k pi\a(k,l))}. (A.12) k zefan*. Because Poisson Model III has the exact same functional form as Poisson Model II, the second derivatives of 3 (p, A) are a 2 2 ( P , A ) _ d \ 2 ~ E • if k = l zefan* — \a(k,l ) if l £ fan*; 0 otherwise E Pia{k,l), zefan fc -£ E fc=1 zefan*. (A.13) From (A.13), we know that #11 — — 2 EZ ^ < ^ zefani #12 = #21 = 0 , and #22 — — 2 EZ ^2Z ^ 0 zefan 2 Therefore, #11 #12 — # n #22 ~ # 12#21 > 0 . #21 #22 So we conclude that 7T is not negative definite and that the objective function of (P3), ^(p, A), is not concave.

PAGE 91

APPENDIX B CRAMER-RAO LOWER BOUNDS To compute the Cramer-Rao lower bound, the Fisher information matrix must be computed. In this appendix, we derive the Fisher information matrix for the sufficiently-parameterized approach associated with each Poisson data model. B.l Poisson Model I For Poisson Model I, the probability density function (PDF) is p ( b = 6; e, a) — n n k ~ l iefan k ,-e k e,\„ lk .n (^'W.O) 6 *' bki\ From (A. 3), we know that the second partial derivatives of logP(B = b; e, A) are d 2 log p(B=b,e,X) dt k dei d 2 log p(B=b,e, X) d\ddt k d 2 logP(ff = b;€,A) dX^dX m E f if A; = / /efan* k if l £ fan fc > 0 otherwise E [eiS(p(k,l)-d)], le fan* (B.l) E [~ a^] if d = m {( k,l):p(k,l)=d } 0 otherwise 84

PAGE 92

85 The expected value of bki is ekeiX p ^,i), so we can calculate the expected values of the second derivatives as follows. E tJ ^ £L if k = l r /Eictiik -A p (fc,i) if l G fan fc 0 otherwise g l aJ '° SF a Sr» 6i6 ' A> l6=g] = E {e,6(p(k,l)-d)], g [ gJg£ agr„ 6i6 ’ A) i6= B ] = /efan* {(fc,/):p(fc,0=d} 1 0 [->?] if d = m otherwise (B.2) Consider a simplified example, where there are only 4 detectors on the ring (Fig. B.l). Detector pairs (1,3) and (2,4) have the same perpendicular distance from the Figure B.l: A simplified example to derive the Fisher information matrix. origin and therefore have the same Ai as the mean of the transmitted photon count. Detector pairs (1,4), (2,1), (3,2), and (4,3) have the same perpendicular distance from the origin and therefore have the same A 2 as the mean of the transmitted photon count. The data {&13, &24) &14, &2i 5 &32, 643} are collected. The problem is to estimate

PAGE 93

86 { e i> ^2, £3) 64, Ai, A 2 }. Then the Fisher information matrix can be found to be A) = where and Tx A 2 Ai A 2 63 A 2 Ai A 2 63 e 2 + €4 T 2 A 2 Ai 64 ex + 63 A 2 T3 A 2 Ci e 2 + €4 Ai A 2 T4 e 2 64 + 63 64 61 62 T5 0 62 + 64 64 + 63 e 2 + 64 61+63 0 Tq 1 Ti — (6 2 A 2 + 63A1 + 64A 2 ), ^1 T 2 — (6iA 2 + 63 A 2 + 64Ai]), 62 T3 — (eiAj + e 2 A 2 + 64A 2 ), ^3 T4 = (eiA 2 + e 2 Ai + 63A 2 ), 64 T5 = T— ( e l e 3 + 6264), Te — t(6164 + e 2 ei + €362 + 6463). A 2 We proved that the determinant of I(e, A) is equal to zero. That is, 1(6, A) is noninvertible. The same conclusion can be extended to all higher-dimensional cases. B.2 Poisson Model II For Poisson Model II, the PDF is P(B = 6; €, A) = n e k=1 (efan k (e k e l \c(k,l)) bkl hi'.

PAGE 94

87 From (A. 11), we know the second partial derivatives of logP(JB = 6; e, A) are a 2 log P(B=b;C, A) de k dti a 2 log p(B=b-,e, A) a\de k a 2 \ogP(B=b-,e,\) a\ 2 E h it /efan* k —Xc(k, l ) if k — l if l G fan* 0 otherwise = E eic(k,l), ie fan*. = -E E iefan*. (B.3) The expected values of the second partial derivatives of logP(J3 = 6; e, A) are _ ^ ejMk& [ik = l Jefan fc k —Xc(k,l) if/ e fan* ’ 0 otherwise (B.4) S r° , '^ e,A) lb=fl] = E f|c(A iefan*. it. B ] = s e fc x iefan*. Consider the example in Fig. B.l again. There are two possible values for {c(k,l)}: ci for detector pairs (1,3) and (2,4), because these two detector pairs have the same perpendicular distance from the origin; c 2 for detector pairs (1,4), (2,1), (3, 2), and (4,3), because these four detector pairs have the same perpendicular distance from the origin. Therefore, the Fisher information matrix is E f ktP & b *' x> i 6=B ] = • Si Ac 2 Acx Ac 2 s 2 X c 2 s 3 Ac 2 Ac 2 S 4 Aci Ac 2 5 5 Ac 2 5 6 Ac 2 Aci Ac 2 S 7 s 2 s 4 s 6 5 8 s 9

PAGE 95

88 where 51 — (c3 c i + €202 + C4 C 2)> ci 5 2 = ^( e 3 c l + e 2C2 + QC 2 ), 53 = — (flC2 + £302 + £ 4 ^ 1 ), ^2 ^4 = — (fiC2 + £302 + £4Ci), S5 ~ (^lCi + £2^2 + £4^-2)) C3 Sq = ^( e l C l + e 2 c 2 + C4C2), S7 = — (^1^2 + £2 Cl + C3C2), £4 ‘S '8 = ^( e l c 2 C 2 Cl C 3 C 2 ), and Sg — — (£i£ 3 Ci + £2^1 + (.\t\C2 + £2^2 + £3^2 + £4^2)It can be readily proved that, regardless of the values Ci and c 2 , the fifth row of the matrix J(e, A) can be expressed as a linear combination of the first four rows. Therefore, the Fisher information matrix I(e, A) has a determinant equal to zero and is non-invertible. The same conclusion can be reached for the general case. B.3 Poisson Model III Because Poisson Model III has the exact same functional form as Poisson Model II, we have the same conclusion for Poisson Model III regarding the Fisher information matrix. In other words, the Fisher information matrix is non-invertible.

PAGE 96

APPENDIX C DERIVATION OF THE E STEP OF THE EM ALGORITHMS In this appendix, we prove that the E step is given by (2.9). Said another way, we determine the conditional expectation of the complete data, X = [B,IV], given the incomplete data, B, and current estimates of e and A. C.l Poisson Model I Using Bayes’ Rule [44], the log-likelihood function of the complete data is given by L c (b, ne, A) = log P{B = 6, N = n; e, A) = log P(B = b \ N = n; e, X)P(N = n; e, A) = lo S II II (?")(cfc€i) 6 ‘ , (l-efcei) Bw " 6w * =1 le fan* + log fi n d=d l {(k,l):p(k,l)=d) 71 «' / = [ lo S ("*/) + hi log(ejt Ci) + ( n k i hi)log(l e*ej)] k=1 zefan*. dp + H 5Z [-A d + n w log(A d ) + log(n fc/ !)]. (C.l) d=d, {{k,l):p(k,l)=d} 89

PAGE 97

90 Let U(e, A; e m , A m ) denote the conditional expectation of the log-likelihood function of the complete data given the incomplete data and estimates e m and A m . Then, U(e, A; e m , A m ) = E[L C {B, N; e, X) \ B = b-, e m , A m ], = E E (^[(S;)|S = 6;e' n ,A m ] + log(e,e / )^|B = 6;e-,A m ] /efan*. ^ _ + log(l-e k e l )E[N kl -B kl \B = b;e m ,\ m ]) + E E (~ + ^og(\ d )E[Nki \B = b\ € m , A m ] d=di {( k,t):p(k,l)=d } + B[log(iVfc/!) | B = 6;e m , A m ]). To maximize U(e, A; e m , X m ) with respect to e, we take into consideration only the terms log (e k e t )E[B k , \ B = 6; e m , A m ] and log(l e k e t )E[N kl B kl \ B = b; e m , A m ] as the remaining terms are independent of e. Similarly, we only need to consider the terms -X d and log(X d )E[N k i | B = 6; e m , A m ] in order to maximize U(e, A; e m , A m ) with respect to A. Let J = £[A^ | B = 6;e m , A m ], Because £[£„|B = 6;6™,A to ] = &„, + C/(e, A; e m , A m ) E E [hi log (cfcCi) + (n£} log(l e*e,)] fefan* do E E d=d! {(fc, /)*(*, <)=d} [-Ad + n$ log AJ + Ci, (C.3) where C\ is independent of e and A. To finish the E step, an expression for must be obtained.

PAGE 98

91 To find n™ t , we need to derive the conditional probability of N k i given the incomplete data. P(N kt = n k i\ B = 6; e, A) = P(N kl = n kl | B ki = b ki ; e, A) P(B k ,=b k ,\N k i=n k i)P(N k ,=n k i) P(B k ,=b k , ) n (C.4) e -‘* e l* p (*.0 < ~ eke ‘ x t>(.i' ; i) ) kl n l(l-tfce t )A p(ti , )]"«-**! ( n k i-b k i)\ The second equality follows from Bayes’ Rule. The conditional expectation of N k i is therefore E[N k i | B ki = b k i\ e, A] + + E n k i>b k i E n k i>b k i E n kl>bkl E n kl>bkl n k iP{N kl = n kl | B kt = b ki ; e, A) ( n ki ~ hi + b k i)P(N k i = n ki | B kl = b ki ; e, A) ( n ki ~ hi) P (N k i = n ki | B kt = b k i\ e, A) hiP(N kl = n k i | B kl = b k i ; c, A) E n *i— 0 p/ g— K* £ <= t d-^p(fc,/)] W Kd ! n' ^.n [( 1 ~t|fet[)Ap(jt i [)] n fcl Ki)* (C.5) where n' kl = n*,; b k i. The first summation is equivalent to the mean of a Poisson random variable with parameter (1 — e k i)A p ( k< ^, and the second summation is the sum of the probability over all possible values of n' kl . Therefore, E[N k i | B k i — b k i ; e, A] — (1 — 6 k Ci)A p ( k ,i) + b k i. (C.6) Given the previous estimates e m and A m , "2 = E[N, ki I B kl = b kl ; e m , A m ] = (1 + b, m, (C.7)

PAGE 99

92 for k = 1,2, e fan*. C.2 Poisson Model II With c(k, l ) known for k = 1, 2, . . . , I, l e fan*, the log-likelihood function of the complete data X = [B, AT] is given by L' c (b, n; e, A) = log P(B — b, N = n, e, A) = log P(B = 6 | iV = n; e, X)P(N = n; e, A) = lo sn II 0(^Ci) 6 *'(l-e it e i ) n ‘'6fc ' * =1 fefan* + log]] n e -Ac(fc, f )( Ac (fc.0) BM fc=1 /efan fc Ukl ' i = E E [ lo gC) + ^log(efcQ) + («w-^/)/o^(l-e jfc e,)] * =1 iefan*. i + 5Z H [-Ac(&,/) + n w log(Ac(M)) + log(n w !)]. (C.8) fc=1 fefan*. Let U'(e, A; e m , A m ) denote the conditional expectation of the log-likelihood function of the complete data given the incomplete data and estimates e m and A m . Then, U'{e, A; e m , A m ) = E[L' C (B , N ; e, \) | B = b\e m , A m j, = E E (£[(%;) I B = + log(t tE ,)£[B t « j B = fc 1 «efan fc + log(l — tkti)E[N k i — B k i | B = 6; c m , A m ]) + E E (— Ac(fc, /) + log(Ac(fc, l))E[N k i \ B = b] e m , A m ] fc 1 iefan/t + E[\og(N k i\) \ B = b: e m , A m ]). To maximize U'(e, A; e m , A m ) with respect to e , we take into consideration only the terms log(e*e,)£[B fc , | B = 6; e m , A m ] and log(l e fc e,)^[^w B ki \ B = 6; e m , A m ]

PAGE 100

93 as the remaining terms are independent of e. Similarly, we only need to consider the terms —A c(k, l ) and log(Ac(A:, l))E[N k i \ B = b\ e m , A m ] in order to maximize u '(e, A; e m , A m ) with respect to A. Let n” k } = E[N kl \ B = b, e m , A m ], Because E[B kl \ B — b\ e m , A m ] = b kt , U'(e, A; e m , A m ) = E E [bu log (cjfcCi) + {n% b k i) log(l e*e,)] fc=1 ie fan* + E E [~Xc(k,l) + rift log \c(k,l)] + C 2 , k=1 ie fan* (C.10) where C 2 is independent of e and A. To finish the E step, an expression for must be obtained. To find n£J, we need to derive the conditional probability of N kl given the incomplete data. P(N k i = n k i | B = 6; e, A) = P{N kl = n kl \ B ki = b ki ; c, A) — p { B kl=bki\N k i=n kl )P(N kl =n k i) P{B k i=b kl ) (”*' )(e fc c/) b *< (1 ekei )n k ,-b k , e —Xc(k,l) = — bi n kV e -e k ' t \c(k,l) (t k ti\c(k,l)) b kl b kl l „-(l-e k e,)\c(k.l) (nkl-bkl)'(C.ll)

PAGE 101

94 The second equality follows from Bayes’ Rule. The conditional expectation of N k i is therefore E[N kt | B kl — b k i ; e, A] = £ n kiP{N k i = n kl | B kl = b kl ; e, A) n kl>bkl + + £ («fei + b k i)P(N kl = n fe/ I e, A) n kl>bkl E b k i)P(N k i = n ki | B kl = b kl \ e, A) n kl>bkt E hiP(N k i = n k i | B kl = b ki \ e, A) nkl>bki £ n ; C~ (1 ~ tfct|)Ac( *’ n K 1 — g*ei)Ac(*r,l)]"fcl bkl T e -('-ekei)\c(k,l) [(l-f fc e i )Ac(fc,Q] n 'fci kl A _ e (n' 1! l kU' (C.12) where n' fci = n fci b kt . The first summation is equivalent to the mean of a Poisson random variable with parameter (1 — e k i)Xc(k,l), and the second summation is the sum of the probability over all possible values of n' kl . Therefore, E[N k i | B k i — b k i ; e, A] — (1 — e k ei)\c(k, l ) + b k \. (C.13) Given the previous estimates e m and A m , "3 = E[Nu I B u = bu; 6 m , A">] = (1 l)+b u , (C.14) for A; = 1, 2, ...,/,/ e fan fc . C.3 Poisson Model III Because Poisson Model III has the exact same functional form as Poisson Model II, the E step associated with Poisson Model III follows from the same derivation as that associated with Poisson Model II. Therefore, the conditional expectation of the log-likelihood function of the complete data X = [B,N] given the incomplete data

PAGE 102

95 and current estimates p m and A m is U"(p, A; p m , A m ) = E E [bid log(p k pi) + (rift b kt ) log(l p kPl )} ^ /EicUi^ + E E [— Ao(A:,/) +n™[ log \a(k,l)] + C 3 , k=1 fefan* where C 3 is independent of p and A, and < = E[JV« I B tl = h,: p"‘. A m ] = (1 ,^ P T)\ m a(kA) + b tl , (C.15) (C.16) for k = 1, 2, e fan*.

PAGE 103

REFERENCES [1] K. Lange, and R. Carson, “EM Reconstruction Algorithms for Emission and Transmission Tomography,” Journal of Computer Assisted Tomography , vol. 8, no. 2, pp. 306-316, Apr. 1984. [2] H. Hiriyannaiah, “X-Ray Computed Tomography for Medical Imaging,” IEEE Signal Processing Magazine , pp. 42-59, March, 1997. [3] K. Lee, K. Jung, Y. Yoon, I. Jeung, and K. Jeong, “The Development of Planar Imaging Technique for Spray Characterization,” Conference on Lasers and ElectroOptics, vol. 2, pp. 163-164, Pacific Rim, 1999. [4] W. Klein, H. Barrett, I. Pang, D. Patton, M. Rogulski, J. Sain, and W. Smith, “FASTSPECT: Electrical and Mechanical Design of a High-Resolution Dynamic SPECT Imager,” IEEE Medical Imaging Conference, vol. 2, pp. 931-933, Oct 1995. [5] E. Garcia, T. Faber, J. Galt, C. Cooke, and R. Folks, “Advances in Nuclear Emission PET and SPECT Imaging,” IEEE Engineering in Medicine and Biology Magazine, vol. 19, no. 5. pp. 21-33, Sep.-Oct. 2000. [6] K. Iwata, M. Wu, and B. Hasegawa, “Design of Combined X-Ray CT and SPECT Systems for Small Animals,” Nuclear Science Symposium, vol. 3, pp. 1608-1612 1999. [7] J. Ollinger, and J. Fessler, “Positron-Emission Tomography,” IEEE Signal Processing Magazine, vol. 14, no. 1, pp. 43-55, Jan. 1997. [8] C. Wu, “Attenuation Correction Methods for Positron Emission Tomography,” Ph. D. Dissertation, University of Florida, pp. 4-5, May, 1998. [9] M. Defrise, D. Townsend, D. Bailey, A. Geissbuhler, C. Michel, and T. Jones, “A Normalization Technique for 3D PET Data,” Phy. Med. Biol., vol. 36, no. 7, pp. 939-952, 1991. [10] E. Hoffman, S. Huang, M. Phelps, and D. Kuhl, “Quantization in PositronEmission Computed Tomography,” Journal of Computer Assisted Tomography, vol. 5, pp. 391-400, 1981. [11] M. Daube-Witherspoon and R. Carson, “Unified Deadtime Correction Model for PET,” IEEE Trans. Medical Imaging, vol. 10, pp. 267-275, 1991. 96

PAGE 104

97 [12] D. Bailey, D. Townsend, P. Kinahan, S. Grootoonk, and T. Jones, “An Investigation of Factors Affecting Detector and Geometric Correction in Normalization of 3-D PET Data,” IEEE Trans, on Nuclear Science , vol. 43, no. 6, pp. 3300-3307 Dec. 1996. [13] S. Derenzo, “Mathematical Removal of Positron Range Blurring in High Resolution Tomography,” IEEE Trans. Nuclear Science , vol. 33, pp. 565-569, 1986. [14] L. Shepp, and Y. Vardi, “Maximum Likelihood Reconstruction for Emission Tomography,” IEEE Trans, on Medical Imaging, vol. 1-1, no. 2, pp. 113-122 Oct. 1982. [15] G. Grimmett, and D. Stirzaker, Probability and Random Processes, pp. 255, New York: Oxford 2001. [16] E. Hoffman, T. Guerrero, G. Germano, W. Digby, and M. Dahlbom, “PET System Calibrations and Corrections for Quantitative and Spatially Accurate Images,” IEEE Trans, on Nuclear Science, vol. 36, no. 1, pp. 1108-1112 Feb 1989. [17] N. Ferreira, R. Trebossen, M. Gregoire, and B. Bendriem, “Influence of Malfunctioning Block Detectors on the Calculation of Single Detector Efficiencies in PET,” IEEE Trans, on Nuclear Science, vol. 46, no. 4, pp. 1062-1069 Aug 1999. [18] Y. Pawitan, S. Kohlmyer, T. Lewellen, and F. O’Sullivan, “PET System Calibration and Attenuation Correction,” IEEE Trans, on Nuclear Science, vol. 44, no. 3, pp. 1249-1253, June 1997. [19] F. Hermansen, T. Spinks, P. Camici, and A. Lammertsma, “Calculation of Single Detector Efficiencies and Extension of the Normalization Sinogram in PET,” Phys. Med. Biol., vol. 42, no. 6, pp. 1143-1154, June 1997. [20] J. Ollinger, “Detector Efficiency and Compton Scatter in Fully 3D PET,” IEEE Trans, on Nuclear Science, vol. 42, no. 4, pp. 1168-1173, Aug. 1995. [21] R. Badawi, M, Lodge, and P. Marsden, “Algorithms for Calculating Detector Efficiency Normalization Coefficients for True Coincidences in 3D PET,” Phy. Med. Biol. , 43. pp. 189-205, 1998. [22] D. Chesler, and C. Sterns, “Calibration of Detector Sensitivity in Positron Cameras,” IEEE Trans, on Medical Science, vol. 37, no. 2, pp. 768-772, 1990. [23] M. Casey, and E. Hoffman, “Quantitation in Positron Emission Computed Tomography: 7. A Technique to Reduce Noise in Accidental Coincidence Measurements and Coincidence Efficiency Calibration,” Journal of Computer Assisted Tomography, vol. 10, no. 5, pp. 845-850, Sep. 1986.

PAGE 105

98 [24] M. Casey, H. Gadagkar, and D. Newport, “A Component Based Method for Normalization in Volume PET,” Proc. 3rd International Meeting on Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine , Aix-IcsBains, pp. 67-71, 1995. [25] Z. Liang, “Detector Response Restoration in Image Reconstruction of High Resolution Positron Emission Tomography,” IEEE Trans, on Medical Imaging , vol. 13, no. 2, pp. 314-321, June 1994. [26] K. Wienhard, L. Eriksson, S. Grootoonk, M. Casey, U. Pietrzyk, and W. Heiss, “Performance Evaluation of the Positron Scanner ECAT EXACT,” Journal of Computer Assisted Tomography , vol. 16, no. 5, pp. 804-813, 1992. [27] Z. Cho, S. Juh, R. Friedenberg, W. Bunney, M. Buchsbaum, and E. Wong, “A New Approach to Very High Resolution Mini-Brain PET Using a Small Number of Large Detectors,” IEEE Trans, on Nuclear Science , vol. 37, no. 2. pp. 842-851, April 1990. [28] R. Huesman, E. Salmeron, and J. Baker, “Compensation for Crystal Penetration in High Resolution Positron Tomography,” IEEE Trans, on Nuclear Science, vol. 36, no. 1, pp. 1100-1107, Feb. 1989. [29] S. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, pp. 157-198, Englewood Cliffs, NJ: Prentice-Hall 1993. [30] A. Dempster, N. Laird, and D. Rubin, “Maximum Likelihood from Incomplete Data via the EM algorithm,” Ann. Roy. Stat. Soc., vol. 39, no. 1, pp. 1-38, Dec. 1977. [31] M. Feder, and E. Weistein, “Parameter Estimation of Superimposed Signals Using the EM Algorithm,” IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. 36, no. 4, pp. 477-489, Apr. 1988. [32] M. Bazaraa, H. Sherali, and C. Shetty, Nonlinear Programming: Theory and Algorithms, pp. 268-275, New York: Wiley 1979. [33] T. Abatzoglou, and B. O’Donnell, “Minimization by Coordinate Descent,” J. of Opt. Theory and Applications, vol. 36, no. 2, pp. 163-174, Feb. 1982. [34] Z. Luo, and P. Tseng, “On the Convergence of the Coordinate Descent Method for Convex Differentiable Minimization,” J. of Opt. Theory and Applications, vol. 72, no. 1, pp. 7-35, Jan. 1992. [35] P. Stoica, and T. Marzetta, “Parameter Estimation Problems with Singular Information Matrices,” IEEE Trans, on Signal Processing, vol. 49, no. 1, pp. 87-90, Jan. 2001. [36] R. Carroll, “PET Reconstruction Algorithm,” Master Thesis, Department of Mathematics, University of Florida, 1998.

PAGE 106

99 [37] Y. Vardi, L. Shepp, and Y. Kaufman, “A Statistical Model for Positron Emission Tomography”’ Journal of the American Statistical Association, vol. 80, no. 389, pp. 8-20, March 1985. [38] M. Madsen, “PET Attenuation Correction Using Mean Attenuation Coefficients: a Simulation Study,” IEEE Trans, on Nuclear Science, vol. 46, no. 6, pp. 21722176, Dec. 1999. [39] A. Papenfuss, G. O’Keefe, and A. Scott, “Segmented Attenuation Correction in Whole Body PET Using Neigh EM Clustering,” IEEE Nuclear Science Symposium Conference Record, vol. 2, no. 2, pp. 13/78-13/81, 2000. [40] A. De Pierro, “On the Relation Between the ISRA and the EM Algorithm for Positron Emission Tomography,” IEEE Trans, on Medical Imaging, vol. 12, no. 2, pp. 328-333, June, 1993. [41] H. Hudson, and R. Larkin, “Accelerated Image Reconstruction Using Ordered Subsets of Projection Data,” IEEE Trans, on Medical Imaging, vol. 13, no. 4, pp. 601-609, Dec. 1994. [42] K. Choi, Personal Communication, University of Florida. [43] B. Noble, and J. Daniel, Applied Linear Algebra, pp. 427, Englewood Cliffs, NJ: Prentice-Hall 1977. [44] Y. Viniotis, Probability and Random Processes, pp. 38, New York: McGraw-Hill 1998.

PAGE 107

BIOGRAPHICAL SKETCH I was born in Kaohsiung, a big modern city in Taiwan. After I earned my B.S. degree from the Electronics Engineering Department of National Jiao Tong University in Hsinchu, I was admitted by the graduate school of University of Florida and studied for Master and Ph. D. degrees in the Department of Electrical and Computer Engineering. My research focuses on digital signal processing. Since I joined the research group led by my advisor, Dr. Anderson, I had been working in the Statistical Signal Processing Lab and doing research on medical image processing. Apart from academics, I enjoy reading Chinese literature, listening to music, and playing basketball. 100

PAGE 108

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. /John M. M. Anderson, Chairman / Associate Professor of Electrical and Computer Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Associate Professor of Electrical Engineering, Hampton University I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Jian Professor o: Engineering lectrical and Computer I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. /— Panagote M. Pardalos Professor of Industrial and Systems Engineering

PAGE 109

This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 2001 Pramod P. Khargonekar Dean, College of Engineering Winfred M. Phillips Dean, Graduate School