
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00039257/00001
Material Information
 Title:
 Maximum likelihood estimation of detector efficiencies in positron emission tomography
 Creator:
 Lee, WenHsiung
 Publication Date:
 2001
 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 ) nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph.D.)University of Florida, 2001.
 Bibliography:
 Includes bibliographical references (leaves 9699).
 General Note:
 Printout.
 General Note:
 Vita.
 Statement of Responsibility:
 by WenHsiung 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 nonprofit 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
WENHSIUNG 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
WenHsiung 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 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.
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 NonCollinearity . . . . .
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 LOGLIKELIHOOD FUNCTIONS
A.1 Poisson Model I
A.2 Poisson Model II . . . . . . . .
A.3 Poisson Model III . . . . . . . .
B CRAMERRAO 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
WenHsiung 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 regionofinterest (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 nonuniformity, 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 maximumlikelihood 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 expectationmaximization algorithms, where the maximization step is solved using two optimization algorithms. As desired, the resulting algorithms have the property that the loglikelihood function is nondecreasing 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, Xray 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], singlephoton emission computed tomography (SPECT) [46], and positron emission tomography (PET) [7].
Planar imaging suffers from relatively low availability of radiopharmaceuticals, poor sensitivity of cameras, and, most importantly, the socalled 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 halflife radionuclides, such as 121 and 201TI. Because of the long halflives
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 aforementioned shortcomings and provides better images in the sense of efficiency, sensitivity and accuracy. In addition, positronemitting 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 halflives (2 minutes, 10 minutes, 20 minutes, and 110 minutes for 1C, 15o, 13N and 18F, respectively), onsite 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 highenergy (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 omnidirectional 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 crosssection 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 80100 cm and an axial extent of 1020 cm. The detectors are shielded from outside the field of view (or ROI) by relatively thick, lead endshields. Most scanners can switch between two operating modes. In the first mode, known as slicecollimated mode or 2dimensional (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 threedimensional (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 bismuthgermanate 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 deadtime [7].
1.3.1 NonCollinearity
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 noncollinearity 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 35 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 mispositioned 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 nonuniformities, detector efficiencies are not equal and less than 100% [12,13]. The efficiency of current PET detectors depends on a number of geometric and nongeometric 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 errorfree 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 wellknown 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)Zy ez(
= Z ! y ! ( . p ) Z . ( )z
Z=Y y!(z  y) Z!
10
_ 2P* (1 P)zY(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~ (lP)Z'(Z)'+Y
P( )=e  E 1
y z'=0
= (Pz)y ((1  P))/ z'=O
y' ___==e 2 (P2) _,4P)
y!
epz(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 24 hours and is usually performed once a month (see Fig. 1.3).
11
Figure 1.3: Blank scan.
Numerous methods [9,1624] 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 fansum 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 fansum algorithm may not provide accurate estimates for the detector efficiencies. Also, the fansum algorithm does not account for the Poisson nature of the blank scan data.
Several methods [9,1619] 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 fixedpoint 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 fansum algorithm, Ferreira's algorithm assumes noisefree 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 noisefree data is the algorithm by Defrise et al. [9]. The algorithm by Defrise et al. parallels the fansum 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 noisefree.
1.6 Outline of Research
All of the estimation methods we have discussed are suboptimal 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 [2528] 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 loglikelihood 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:
" Overparameterized approach, where the unknowns are the efficiencies of the
detector pairs, {Ekl}, and certain emission means.
" Sufficientlyparameterized approach, where unknowns are the efficiencies of individual detectors, {Ek}, and certain emission means.
Because the overparameterized 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 sufficientlyparameterized approach. Closed form expressions for the E step and M step of the EM algorithm are derived for our first model using the overparameterized formulation. However, simulations show that the overparameterized approach has poor performance. Thus, the results for the overparameterized 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 sufficientlyparameterized 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 fixedpoint iteration or the method of coordinate descent. All the resulting (EM) algorithms have
15
the property that the loglikelihood 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 fixedpoint 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 fansum 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 fansum 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 fieldofview. 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 midpoints 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 SiemensCTI 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
"71 a
"\ 
18
Table 2.1: Perpendicular distances of detector pairs for the first projection of the SiemensCTI 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 loglikelihood function
I
1 (E, A) A log P(B = b; E, A) = log ] ekAp(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 loglikelihood 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.
CramerRao Lower Bound
To compute the CramerRao 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 noninvertible 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 loglikelihood 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 loglikelihood 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 closedform 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 closedform 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 283287] or a fixedpoint iteration. See Abatzoglou et al. [33] or Luo et al. [34] for information on the convergence properties of CD method.
25
EMCD 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 1D optimization algorithms (e.g., dichotomous search method, golden section method, and Fibonacci search [32, pages 268275]). 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 1D 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 EMCD algorithm. To solve the maximization problem in (2.18), we use the dichotomous search method [32].
EMFP Algorithm
The second way we solve the problem (PE) is to use a fixedpoint iteration. Although it guarantees nonnegative detector efficiency estimates, the fixedpoint 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 fixedpoint iteration did have this property. An advantage of the fixedpoint iteration is that it converges faster in practice than the coordinate descent algorithm.
Ignoring the constraint in (PE), our approach is to use a fixedpoint 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'kjbk) ] = 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 lefthand side of (2.25) yields
1k bk l _ _ _ _ , k = 1 ,2 ....I1. (2 .26 )
1cfan, lEfank
The following fixedpoint 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 EMCD 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 fixedpoint iteration as the EMFP algorithm.
28
OverParameterized Approach
We now assume the unknowns are the set of parameters {6ki}. This case is referred to as overparameterized 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 overparameterized formulation is the same as the E step for the sufficientlyparameterized 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 (238) 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 noninvertibility 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 angleofvies (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 SiemensCTI 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 angleofview 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 Ik
Figure 2.9: The parameters {c(k, l)} for the SiemensCTI ECAT EXACT 921 scanner calculated according to (2.42). where 02(E, A) is the loglikelihood function:
#2(E, A) = log P(B = b; E, A) = log 1i fi ekL 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 noninvertible (see Appendix B).
2.2.2 Estimation Algorithms Associated with Poisson Model II
We show in Appendix C that the conditional expectation of the loglikelihood 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 loglikelihood 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) E1 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 EMCD2 algorithm. We refer to the estimation algorithm using (2.27) and (2.55) for Poisson Model II as the EMFP2 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 eLk' [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
PkeLk '(  eLky), (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  eLP) 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 cm1 [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 eL,,(,)ju (1  ekI))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  eLk (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
1k
160 170 180 190
Figure 3.6: The parameters {a(k, l)} for the SiemensCTI 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 loglikelihood function:
03(p, A) = log P(B = b; p, A) = log f1 11 ekpIa,)(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 noninvertible. Therefore, the CramerRao 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 loglikelihood 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)+nlogAa(k,1)]+C3,
k=1 IEfan k
(3.14)
where L'(B, N; p, A) is the loglikelihood 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 CDbased 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 fixedpoint 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 EMCD3 algorithm. We refer to the estimation algorithm using (3.21) and (3.22) for Poisson Model III as the EMFP3 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 fansum 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 EMCD and EMFP 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 EMCD2 and EMFP2 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 EMCD3 and EMFP3 algorithms is
b
A0 = k=1 Efank a (k, l)
k=1 lEfank
48
49
The steps of the EMCD and EMFP 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 EMFP algorithm converged much faster than the EMCD algorithm. Consequently, we used 6, = 6 = 101 for the EMFP algorithm and 6, = 6 = 10' for the EMCD algorithm. As for both the EMCD2 and EMFP2 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 EMFP2 algorithm and 6, = 6 = 102.5 for the EMCD2 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 SiemensCTI 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 2hour scan. Uniform, piecewise 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 EMCD3 and EMFP3 algorithms are based on Poisson Model III, which incorporates the effect of detector penetration, the EMCD3 and EMFP3 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 EMFP3 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 fansum, Ferreira's, EMCD, EMFP, EMCD2, and EMFP2 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 fansum 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 piecewise 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 fansum, Ferreira's, EMCD, EMFP, EMCD2, and EMFP2 algorithms are plotted in Figs. 4.2(b)4.2(g), respectively. The fansum 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 meansquare 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 fansum 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 Fansum Ferreira's EMCD EMFP EMCD2 EMFP2
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 Fansum Ferreira's EMCD EMFP EMCD2 EMFP2
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 fansum 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, EMCD, EMFP, EMCD2, and EMFP2 algorithms, respectively. Based on Fig. 4.5, we confer that all four of the proposed algorithms provide much better detector efficiency estimates than the fansum 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 Fansum Ferreira's EMCD EMFP EMCD2 EMFP2
MSE 100.8640 100.7207 100.2472 100.2704 100.2215 100.2070
tor efficiency estimates with smaller MSE than the fansum 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 fansum algorithm and Ferreira's algo
54
rithm in the sense that they provide estimates with smaller VR in every MC run than the fansum 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 fansum and Ferreira's algorithms. In terms of mean relative error, the EMFP algorithm has the best performance, and both the EMCD and EMCD2 algorithms are better than the fansum 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 fansum and Ferreira's algorithms. Specifically, the proposed algorithms provide estimates with smaller maximum and average absolute errors than the fansum 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, fansum 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 EMCD, EMFP, EMCD2, and EMFP2 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
Fansum 786.1761 785.5547 788.4473 794.3706 802.9343
Ferreira's 786.1472 785.5140 788.3954 794.3079 802.8612
EMCD 786.1328 785.4880 788.3587 794.2610 802.8044
EMFP 786.1363 785.4823 788.3440 794.2375 802.7723
EMCD2 786.2022 785.3626 788.0479 793.7720 802.1426
EMFP2 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 1620 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 SiemensCTI 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 fansum 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 fansum, Ferreira's, EMCD, EMFP, EMCD2, EMFP2 algorithms, respectively, are shown. The reconstructed emission images are indistinguishable. As the EMCD3 and EMFP3 algorithms are based on Poisson Model III which incorporates the effect of detector penetration, when applying the EMCD3 and EMFP3 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 EMCD3, EMFP3 algorithms and no normalization (i.e., {k = 1}) are plotted in Figs 4.17(a)(c). The reconstructed emission images for the EMCD3 and EMFP3 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 EMCD3 and EMFP3 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) fansum algorithm; (c) Ferreira's algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 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) fansum algorithm; (c) Ferreira's algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 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 104
  (a) Fansum algorithm x (b) Ferreira's algorithm
(c) EMCD
0 (d) EMFP
o (e) EMCD2
(f) EMFP2
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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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, piecewise uniform detector efficiencies, and:
(a) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm.
62
2
0
2
   (a) Fansum algorithm
x (b) Ferreira's algorithm
(c) EMCD
o (d) EMFP
o (e) EMCD2
A (f) EMFP2
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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm.
   (a) Fansum algorithm
x (b) Ferreira's algorithm
(c) EMCD
o (d) EMFP
o (e) EMCD2
(f) EMFP2
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) Fansum algorithm
(b) Ferreira's algorithm
(c) EMCD
0 (d) EMFP
o (e) EMCD2
A (f) EMFP2
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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) Fansum algorithm x (b) Ferreira's algorithm
1.016 (c) EMCD
0 (d) EMFP
3 (e) EMCD2
1.0159 A (f) EMFP2
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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) Fansum algorithm
x (b) Ferreira's algorithm
(c) EMCD
o (d) EMFP
o (e) EMCD2
A (f) EMFP2
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) fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm. Note that the detector efficiency estimates are normalized so that their mean equals the mean of true efficiencies.
68
   (a) Fansum algorithm x (b) Ferreira's algorithm
(c) EMCD o (d) EMFP
0 (e) EMCD2
A (f) EMFP2
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) fansum 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) EMCD algorithm; (b) EMFP algorithm;
(c) EMCD2 algorithm; and (d) EMFP2 algorithm.
72
12000
0 (a) No normalization
   (b) Fansum algorithm x (c) Ferreira's algorithm
10000 0 (d) EMCD
o (e) EMFP
0 (f) EMCD2
(g) EMFP2
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) fansum algorithm;
(c) Ferreira's algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 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 fansum algorithm; (b) Ferreira's algorithm; (c) EMCD algorithm; and (d) EMFP 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) EMCD2 algorithm; and (b) EMFP2 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) EMCD3 algorithm; (b) EMFP3 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 fixedpoint 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 fixedpoint 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 loglikelihood 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 fansum and Ferreira's algorithms. However, the reconstructed emission images for all of the estimation algorithms were visually iden
75
76
tical. In addition, the meansquare 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 EMCD3 and EMFP3 algorithms to have a large impact because Poisson Model III takes the effect of detector penetration into consideration. The reconstructed emission images for the fansum, Ferreira's, EMCD, EMFP, EMCD2, and EMFP2 algorithms are indistinguishable. When applying the EMCD3 and EMFP3 algorithms to reconstruct emission images, a different probability matrix P is used. The reconstructed emission images for the EMCD3 and EMFP3 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 EMCD3 and EMFP3 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 LOGLIKELIHOOD FUNCTIONS
A.1 Poisson Model I The loglikelihood 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 lefthand 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 loglikelihood 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 overparameterized approach, the same conclusion can be reached by a similar derivation. The loglikelihood 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 loglikelihood 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 loglikelihood function is not concave.
A.2 Poisson Model II
The loglikelihood 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
CRAMERRAO LOWER BOUNDS To compute the CramerRao lower bound, the Fisher information matrix must be computed. In this appendix, we derive the Fisher information matrix for the sufficientlyparameterized 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 ek(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 bA)  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 higherdimensional 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(9fAc(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 noninvertible. 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 noninvertible.
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 loglikelihood 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 loglikelihood 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 (1EkEl)"klbkle p(kI) nkl
ek l) p(k,L) (CkeLAp(k,I))bkl
e bkl!
(1fkl)Ap(k,) [(1 EkEL)Ap(k,l) Ikibkl (nkIbkl)!
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(1EkfI)Ap(k,) [(1kf)Ap(k,)]"k nk, > >0 nk
+ bkL Z: e( 1kE)Apk,1 [(1kE)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 loglikelihood 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(1ekel)kI
k=1 IEfank
+ log y eAc(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 loglikelihood 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) (C11)
eekel'Wck,') (e'f"\c(k,'))bk
 ~bkI,~E~E/b nkl
(1kEL)Ac(k,) [(1eke)Ac(kL)"kIbk1
(nkjbkj)!

Full Text 
PAGE 1
MAXIMUM LIKELIHOOD ESTIMATION OF DETECTOR EFFICIENCIES IN POSITRON EMISSION TOMOGRAPHY i i . ) !  By AQr >  WENHSIUNG 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 WenHsiung 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 NonCollinearity 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 LOGLIKELIHOOD FUNCTIONS 78 A.l Poisson Model I 78 A. 2 Poisson Model II 81 A. 3 Poisson Model III 82 B CRAMERRAO 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 WenHsiung 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 regionofinterest (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 nonuniformity, 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 maximumlikelihood 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 expectationmaximization algorithms, where the maximization step is solved using two optimization algorithms. As desired, the resulting algorithms have the property that the loglikelihood function is nondecreasing 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, Xray 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], singlephoton emission computed tomography (SPECT) [46], and positron emission tomography (PET) [7]. Planar imaging suffers from relatively low availability of radiopharmaceuticals, poor sensitivity of cameras, and, most importantly, the socalled 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 halflife radionuclides, such as 123 I and 201 TI. Because of the long halflives 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 aforementioned shortcomings and provides better images in the sense of efficiency, sensitivity and accuracy. In addition, positronemitting 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 halflives (2 minutes, 10 minutes, 20 minutes, and 110 minutes for 11 C, ls O, 13 N and 18 F, respectively), onsite 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 highenergy (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 omnidirectional 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 crosssection 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 80100 cm and an axial extent of 1020 cm. The detectors are shielded from outside the field of view (or ROI) by relatively thick, lead endshields. Most scanners can switch between two operating modes. In the first mode, known as slicecollimated mode or 2dimensional (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 threedimensional (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 bismuthgermanate 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 deadtime [7]. 1.3.1 NonCollinearity 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 noncollinearity 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 35 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 mispositioned 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 nonuniformities, detector efficiencies are not equal and less than 100% [12, 13]. The efficiency of current PET detectors depends on a number of geometric and nongeometric 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 errorfree 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 wellknown 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 ^ (ip) z y (zy = e y' 2 / ! ^ (* ~ y)\ Using the change of variable z' = zy, 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 (ip)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 24 hours and is usually performed once a month (see Fig. 1.3).
PAGE 18
11 Figure 1.3: Blank scan. Numerous methods [9, 1624] 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 fansum 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 fansum algorithm may not provide accurate estimates for the detector efficiencies. Also, the fansum algorithm does not account for the Poisson nature of the blank scan data. Several methods [9,1619] 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 fixedpoint 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 fansum algorithm, FerreiraÂ’s algorithm assumes noisefree 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 noisefree data is the algorithm by Defrise et al. [9]. The algorithm by Defrise et al. parallels the fansum 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 noisefree. 1.6 Outline of Research All of the estimation methods we have discussed are suboptimal 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 [2528] 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 loglikelihood 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: Â• Overparameterized approach, where the unknowns are the efficiencies of the detector pairs, {e*,;}, and certain emission means. Â• Sufficientlyparameterized approach, where unknowns are the efficiencies of individual detectors, {ejt}, and certain emission means. Because the overparameterized 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 sufficientlyparameterized approach. Closed form expressions for the E step and M step of the EM algorithm are derived for our first model using the overparameterized formulation. However, simulations show that the overparameterized approach has poor performance. Thus, the results for the overparameterized 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 sufficientlyparameterized 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 fixedpoint iteration or the method of coordinate descent. All the resulting (EM) algorithms have
PAGE 22
15 the property that the loglikelihood 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 fixedpoint 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 fansum 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 fansum 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 fieldofview. 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 midpoints 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 SiemensCTI 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 SiemensCTI ECAT EXACT 921 scanner. 1st projection member pair ( k , l) lk 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 loglikelihood 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 loglikelihood 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. CramerRao Lower Bound To compute the CramerRao 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 noninvertible 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 loglikelihood 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 loglikelihood 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 closedform 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 283287] or a fixedpoint iteration. See Abatzoglou et al. [33] or Luo et al. [34] for information on the convergence properties of CD method.
PAGE 32
25 EMCD 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 fcl ) 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 , Â• , zii, 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 1D optimization algorithms (e.g., dichotomous search method, golden section method, and Fibonacci search [32, pages 268275]). 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 1D 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 EMCD algorithm. To solve the maximization problem in (2.18), we use the dichotomous search method [32]. EMFP Algorithm The second way we solve the problem (P e ) is to use a fixedpoint iteration. Although it guarantees nonnegative detector efficiency estimates, the fixedpoint 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 fixedpoint iteration did have this property. An advantage of the fixedpoint iteration is that it converges faster in practice than the coordinate descent algorithm. Ignoring the constraint in (Pe), our approach is to use a fixedpoint 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 + (Â«SW 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 lefthand side of (2.25) yields Â£ Â£ T^7T> k = 1,2, . . . , i iefan*. iefan fc (2.26) The following fixedpoint iteration is then used to find an approximate solution to the system of equations in (2.26): m+l,i+l Â£ fefan* rHTft 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 EMCD 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 fixedpoint iteration as the EMFP algorithm.
PAGE 35
28 OverParameterized Approach We now assume the unknowns are the set of parameters {e fci }. This case is referred to as overparameterized, 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 overparameterized formulation is the same as the E step for the sufficientlyparameterized 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. (234) 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 noninvertibility 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 Â’ (240) 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 angleofvies (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 SiemensCTI 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 angleofview 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 SiemensCTI ECAT EXACT 921 scanner calculated according to (2.42). where 0 2 (e> A) is the loglikelihood 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 noninvertible (see Appendix B). 2.2.2 Estimation Algorithms Associated with Poisson Model II We show in Appendix C that the conditional expectation of the loglikelihood 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 , (248) fc=1 iefan fc where L' C (B, iV; e, A) is the loglikelihood 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 EMCD2 algorithm. We refer to the estimation algorithm using (2.27) and (2.55) for Poisson Model II as the EMFP2 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^ _ gL; (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 SiemensCTI 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 : Np,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 loglikelihood 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 pp)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 CDbased 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 EMCD3 algorithm. We refer to the estimation algorithm using (3.21) and (3.22) for Poisson Model III as the EMFP3 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 fansum 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 EMCD and EMFP 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 EMCD2 and EMFP2 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 EMCD3 and EMFP3 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 EMCD and EMFP 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 EMFP algorithm converged much faster than the EMCD algorithm. Consequently, we used S e = 8 = 10 11 for the EMFP algorithm and 8 e = 8 = 10 3 for the EMCD algorithm. As for both the EMCD2 and EMFP2 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 35 for the EMFP2 algorithm and 8 t = 5 = 10 2 5 for the EMCD2 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 SiemensCTI 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 2hour scan. Uniform, piecewise 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 EMCD3 and EMFP3 algorithms are based on Poisson Model III, which incorporates the effect of detector penetration, the EMCD3 and EMFP3 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 EMFP3 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 fansum, FerreiraÂ’s, EMCD, EMFP, EMCD2, and EMFP2 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 fansum 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 piecewise 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 fansum, FerreiraÂ’s, EMCD, EMFP, EMCD2, and EMFP2 algorithms are plotted in Figs. 4.2(b)4.2(g), respectively. The fansum 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 meansquare 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 fansum 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 Fansum FerreiraÂ’s EMCD EMFP EMCD2 EMFP2 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 Fansum FerreiraÂ’s EMCD EMFP EMCD2 EMFP2 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 fansum 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, EMCD, EMFP, EMCD2, and EMFP2 algorithms, respectively. Based on Fig. 4.5, we confer that all four of the proposed algorithms provide much better detector efficiency estimates than the fansum 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 Fansum FerreiraÂ’s EMCD EMFP EMCD2 EMFP2 MSE 100.8640 100.7207 100.2472 100.2704 100.2215 100.2070 tor efficiency estimates with smaller MSE than the fansum 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 fansum 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 fansum 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 fansum and FerreiraÂ’s algorithms. In terms of mean relative error, the EMFP algorithm has the best performance, and both the EMCD and EMCD2 algorithms are better than the fansum 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 fansum and FerreiraÂ’s algorithms. Specifically, the proposed algorithms provide estimates with smaller maximum and average absolute errors than the fansum 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, fansum 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 Fansum 786.1761 785.5547 788.4473 794.3706 802.9343 FerreiraÂ’s 786.1472 785.5140 788.3954 794.3079 802.8612 EMCD 786.1328 785.4880 788.3587 794.2610 802.8044 EMFP 786.1363 785.4823 788.3440 794.2375 802.7723 EMCD2 786.2022 785.3626 788.0479 793.7720 802.1426 EMFP2 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 1620 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 SiemensCTI 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 fansum 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 fansum, FerreiraÂ’s, EMCD, EMFP, EMCD2, EMFP2 algorithms, respectively, are shown. The reconstructed emission images are indistinguishable. As the EMCD3 and EMFP3 algorithms are based on Poisson Model III which incorporates the effect of detector penetration, when applying the EMCD3 and EMFP3 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 EMCD3, EMFP3 algorithms and no normalization (i.e., {p k = 1}) are plotted in Figs 4.17(a)(c). The reconstructed emission images for the EMCD3 and EMFP3 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 EMCD3 and EMFP3 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) fansum algorithm; (c) FerreiraÂ’s algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 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) fansum algorithm; (c) FerreiraÂ’s algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 algorithm.
PAGE 68
61 x 10 i r 5 Â— (a) Fansum algorithm X (b) FerreiraÂ’s algorithm (c) EMCD o (d) EMFP (e) EMCD2 A (f) EMFP2 4x 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm.
PAGE 69
62 o 1 1 1 1 Â— (a) Fansum algorithm _ X (b) FerreiraÂ’s algorithm (c) EMCD o (d) EMFP (e) EMCD2 A (f) EMFP2 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, piecewise uniform detector efficiencies, and: (a) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm.
PAGE 72
65 i.ii i.i 1.09 1.08 "i r ~i r Â— (a) Fansum algorithm X (b) FerreiraÂ’s algorithm (c) EMCD o (d) EMFP (e) EMCD2 A (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 algorithm.
PAGE 73
66 1.0161 1.016 1.0159 o 1.0158 <5 1 r Â”i r Â— (a) Fansum algorithm X (b) FerreiraÂ’s algorithm (c) EMCD o (d) EMFP (e) EMCD2 A (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; (d) EMFP algorithm; (e) EMCD2 algorithm; and (f) EMFP2 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) EMCD algorithm; (b) EMFP algorithm; (c) EMCD2 algorithm; and (d) EMFP2 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) fansum algorithm; (c) FerreiraÂ’s algorithm; (d) EMCD algorithm; (e) EMFP algorithm; (f) EMCD2 algorithm; and (g) EMFP2 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 fansum algorithm; (b) FerreiraÂ’s algorithm; (c) EMCD algorithm; and (d) EMFP algorithm. (a) (b) >* lAr Â§ . . Â§ v' Â’J. > t * Â«ti *L.. v. * * ;V f% ;u : $:?<* K Hf' Figure 4.16: Reconstructed emission images of plane 10 of a thorax phantom at the 20 th iteration of EMML algorithm using: (a) EMCD2 algorithm; and (b) EMFP2 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 fixedpoint 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 fixedpoint 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 loglikelihood 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 fansum 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 meansquare 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 EMCD3 and EMFP3 algorithms to have a large impact because Poisson Model III takes the effect of detector penetration into consideration. The reconstructed emission images for the fansum, FerreiraÂ’s, EMCD, EMFP, EMCD2, and EMFP2 algorithms are indistinguishable. When applying the EMCD3 and EMFP3 algorithms to reconstruct emission images, a different probability matrix P is used. The reconstructed emission images for the EMCD3 and EMFP3 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 EMCD3 and EMFP3 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 LOGLIKELIHOOD FUNCTIONS A.l Poisson Model I The loglikelihood 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 lefthand 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 loglikelihood 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 + tJ) 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 overparameterized, approach, the same conclusion can be reached by a similar derivation. The loglikelihood 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 loglikelihood function are dlogP(B = 6;e,A) , , hi , , Â„ T , _ , jr = 'W0 + Â— , 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) dtkld^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 loglikelihood function is not concave. A. 2 Poisson Model II The loglikelihood 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* (A10)
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 CRAMERRAO LOWER BOUNDS To compute the CramerRao lower bound, the Fisher information matrix must be computed. In this appendix, we derive the Fisher information matrix for the sufficientlyparameterized 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 higherdimensional 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 fc(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 noninvertible. 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 noninvertible.
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 loglikelihood 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 Â‘ , (lefcei) 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 loglikelihood 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(le 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(ltfce t )A p(ti , )]"Â«**! ( n k ib 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 ~tfet[)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 loglikelihood 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 *'(le 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^(le 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 loglikelihood 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 Â„(le k e,)\c(k.l) (nklbkl)'(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) [(lf 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 loglikelihood 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. 306316, Apr. 1984. [2] H. Hiriyannaiah, Â“XRay Computed Tomography for Medical Imaging,Â” IEEE Signal Processing Magazine , pp. 4259, 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. 163164, 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 HighResolution Dynamic SPECT Imager,Â” IEEE Medical Imaging Conference, vol. 2, pp. 931933, 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. 2133, Sep.Oct. 2000. [6] K. Iwata, M. Wu, and B. Hasegawa, Â“Design of Combined XRay CT and SPECT Systems for Small Animals,Â” Nuclear Science Symposium, vol. 3, pp. 16081612 1999. [7] J. Ollinger, and J. Fessler, Â“PositronEmission Tomography,Â” IEEE Signal Processing Magazine, vol. 14, no. 1, pp. 4355, Jan. 1997. [8] C. Wu, Â“Attenuation Correction Methods for Positron Emission Tomography,Â” Ph. D. Dissertation, University of Florida, pp. 45, 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. 939952, 1991. [10] E. Hoffman, S. Huang, M. Phelps, and D. Kuhl, Â“Quantization in PositronEmission Computed Tomography,Â” Journal of Computer Assisted Tomography, vol. 5, pp. 391400, 1981. [11] M. DaubeWitherspoon and R. Carson, Â“Unified Deadtime Correction Model for PET,Â” IEEE Trans. Medical Imaging, vol. 10, pp. 267275, 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 3D PET Data,Â” IEEE Trans, on Nuclear Science , vol. 43, no. 6, pp. 33003307 Dec. 1996. [13] S. Derenzo, Â“Mathematical Removal of Positron Range Blurring in High Resolution Tomography,Â” IEEE Trans. Nuclear Science , vol. 33, pp. 565569, 1986. [14] L. Shepp, and Y. Vardi, Â“Maximum Likelihood Reconstruction for Emission Tomography,Â” IEEE Trans, on Medical Imaging, vol. 11, no. 2, pp. 113122 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. 11081112 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. 10621069 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. 12491253, 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. 11431154, June 1997. [20] J. Ollinger, Â“Detector Efficiency and Compton Scatter in Fully 3D PET,Â” IEEE Trans, on Nuclear Science, vol. 42, no. 4, pp. 11681173, 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. 189205, 1998. [22] D. Chesler, and C. Sterns, Â“Calibration of Detector Sensitivity in Positron Cameras,Â” IEEE Trans, on Medical Science, vol. 37, no. 2, pp. 768772, 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. 845850, 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 , AixIcsBains, pp. 6771, 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. 314321, 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. 804813, 1992. [27] Z. Cho, S. Juh, R. Friedenberg, W. Bunney, M. Buchsbaum, and E. Wong, Â“A New Approach to Very High Resolution MiniBrain PET Using a Small Number of Large Detectors,Â” IEEE Trans, on Nuclear Science , vol. 37, no. 2. pp. 842851, 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. 11001107, Feb. 1989. [29] S. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, pp. 157198, Englewood Cliffs, NJ: PrenticeHall 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. 138, 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. 477489, Apr. 1988. [32] M. Bazaraa, H. Sherali, and C. Shetty, Nonlinear Programming: Theory and Algorithms, pp. 268275, 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. 163174, 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. 735, 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. 8790, 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. 820, 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/7813/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. 328333, 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. 601609, Dec. 1994. [42] K. Choi, Personal Communication, University of Florida. [43] B. Noble, and J. Daniel, Applied Linear Algebra, pp. 427, Englewood Cliffs, NJ: PrenticeHall 1977. [44] Y. Viniotis, Probability and Random Processes, pp. 38, New York: McGrawHill 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

