University of Southern California Dissertations and Theses

Direct wholebody Patlak and Logan image estimation from listmode PET data

(USC Thesis Other)

Direct wholebody Patlak and Logan image estimation from listmode PET data

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DIRECTWHOLEBODYPATLAKANDLOGANIMAGEESTIMATIONFROM
LISTMODEPETDATA
by
WentaoZhu
ADissertationPresentedtothe
FACULTYOFTHEUSCGRADUATESCHOOL
UNIVERSITYOFSOUTHERNCALIFORNIA
InFulfillmentofthe
RequirementsfortheDegree
DOCTOROFPHILOSOPHY
(ELECTRICALENGINEERING)
August2014
Copyright 2014 WentaoZhu
Contents
Acknowledgements xii
Abstract xiv
Chapter1
Introduction 1
1.1 PositronEmissionTomography . . . . . . . . . . . . . . . . . . . . . . 1
1.2
18
F-FDGPET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 StandardizedUptakeValue . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 WholebodyDynamicPET . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 TracerKineticModeling . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 GraphicalAnalysisofDynamicPETData . . . . . . . . . . . . . . . . 14
1.7 DirectEstimationofKineticParametricImagefromDynamicPETData 18
1.8 PETattenuationcorrectionwithcontrastCT . . . . . . . . . . . . . . . 22
1.9 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.10 Organizationofthedissertation . . . . . . . . . . . . . . . . . . . . . . 28
Chapter2
DirectPatlakestimationfromlistmodePETdata 29
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
ii
2.2.1 DirectPatlakEstimationfromtwoFrameList-ModeData . . . 33
2.2.2 TheCramer-Raoboundofthealgorithm . . . . . . . . . . . . . 37
2.2.3 EnforcingCountIndependentResolutioninPatlakImages . . . 39
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.3.1 CramerRaoAnalysis . . . . . . . . . . . . . . . . . . . . . . . 45
2.3.2 ROCAnalysis: DualTime-PointPatlakvs. FractionalSUV . . . 48
2.3.3 SimulatedImagingStudies . . . . . . . . . . . . . . . . . . . . 52
2.3.4 Count-InvariantResolution . . . . . . . . . . . . . . . . . . . . 55
2.3.5 MatchingResolutioninPatlakand%DSUVimages . . . . . . . 59
2.4 Discussionandconclusion . . . . . . . . . . . . . . . . . . . . . . . . 61
Chapter3
ClinicalapplicationofdirectwholebodyPatlakestimation 64
3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.3 Materialsandmethods . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.3.1 WholebodyPatlakAnalysis . . . . . . . . . . . . . . . . . . . 68
3.3.2 DataProcessing . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3.3 Patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.3.4 ImagingProtocol . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.5 StatisticalAnalysis . . . . . . . . . . . . . . . . . . . . . . . . 75
3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.4.1 InputFunctionEstimation . . . . . . . . . . . . . . . . . . . . 76
3.4.2 RandomandScatteredEvents . . . . . . . . . . . . . . . . . . 77
3.4.3 PatientMotionCompensation . . . . . . . . . . . . . . . . . . 77
3.4.4 Patlakv.s. SUVand%DSUV: . . . . . . . . . . . . . . . . . . 80
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
iii
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Chapter4
DirectLoganparameterestimationfromlistmodedataforreversibletracer 87
4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3 Materialsandmethods . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.3.1 DirectLoganparameterEstimation . . . . . . . . . . . . . . . . 90
4.3.2 Overcomingill-conditioning . . . . . . . . . . . . . . . . . . . 92
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.4.1 LoganandREplotvalidation . . . . . . . . . . . . . . . . . . . 96
4.4.2 Cramer-Raoanalysis . . . . . . . . . . . . . . . . . . . . . . . 96
4.4.3 Monte-Carlosimulation . . . . . . . . . . . . . . . . . . . . . . 99
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Chapter5
JointestimationofactivityandattenuationforcontrastenhancedPET/CT106
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.2 Materialsandmethods . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.2.1 BilinearmappingfromCTtoattenuationcoefficient . . . . . . . 110
5.2.2 Jointestimationofemissionandattenuation . . . . . . . . . . . 111
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.3.1 2DhistogramofnoncontrastedandcontrastCT . . . . . . . . . 113
5.3.2 3Dsimulationforestimationofemissionandattenuationmap . 114
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
iv
Chapter6
ConclusionandFuturework 119
6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.2 Futurework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
ReferenceList 125
v
ListofFigures
1.1 Physical principle of PET: positron generation, the annihilation of a
positionwithanelectron,andthegenerationofaphotonpair. . . . . . . 2
1.2 PET system 2D illustration. The annihilation happened in the object
generates a pair of photons traveling at the opposite direction. The two
photonsarethendetectedbythedetectorpair. . . . . . . . . . . . . . . 3
1.3 lineintegralsofanimageindifferentviewsandthecorrespondingRadon
transform,witharrowindicatingthe3lineprofilesshownontheleft. . . 4
1.4 WholebodyPETCTscanillustration. Thepatientissteppedthroughthe
scanneratmultiplebedpositions. . . . . . . . . . . . . . . . . . . . . . 9
1.5 2-tissue compartment model for kinetic analysis. C
p
(t) is the activity
in the plasma. C
1
(t) and C
2
(t) correspond to the activity in 2 tissue
compartmentsattimet. . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1 HistogramoftheestimatedPatlakslopefittingaGaussiancurve . . . . 47
2.2 Relative standard deviation of Patlak slope and intercept for list-mode
andsinogrambasedestimationmethodwithdifferentdurationofacqui-
sitionframes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.3 An example of the distribution of Patlak slope parameter (upper) and
%DSUV(lower)fortumorandbackground(BG) . . . . . . . . . . . . 50
vi
2.4 Detection of tumors relative to background using rate parameters from
Table 2.1. We compare area under ROC curves (AUC) for Patlak slope
(sinogram and list-mode based estimation) and %DSUV as a function
of (a) count rate, (b) duration of first frame, and (c) first frame start
time. Panel (d) shows the area under the ROC curves averaged over
simulationsfordifferentparametersetslistedinTable2.1 . . . . . . . . 51
2.5 (a)Uniformcylinder(background)withfivecylinders(tumor)ofdiffer-
entsizes. (b)Themean(over100MonteCarlotrials)ofthePatlakslope
image reconstructed from two 5 min frame list-mode data starting at 40
and80minwithTACsasinFig. 2.6. . . . . . . . . . . . . . . . . . . . 53
2.6 Simulated time-activity-curves for tumor, background, and the blood
inputfunctionwithkineticparametersfromrow1inTable2.1. . . . . . 54
2.7 Mean and s.d. (error bars) for each of five ROIs shown in Fig. 2.5(a)
for list-mode Patlak and %DSUV. %DSUV values were computed by
applying(2.39)toactivityineachframeafteraveragingovertheROI. . 55
2.8 Relative s.d. of Patlak slope estimate computed from 100 Monte-Carlo
trialsforthebiggesttumor. . . . . . . . . . . . . . . . . . . . . . . . . 56
2.9 AreaundertheROCcurveforkineticparametersets1,2,3inTable2.1
forPatlakand%DSUVbasedtumordetection. Seetextfordetails. . . . 57
2.10 (left) Digital phantom for validating count-independent resolution and
(b)profilethroughthe2pointsourcesinthephantom . . . . . . . . . . 58
2.11 Profiles through one reconstructed Patlak slope image with (left) spa-
tially invariant smoothing and (right) spatially variant smoothing for
count-independentresolutionusing(2.37). . . . . . . . . . . . . . . . . 59
2.12 ResolutionversuslogofsmoothingparameterinstaticMAPreconstruc-
tionscomputedatthecenterofhottumorforthephantominFig. 2.5. . 60
vii
3.1 (A)theinputfunctionfrom7liverscanpatientsand(B)wholetimeesti-
mationoftheinputfunctionforoneofthewholebodypatients. . . . . . 77
3.2 The2integratedscattersinogramcorrespondingtothesecondbedposi-
tioninthe50thdirectplane.. . . . . . . . . . . . . . . . . . . . . . . . 78
3.3 Difference of 2 MAP reconstructed images of frame 1 and frame 2
before(A,B)andafterregistration(C,D) . . . . . . . . . . . . . . . . 79
3.4 Patlakinterceptimages: (A,B)scattercorrectiononlyforthefirstframe
data, no motion correction. (C, D) scatter, random correction, motion
corrected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.5 ExampleofreconstructedSUV1,SUV2,Patlakslopeandinterceptimages
80
3.6 Coronal view of SUV1 (A), SUV2 (B), Patlak slope (C), and Patlak
intercept(D)images. Thearrowsindicatesthetumorintheleftbreast. . 83
3.7 Coronal view of SUV1 (A), SUV2 (B), Patlak slope (C), and Patlak
intercept(D)images. Thearrowsindicatethetumorintheliver. . . . . 84
4.1 ThedynamicchangeofC
2
(t)=(C
1
(t)+C
2
(t))andC
1
(t)=C
p
(t),C
2
(t)=C
p
(t)
asafunctionoftime. . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.2 Thebloodinputfunctionandthetime-activity-curveforthetumorregion
forthefirstrowinTable 4.1. . . . . . . . . . . . . . . . . . . . . . . . 98
4.3 standarddeviationcomputedfromCramerRaoanalysisandsimulation,
bothwithourmethodandthesinogrammethod. . . . . . . . . . . . . . 100
4.4 (a) Uniform cylinder (background) with five cylinders (tumor) of dif-
ferent sizes. (b) The mean (over 100 Monte Carlo trials) of the Logan
slope image reconstructed with our method with TACs as in Fig. 4.5.
(c) The mean (over 100 Monte Carlo trials) of the Logan slope image
reconstructedwiththesinogrambasedapproachwithTACsasinFig. 4.5.102
viii
4.5 Simulated time-activity-curves for tumor, background, and the blood
inputfunction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.6 profilethroughthesmallesttumorrodforthemeanLoganslopeimages
computedwithourmethodandthesinogrambasedmethod. . . . . . . . 104
5.1 thestandardconvertingcurvefromSiemens(blue)andourimplemented
curve for contrast enhanced CT (red color indicating the change made
uponoriginalmappingcurve). . . . . . . . . . . . . . . . . . . . . . . 111
5.2 the 2D histogram of pre-phase with no contrast agent and the venous
phaseafterthepatientconsumescontrastagent. The2peakscorrespond
totheincreaseofCTvaluesinsofttissuessuchasliverandheart. . . . . 114
5.3 Thereconsructedemissionimagewithcorrectattenuationmap(top)and
theoriginalnoncontrastedattenuationmap(bottom). . . . . . . . . . . 115
5.4 Thereconstructedemissionandattenuationimagewiththesimplerescal-
ingmethod(top)andtheoneswithourjointestimationalgorithm(bottom).117
ix
ListofTables
2.1 Pairs of kinetic model parameters for background and tumors identified
intheliterature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.2 Comparisonofthes.d. ofPatlakslopeobtainedfromtheCRboundwith
mean and s.d. computed by Monte Carlo simulation with maximum
likelihoodparameterestimation . . . . . . . . . . . . . . . . . . . . . . 47
2.3 CRC at 2 point sources for the low count level and the uniform spatial
penalty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.4 CRC at 2 point sources for low and high count levels and our modified
spatialpenalty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.1 Patlak parameters, SUV and %DSUV values from liver, thigh muscle,
andlumbarspine(L2)ROIs . . . . . . . . . . . . . . . . . . . . . . . . 81
3.2 Contrast of breast tumor and liver tumor against background for SUV
andPatlakmethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.1 reversible kinetic parameters for several brain or body regions reported
inliterature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.2 Comparison of the relative s.d. (s.d. divided by the mean) of Logan
parameters obtained from the CR bound with and without cumulated
sinogramdatafromtime=0-mintothesteadystate60-min . . . . . . . 99
4.3 mean and s.d. of the 5 tumor rod ROIs in 100 Monte-Carlo simulation
forourmethodandthesinogrambasedmethod . . . . . . . . . . . . . 103
x
5.1 Comparisonoftherelativeerrorofthesimplerescalingmethodandour
jointestimationmethodinestimatingtheattenuationcoefficientandthe
emission activity. The numbers are relative difference of the mean of
ROIsforestimatedvaluesandthegroundtruth. . . . . . . . . . . . . . 116
xi
Acknowledgements
IwouldliketofirstthankmyadvisorDr. RichardM.Leahyforhisguidanceandsupport
inmyresearchthroughouttheseyears. Hisinsightfulperspectiveandrichexperiencein
PET image reconstruction as well as other research directions in the signal and image
processing field made my experience in University of Southern California very benefi-
cial. Heistrulyagoodmentorineducatingspecializedresearcherswithseriousattitudes
towardsresearch. Ilearnedmuchfromhim.
I also want to thank Dr. Peter S. Conti, Dr. Bing Bai, and Dr. Quanzheng Li for
their invaluable comments and suggestions to my research and their involvement in my
research work. I would like to thank Dr. Conti particularly for his collaboration on
the wholebody data acquisition protocol and providing clinical data. Without his help
the application of my research on clinical data is not possible. Dr. Bai and Dr. Li are
also very generous in having discussion with me on my research projects and actively
involvinginthedetailedprojects.
During my internship in Siemens in 2012, Dr. Michael Casey, Dr. Christian Michel
and many other researchers offered me great help in my participated research, and con-
tinuouslyhelpedmethroughtherestofmyPh.D.career.
IntheDepartmentofElectricalEngineeringattheUniversityofSouthernCalifornia,
Dr. Krishna Nayak, Dr. Justin Haldar, and Dr. Urbashi Mitra all gave their very helpful
suggestionsandideastomyresearch.
xii
Finally, special thanks to my parents and my wife Shuran Zhang. I can not achieve
thesuccesswithoutyourmoralsupport.
TheresearchissupportedbyNIHgrantsR01EB013293andR01EB010197.
xiii
Abstract
We investigate using list-mode PET data to perform Patlak and Logan modeling. We
first propose an estimation method for irreversible tracers (where Patlak modeling is
applied) for dynamic PET studies in which we compute voxel-wise estimates of Pat-
lak parameters using two frames of data. Our approach directly uses list-mode arrival
time for each event to estimate the Patlak parametric image. We use a penalized like-
lihood method in which the penalty function uses spatially variant weighting to ensure
a count independent local impulse response. We evaluate performance of the method
in comparison to fractional changes in standardized uptake value (%DSUV) between
thetwoframesusingCramer-RaoanalysisandMonteCarlosimulation. Receiveroper-
ating characteristic (ROC) curves are used to compare performance in differentiating
tumors from background based on the dynamic data sets. Using area under the ROC
curve as a performance metric, we show superior performance of Patlak compared to
%DSUV over a wide range of dynamic data sets and parameters. These results suggest
thatPatlakanalysiscanbeusedfordynamicdual-framePETdataandmayleadtobetter
quantitativeresultsrelativeto%DSUVmethod.
Wethenextendourmethodtoclinicalwholebodypatientdatawith
18
F-FDG,which
is a commonly used irreversible tracer for oncologic applications. Several practical
issues such as estimation of dynamic randoms and scatters rate functions, motion cor-
rection,andtheestimationofbloodinputfunctionarediscussedandeffectivelyhandled.
xiv
Anewwhole-bodyPETscanprotocolisproposed,wherethepatientissteppedthrough
the PET scanner twice, after which the Patlak parametric image for each bed position
is estimated from the listmode data of the two frames and stitched together to form a
wholebodyPatlakparametricimage. Tenpatientswerescannedusingthisprotocol. We
then draw liver, muscle, and spine regions of interest (ROIs) for these patients and the
meanandstandarddeviation(s.d.) forthePatlakslope(netinfluxrate)andinterceptare
computed. WealsocomputetheSUVvaluesandthes.d. forthesameROIsatthe2cor-
respondingframestime,afterwhichthefractionalSUV(%DSUV)iscalculated. Results
showthatPatlakmethodhaslessinter-subjectvariationandhighertumor-to-background
contrast. Overall,thePatlakestimatesprovidequantitativeparametersbetterdescribing
tracer uptake than SUV. This may be useful for cancer staging, treatment planning and
assessingresponsetotherapy.
We also investigate the direct estimation of Logan parameters (where a reversible
tracer is used) from listmode PET data. The ill-conditioning nature of listmode recon-
structionforPETdatafromareversibleradioactivetracerisexploredanddiscussed. We
derive a modified scheme improving the conditioning of the problem by utilizing data
beforethesteadystate. Thismethoddiffersfromapreviouslyproposedsinogrambased
reconstruction methods by other researchers in that it follows the true likelihood of the
data. Cramer-Rao analysis and Monte-Carlo simulation demonstrate the effectiveness
ofourmethod.
NowadaysmanyPET/CTscansareperformedwithcontrastagentsaimingatenhanc-
ingsofttissuecontrastinCTimages. ThecontrastagentwillincreasetheCTvaluesig-
nificantly especially in the blood rich soft tissues such as liver and heart. However the
liner attenuation coefficient of contrast agent is only slightly higher than water at 511
keV.Asaresult,usingthesamemappingfromlowtohighenergiesfornon-contrastCT
willresultinover-estimationofattenuationinregionscontainscontrastregionforPET.
xv
ThisleadstoinaccuraciesinPETimages. Wedescribesajointestimationalgorithmfor
emission and attenuation maps using contrast enhanced CT as a prior. A comparison
with simple rescaling scheme shows that our method can improve quantitative image
accuracywhenusingCTwithcontrasttocomputeattenuationcoefficients.
xvi
Chapter1
Introduction
1.1 PositronEmissionTomography
Positron emission tomography (PET) is a powerful technique for staging cancer and
assessingresponsetotherapy. ItdiffersfromanatomicalimagingtechniquessuchasX-
ray computed tomography (X-ray CT) and magnetic resonance imaging (MRI) in that
it reflects different tracer uptake characteristics of various organs and tissues, thus it is
a functional imaging technique. The radioactive tracer is injected before a PET scan
starts and the changes in the spatial distribution of the tracer are then imaged over one
or more time intervals by detecting the photon pairs produced as a result of positron
emission and annihilation. One of the most commonly used tracer, and that on which
we focus for most of this thesis, is 2-deoxy-2-(
18
F)fluoro-D-glucose (
18
F-FDG), which
is a glucose analog that accumulates in cells in proportion to their metabolic rate.
18
F-
FDGPETisanimportantclinicaltoolfordetecting,monitoringandstagingcancerand
assessingtherapeuticresponse[1],[2].
ThephysicalprincipleofPETisillustratedinFig. 1.1. Positronsareproducedfrom
the nuclei of the radioisotopes. The annihilation of the positron with an electron after it
travelsashortrange(positronrange)willproduceapairofhighenergy(511keV)pho-
tons. Thispairofphotonstravelsinapproximatelyoppositedirections,andaredetected
in rings of scintillation-based or solid-state detectors that surround the patient. Usually
a narrow time window is set around 10ns. Two detection events occurring within this
1
time window are called coincidence and are assumed to come from the same annihi-
lation. The annihilation, as a result, is assumed to happen along the line of response
(LOR) connecting the two detectors (Fig. 1.2). Recently Time-of-Flight (TOF) PET
systems have been introduced which can better localize the annihilation position along
the LOR by measuring the difference in arrival times at the two detectors. The finer
spatial localization of the events improves the noise properties and consequently can
enhancelesiondetectionperformance[3].
Figure1.1: PhysicalprincipleofPET:positrongeneration,theannihilationofaposition
withanelectron,andthegenerationofaphotonpair.
ImagereconstructionisperformedafterdatacollectiontoformthePETimagewhich
usuallyrepresentsthedistributionofthetracer. Thesetoflineintegralsofa2Dimageat
differentanglesofviewsformits2DRadontransform[4]. Theselineintegralsaregen-
erated for each LOR in the PET scanner by summing (or ”histogramming”) all events
detected along that LOR. Figure 1.3 shows the line integrals of an image at 3 different
2
Figure 1.2: PET system 2D illustration. The annihilation happened in the object gen-
erates a pair of photons traveling at the opposite direction. The two photons are then
detectedbythedetectorpair.
viewsandthecorrespondingRadontransform,witharrowindicatingthe3profiles. The
Radon transform of image is also referred to as a sinogram in the context of medical
imaging. The classical analytical image reconstruction methods developed for 2D data
reconstruction is filtered backprojection (FBP) algorithm [4]. It can be shown that an
unfiltered backprojection yields a blurry image from the original image’s convolution
withaconefilter. ThereforetheFBPmethodapplyarampfilterinthesinogramdomain
to deconvolve with the cone filter, and then performs the backprojection. For 3D data,
analytical inversion formula for the projection data exist, however there is missing data
in oblique directions. Reconstruction can be performed with 3D reconstruction algo-
rithm such as 3DRP [5]. Conventionally, an initial image estimate is first obtained by
3
2D image reconstruction with direct planes only. The missing oblique lines of integrals
arethenestimatedbyforwardprojectionoftheinitialimage,andmergedwiththemea-
suredobliqueLORstoformthecomplete2Dprojections. Finallythe3Dbackprojection
isperformedafterfilteringdatawiththe2DColsher’sfilter. Anotherwaytoreconstruct
3D dataset is to transform the 3D dataset into a 2D dataset using one of the rebinning
methodssuchassingle-slicerebinning(SSRB)[6]andFourierrebinning(FORE).
Figure 1.3: line integrals of an image in different views and the corresponding Radon
transform,witharrowindicatingthe3lineprofilesshownontheleft.
Although the analytical approaches are fast, the line integral model may limit the
accuracy because in reality the model relating the image to the sinogram (the projec-
tion model) is more complex than simple mathematical line integration and the data is
corrupted by statistical noise. Researchers therefore developed statistical approaches
to image reconstruction based on the combination of an improved model of the system
response and a statistical model that reflects the photon-limited nature of the data [7].
The image is reconstructed by optimizing a cost function such as the log-likelihood
[8] or penalized log-likelihood (or equivalently, maximum a posteriori or MAP)
4
[9]. The cost function is optimized using numerical algorithms such as Expectation-
Maximization (EM) or Ordered-Subsets-Expectation-Maximization (OSEM) [10] (for
maximumlikelihood)andgeneralizedEMorconjugategradient(forMAP)[9].
Ordered-Subsets algorithm divides the whole data to multiple subsets and updates
the image using data from one subset at a time. A cyclic pass through all the subsets
constitutes one iteration. It is usually faster than the algorithms using all the data in
providing quantitatively meaningful images but does not maximize the likelihood of
whole dataset, and will lead to limit cycles and noisy reconstructed image especially
when the iteration number increases. This is partially because the gradient calculation
based on only one subset is not exact in all iterations. Effective methods have been
developedtohandletheconvergenceissueforordered-subsetsmethods,bydiminishing
the stepsize in the OS iterations to form a relaxed OS algorithm [11], or switching to
nonOSalgorithmsafterseveralOSiterations.
In this thesis, we develop our image reconstruction method based on penalized
maximum-likelihood, and use conjugate gradient method to optimize the cost function.
We adopt the penalized-ML instead of ML because the true ML solution can be very
noisy, therefore a regularization can be helpful in generating more stabilized images.
Thepenaltyaimsatreducingdramaticchangeofintensitiesinneighboringpixels.
1.2
18
F-FDGPET
18
F-FDG is a commonly used irreversible tracer in PET imaging. It is a glucose ana-
log with the positron-emitting radioactive isotope
18
F substituting the normal hydroxyl
group in the glucose molecule. It is utilized by cells similarly to glucose, and therefore
reflects glucose metabolic rate with high uptake in highly active regions such as brain,
heart, and tumors. Unlike glucose which is completely metabolized in cells, FDG is
5
only partially metabolized to form
18
F-FDG-6-phosphate after entering the cell, where
itwillbetrappedandwillnotbefurthermetabolized. Therefore,thespatialdistribution
of 18F-FDG is a good indicator of the distribution of glucose uptake or metabolic rate.
Because malignant tumors are usually more active metabolically than normal tissues,
they trap more
18
F-FDG-6-phosphate and consequently show higher intensity in PET
images. For this reason,
18
F-FDG is now widely used in the clinical setting for tumor
detectionandcancerstaging.
1.3 StandardizedUptakeValue
Currentlythequantificationof
18
F-FDGPETinclinicaluseprimarilyreliesonthemea-
surement of standardized uptake value (SUV) [12]. The formula for computing SUV is
giveninEquation 1.1,whereC(t)isthetissueradioactivityconcentration. Thedoseand
body weight are the normalization factors. Over the last decade, developments of PET
instrumentation and clinical practice have made SUV a useful tool to quantify tumor
uptake. For example, criteria for evaluating response to therapy in solid tumors based
onSUVvalueshavebeenproposed[13].
SUV =
C(t)
dose=weight
(1.1)
Nonetheless, there are limitations to the application of SUV and it should be used
withcaution. SUVmeasuresthetumoruptakeandmaynotreflectthetruepharmacoki-
netics of the tracer. For example, the uptake of
18
F-FDG is affected by many factors.
The uptake time prior to scanning may significantly influence SUV values and thus the
ability to reliably determine response to therapy [14]. SUV measurements also change
withrespecttodifferentimagereconstructionalgorithms,amongdifferentclinicalsites.
6
Masa-Ah et al [15] found that SUV values correlate with factors such as dose extrava-
sations, attenuation parameters, partial volume effect, plasma glucose level in blood.
Wiyaporn et al showed that SUV values increase with the number of iterations of the
reconstructionalgorithmandarealsoaffectedbytumorsize[16]. Otherfactorssuchas
thebloodflow,andrateofclearancemayalsoaffectthemeasuredSUV.
Among the various SUV based biomarkers, one of the most commonly used ones
is SUV
mean
computed as the average SUV value in a ROI. However some cancer treat-
mentresponseisusuallyassessedwithSUV
max
insteadofthemeanSUVvalue,because
this provides lower interobserver variability. On the other hand, there are also contra-
dictory findings that SUV
mean
is more reliable than SUV
max
. For example, Nahmias
[17]reportedthatSUV
max
hasworsereproducibilitythantheSUV
mean
. AnotherSUV-
oriented biomarker is SUV
peak
, which is calculated as the ”maximal average”, by mov-
ingavolumeof 1mLintheROIuntilthesmallvolumehasthemaximalVOIaverage.
Boellaard[18]suggestedthatthecombinationofSUV
max
andSUV
peak
shouldbecome
the standard approach in multicentre FDG PET/CT studies. However, Lodge showed
thatSUV
peak
islesssensitivetoimagecharacteristics[19].
Researchers have also looked at changes in SUV value between two acquisition
frames, expressed as percentage fractional difference: %DSUV = 100 (SUV
2

SUV
1
)=SUV
1
%,whereSUV
1
andSUV
2
aretheSUVvaluesat2timepoints. Byusing
fractional changes, these measures become invariant to the scaling by dose and body-
weightusedinstandardSUV.Forthisreason%DSUVreflectsonlydynamicchangesin
uptakeandisrobusttotheuncertaintiesgoverningtherelationshipbetweenFDGuptake
andpatientweightthatmakestandardSUVvaluesonlysemi-quantitative. Alkhawaldeh
et al. [20] found that dual time-point PET can improve diagnostic accuracy relative
to standard single frame SUV, increasing sensitivity and specificity for malignant lung
nodules diagnostic, especially for small lung lesions that have low SUVs. Prieto et al.
7
[21] showed that dual time-point PET can improve sensitivity for the identification and
volume delineation of high-grade brain tumors compared with standard PET studies.
Zhuang et al. [22] and Kumar et al. [23] found that malignant lesions showed a sig-
nificant increase in SUV over time, while benign lesions showed a decrease. However,
negativeresultshavealsobeenreported: Cloranetal. [24]reportedthatdualtime-point
FDGPETmaynotbeofbenefitintheassessmentofpulmonarynoduleswithmaximum
SUV of less than 2.5 on initial imaging Chan et al. [25] suggested that dual time-point
FDG-PETdoesnotimprovetheoverallevaluationofpulmonarylesions.
Considering the various biomarkers related to SUV and the controversy among dif-
ferentfindingsaboutwhicharesuperior,thestaticPETimageandthesemi-quantitative
biomarker SUV may not be an optimal solution in clinical applications such as can-
cer staging and therapy assessment. Dynamic PET provides another solution for this
purpose.
1.4 WholebodyDynamicPET
Wholebody PET imaging is used clinically for tumor staging in oncology [26]. Whole-
body scans step the patient through the scanner, reconstructing one image for each bed
position and then stitching the results together to form the final ”whole body” scan. A
graphical illustration of wholebody PET/CT is shown in Fig. 1.4. In this figure the
patient is being stepped through the scanner and is in the bed position covering the
abdomen.
Wholebody PET is used extensively in the clinic for diagnosing and monitoring the
treatment response of primary and metastatic tumors. It has been shown to enhance
lesion detectability for primary malignancies and metastatic lesions [27]. It is also
shown to be of more clinical value in posttreatment evaluation in Hodgkin’s disease
8
and Non-Hodgkin’s lymphoma than CT imaging [28]. Due to the limited axial field
of view (FOV), the patient is usually stepped through the scanner gantry such that the
bodybetweenbrainandthighiscovered. AwholebodystaticPETimageisobtainedby
stitchingseveralreconstructedPETimagesofdifferentbedpositions.
Clinically, only static wholebody PET data are routinely acquired. Static PET
imagesrepresentcumulativetraceruptakeatthetimeofscanning(typically1hourpost
injection)andareofclinicalvaluefordetectionandsemiquantitativeanalysisofprimary
tumors and metastatic lesions. However, static PET images do not directly reflect the
changeoftraceruptakeovertime. Therefore,thisprotocolmayloseimportanttemporal
information that may help to characterize and differentiate different tissues and disease
processes(e.g. metastaticlesionsversusinflammation).
Figure 1.4: Wholebody PETCT scan illustration. The patient is stepped through the
scanneratmultiplebedpositions.
9
Dynamic PET monitors changes in tracer uptake over time to allow the fitting of
quantitative pharmacokinetic models. In dynamic imaging, the data is acquired con-
tinuously from the time the tracer is administered and then fit to a kinetic model [29].
These models allow us to represent tracer uptake in terms of quantitative rate parame-
terswhichhavethepotentialforimprovedtumorcharacterizationcomparedtothemore
widely used semi-quantitative measures such as the SUV [30]. It has been shown that
dynamicPETimagingcanimprovelesiondetectionsensitivityandspecificityforcertain
tumors[31][22].
Traditionally,thedynamicanalysisisperformedafterstaticimagereconstructionof
multipleframes. Theseframesofdataareobtainedbydividingthewholescanningtime
intoseveralsegmentsandthenhistogrammingdatawithineachofthemintosinograms.
Typicalprotocolsuseveryshortframesatthebeginningofscan(10-30seconds)inorder
to capture the rapid change of tracer distribution immediately after injection and longer
frames (10-30 minutes) towards the end of the scan. The PET image for each frame is
thenreconstructed. FinallythedynamicanalysisisperformedtothestackofstaticPET
images.
Anumberofchallengesariseinusingthestandardmulti-frameapproachtodynamic
PET. First of all, the choice of frame numbers and durations is usually difficult to opti-
mize. An insufficient number of frames will lose temporal information in the dynamic
process. On the other hand, too many frames will significantly reduce the counts in
each frame, potentially causing bias in estimated activity in iterative image reconstruc-
tion method [32] and much noisier images. Besides, statistical analysis of the static
images will require that the resolution of these static images match (especially when a
voxel-wise analysis is needed), because it is usually favorable to have a constant tem-
poralresolutionthroughouttime. Thisisadifficultproblemwhenpenalizedmaximum-
likelihood method (MAP) is used for image reconstruction, because the resolution is
10
data-dependent [33]. The method proposed in [34] achieves a count-independent res-
olution, however matching the resolution for all the static images by selecting differ-
ent data-dependent smoothing parameters is still time-consuming. In Chapter 2 we
will show that our direct dynamic analysis without multiple-frame image reconstruc-
tionobviatesthislimitation.
Currentwholebodyimagingprotocolsareincompatiblewithdynamicstudies,which
requires continuous data acquisition for multiple bed positions. In Chapter 3, we will
describeatwo-passdynamicimagingprotocolforwholebodyPET.Wethenreconstruct
the Patlak parametric images from the listmode data. Recently Karakatsanis et al [35]
[36]proposedaclinicalprotocolforwholebodyPatlakestimation. Theirworkalsoaims
at the quantification of tracer dynamics using Patlak model. However their method is
basedonthetraditionalmulti-passsinogramapproach. Nonethelesstheirresultsdemon-
stratedadvantageofPatlakmodelingoverSUV.
1.5 TracerKineticModeling
TracerKineticModelingisfrequentlyusedtoderivephysiologicalandpharmacological
parametersfromPETdata[37].
One of the most commonly used models is the two-tissue compartment model as
showninFig. 1.5with4kineticparameters,K
1
,k
2
,k
3
,k
4
,denotingrateoftransferfrom
plasma to first compartment, first compartment to plasma, first to second compartment,
and second compartment to the first respectively. The two-tissue compartment model
is appropriate for
18
F-FDG modeling, where the first tissue compartment represents the
FDG concentration and the second tissue compartment represents the phosphorylated
FDG (FDG-6-P) concentration in the tissue respectively. They are different in that the
tracerinthefreecompartmentisnotmetabolizedbutwillbemetabolizedinthesecond.
11
This model describes the dynamic tracer distribution in the plasma and tissue with the
assumption that the rate of transfer between each compartment is constant. Whenk
4
=
0,thetracerisreferredas”irreversible”sinceitcannotgobacktothefirstcompartment
from the second. For
18
F-FDG,k
4
is approximately 0 for most tissue types and tumors.
For the two-tissue compartment model the equations for tracer concentration in tissue,
ROI(t), are given in Eqs 1.2, 1.3, 1.4 where C
p
(t) is the concentration of tracer in
arterial plasma, C
B
(t) is the total radioactivity concentration in tissue vasculature, and
V
B
isaknownconstantreflectingvolumetricfractionofbloodintissue. C
1
(t)andC
2
(t)
aretheconcentrationsofnon-metabolizedandmetabolizedtracerintissuecompartment
1and2respectively.
Figure 1.5: 2-tissue compartment model for kinetic analysis. C
p
(t) is the activity in the
plasma. C
1
(t)andC
2
(t)correspondtotheactivityin2tissuecompartmentsattimet.
dC
1
(t)
dt
=K
1
C
p
(t)(k
2
+k
3
)C
1
(t)+k
4
C
2
(t) (1.2)
dC
2
(t)
dt
=k
3
C
1
(t)k
4
C
2
(t) (1.3)
12
ROI(t) =C
1
(t)+C
2
(t)+V
B
C
B
(t) (1.4)
The two-tissue compartment model is a quantitative analysis of the transfer rates of
tracerinandoutofthetissue.
It has been shown that kinetic analysis provides clinical benefits. Takesh [38]
showed the benefits of kinetic analysis in improving the diagnostic accuracy in vari-
oustumorsanditsapplicationontherapymonitoring. Sugawara[39]foundstatistically
significant differences in the kinetic parameters between mature teratoma and necrosis
orscar,althoughvisualinterpretationandSUVmetricisnotcapableofmakingadiffer-
entiation. HeconcludedthatFDGPETwithkineticanalysiscanbeapromisingmethod
formanagementofdiseaseinpatientswithgermcelltumoraftertreatment.
In full kinetic analysis using the two-tissue compartment model, the optimization
problemisnonconvexforthekineticparametersbecausetherateparametersK
1
,k
2
,k
3
,
k
4
are nonlinearly related to the data. Furthermore, the estimates are unstable and the
computational cost is high. Therefore, we use a simplified graphical method instead of
thefullkineticanalysisasanalternate,asintroducedinthenextsection.
Traditional dynamic PET image reconstruction histograms photon counts into mul-
tiple sinograms and then reconstructs a sequence of images. However, this histogram-
ming operation loses the temporal information of counts arriving at different time but
histogrammed into one sinogram. Furthermore, this will be a quite inefficient way to
store data especially for low count dynamic studies (often the case in modern three-
dimensional (3-D) PET systems [40]). For example, a dynamic study with N-frame
sinogramdataneedsNtimesthespaceofasinglesinogram.
It is therefore beneficial to store the data in a different format, i.e. the listmode data
format. Thelistmodefilestoreseacheventanditsinformationsuchasarrivingtimeand
the 2 corresponding detector indexes. Unlike sinogram data format that allocates space
13
for each possible LOR and for each frame, the listmode data saves the arrival events
onlyanddoesnotstorethezero-count LORs. Asma et al. [41] has shownthat the sino-
gram/timogram format (the sinogram stores the number of coincidences for each LOR
in the whole acquisition time and the timogram stores the arrival time of the events)
significantlyreducethedatastoragecomparedwithmulti-frameandcanachievehigher
temporalresolution. Meanwhile,thelistmodedatacontainsmoreinformationthansino-
gram data. We have shown from Cramer-Rao bound analysis that for direct Patlak esti-
mation,listmodemethodyieldssmallervariancethansinogrambasedmethod[42].
Considering these benefits of the listmode data format, the direct dynamic analysis
methodfromlistmodePETdataisattractive. Thisdataformatallowsanevent-by-event
calculation in parametric image reconstruction with inhomogeneous Poisson modeled
data,followingSnyder’s[43]derivationofthePoissonlikelihoodforlistmodePETdata.
Also, we can analyze the behavior of temporal and spatial resolution and the influence
of the penalty so as to enforce that the resolution of the final images is desirable [42],
which is more difficult with multiple static images. Furthermore, the estimation finally
results in only one reconstruction instead of multiple static reconstructions, thus it has
thebenefitsofcomputationalefficiencyandsimplicityinpracticalprocessing.
1.6 GraphicalAnalysisofDynamicPETData
Graphical analysis methods have been applied frequently to dynamic PET data. The
basicideaofgraphicalanalysisistotransformthemeasuredPETdatasuchthatalinear
relationshipcanbeestablishedforthetransformedvariables,wheretheslopeandinter-
cept have physiological meanings. Patlak graphical analysis was proposed by Patlak
[44],whosimplifiedthetwo-tissuecompartmentkineticanalysisforirreversibletracers
14
byfindingalinearrelationshipbetweenthetotaltissueactivityandthebloodinputfunc-
tionaftersomedatatransformationinthesteadystate. Herethetermsteadystatemeans
thatallfluxesandvelocitiesareconstantwithrespecttotime[45]. ForthecaseofFDG,
the steady state for most of the tissue types begins around 30-45 min after injection.
Rewriting ( 1.2) withk
4
= 0 yields the equation (1.5), and can be further simplified as
(1.6).
ROI(t)
C
p
(t)
=
K
1
k
3
k
2
+k
3
R
t
0
C
p
()d
C
p
(t)
+
k
2
k
2
+k
3
C
1
(t)
C
p
(t)
+
V
B
C
B
(t)
C
p
(t)
(1.5)
ROI(t)
C
p
(t)
=
R
t
0
C
p
()d
C
p
(t)
+q (1.6)
Afterthesteadystatebegins,thetissueFDGconcentrationC
1
(t)followstheplasma
concentration. This leads to C
1
(t)=C
p
(t) and C
B
(t)=C
P
(t) be constant, which ensures
that the intercept is constant and the Patlak equation is valid. The term
K
1
k
3
k
2
+k
3
is known
asthe Patlakslope andthe term
k
2
k
2
+k
3
C
1
(t)
Cp(t)
+
V
B
C
B
(t)
Cp(t)
is the Patlak intercept. The Patlak
equation provides a much simpler relation between the tissue activity and the blood
input function. Particularly, the Patlak slope is physiologically meaningful because it is
thenetinfluxrateofthetracer.
The blood input function can be obtained invasively with arterial blood sampling,
known as an intrusive and risky process [46]. It is usually difficult to perform in clinic.
Severalmethodshavebeenproposedfornoninvasiveestimationofthebloodinputfunc-
tion. Wang[47]proposedamethodtoestimatetheinputfunctionandthekineticparam-
eters simultaneously. The method first runs the analysis of space intersections to obtain
an initial estimate of the input functions, and then iteratively updates the input function
and the receptor parameters. Naganawa [48] used an intersectional searching algorithm
and obtained a robust estimation of the arterial input function. Beside the approach of
15
estimating blood input function, the reference region also provides a possible solution
to obviate the usage of blood input function. The reference region corresponds to the
tissue where the tracer is not accumulated because of specific binding or metabolism,
i.e. k
3
= 0. Thedifferentialequationcanthereforebesimplifiedas:
dREF(t)
dt
=K
1
REF
C
p
(t)k
2
REF
REF(t) (1.7)
With the integral of this equation REF(t) = K
1
REF
R
C
p
()d
k
2
REF
R
REF()d and substitute the expression of C
p
(t) into ( 1.5), we can
obtainthefollowingPatlakequationwithreferenceregions:
ROI(t)
REF(t)
=
K
1
K
1
REF
k
2
REF
k
3
k
2
+k
3
R
t
0
REF ()d
REF(t)
+
K
1
K
1
REF
k
3
k
2
+k
3
+
k
2
k
2
+k
3
C
1
(t)
REF(t)
(1.8)
The Patlak model only applies to irreversible (k
4
=0) tracers. For reversible tracers
withk
4
6=0,thetracercangobackfromthesecondtissuecompartmentbacktothefirst
and different models are needed. By solving the same differential equations 1.2, 1.3,
1.4withoutdroppingthek
4
term,Logan[49]proposedtheLogangraphicalequationas
in (1.9), which is also a linear approximation equation after the steady state, when all
fluxes and velocities becomes constant with respect to time. Here DV =
K
1
k
2
(1 +
k
3
k
4
)
is the distribution volume. The Logan method is simple and more robust than the full
kineticanalysis.
R
t
0
ROI()d
ROI(t)
=DV
R
t
0
C
p
()d
ROI(t)

1
k
2
(1+
k
3
k
4
)
1
k
4
C
2
(t)
C
1
(t)+C
2
(t)
+
V
B
R
t
0
C
B
()d
ROI(t)
(1.9)
16
Neglecting the last term on the assumption that the blood volume fraction is neg-
ligible, the condition for linearity in Logan equation is that the intercept(1=k2)(1 +
k3=k4) (1=k4)(C
2
(t)=(C
1
(t) +C
2
(t))) is constant. After the steady state, the com-
partment concentrations gradually follow the plasma concentration and the equilibrium
betweenthe2tissuecompartmentsarereachedsothatC
2
(t)=(C
1
(t)+C
2
(t))isconstant,
whichensuresthattheinterceptisconstant.
SimilarasinPatlakequation,theLoganequationcanalsousethereferenceregionto
obviatethenecessityofbloodinputfunction. Skippingthesimilarderivationasin[49],
we present the final formula of the reference region based Logan equation as follows.
Here DVR = DV=DV
REF
is the ratio of distribution volume for the tissue and the
reference region. REF(t) is the activity function for the reference region andDV
REF
is
thedistributionvolumeforthereferenceregion.
R
t
0
ROI()d
ROI(t)
=DVR
R
t
0
REF ()d +REF(t)=k
2
REF
ROI(t)
+q (1.10)
As shown above, the corresponding Patlak and Logan graphical analysis methods
both reduce to transforming the data and fitting to a linear regression model. Unlike
kineticanalysisinvolvingill-conditionednonlinearoptimizationwith3or4parameters,
the Patlak and Logan analysis only need a linear regression to estimate the slope and
intercept. So compared with the original compartmental analysis, the Patlak and Logan
analysisismorerobustandlesscomputationallyintensive. Therefore,ourgoalindirect
dynamicanalysisistoestimatePatlakandLoganparameters.
One problem for the Logan parameter estimation is that the noise in the time-
activity-curve will introduce a bias in the distribution volume (the slope), especially
for large distribution volume [50] [51]. Feng [52] [53] has proposed a generalized least
17
squaremethodfornonuniformlysampledsystemtoremovethebias,whenthedatacan
be simplified to match a one-tissue compartment model. This is potentially useful in
removingthebiasinLoganestimationfromnoisydata. Loganmethodalsodiffersfrom
Patlak method in that the data has to be reorganized so that the dependent variable rep-
resents the integration of activity from time = 0 to each time point. This requires that
thedataacquisitionstartsfromthetimeofinjection.
1.7 DirectEstimationofKineticParametricImagefrom
DynamicPETData
Many researchers have investigated direct estimation of kinetic parameters from
dynamicPETdata,asopposedtothetraditionaldynamicanalysis,whichfitsthemodel
from multiple reconstructed images. This idea was first proposed by Snyder [43] in
1984 for the listmode data (the listmode data is a list saving the events information,
storingthespatialdetectorpairindexesandthetiminginformationassociatedwitheach
arriving photon pair). Since then several direct estimation method have been described,
includingmaximumlikelihoodoriented[54],[55],[56]andpenalizedMLmethods[57],
[58].
However, the rate parameters in these kinetic models are nonlinearly related to the
data and result in nonconvex optimization problems, which can be difficult to solve.
Furthermore, the relatively high dimensionality of the search (four per voxel for the
two compartment model for reversible tracer and three for irreversible tracer such as
FDG) can result in instability of the estimates, limiting the use of these direct estima-
tors. Therefore, we perform direct Patlak and Logan analysis instead of the full kinetic
analysis.
18
For the Patlak equation, the rate function in the tissue after the steady state can be
represented as the linear combination of two basis functions, the blood input function
and its integral. These 2 basis functions are distinct because after the steady state, the
blood input function is decreasing and its integral is increasing. Therefore theoretically
using2disjointframesafterthesteadystate,weareabletoestimatethe2Patlakparam-
eters(slopeandintercept). ThisindicatesthatthePatlakanalysiscanbeperformedwith
limiteddataafterthesteadystate.
Before our work, Wang et al. [59] has described a penalized ML method for direct
voxel-wise reconstruction of Patlak parameters using a preconditioned conjugate gra-
dient method applied to a sequence of sinogram frames. The key formula relating the
observeddatatothePatlakparametersis(1.11)
y = (A
P) +r (1.11)
where y is the vector representing mean counts in sinogram bins, is a vector con-
structed by Patlak slope and intercept for each image voxels. P is the PET system
matrix, Aisthetemporalmatrix withA = [

S
p

C
p
] and

S
p
= [

S
p
(1)

S
p
(2):::

S
p
(N)]
T
and
that

C
p
= [

C
p
(1)

C
p
(2):::

C
p
(N)]
T
. r denotes the randoms and scatters.

S
p
and

C
p
are
the integral of 2 temporal basis functions at N sinogram frame time. The C
p
(t) is the
bloodinputfunctionandS
p
(t)istheintegralofthebloodinputfunction. Thissinogram
based method potentially still loses information in the data, because multiple counts
arrivingatclosetimearehistogrammedintothesamesinogram. Karakatsanisetal[60]
[36] proposed a clinical protocol for wholebody Patlak estimation composed of an ini-
tial 6min dynamic scan over the heart followed by a sequence of PET scans (6-pass x 7
bed position each). The Patlak image estimate is obtained by applying a hybrid linear
regression model to the temporal static images. Their approach uses an ordinary least
square (OLS) regression to estimate the parameters, and also aims at the quantification
19
of tracer dynamics with Patlak model, with results indicating the advantage of Patlak
modelingoverSUV.Nevertheless,againthetemporalinformationofcountsisnotfully
used. Furthermore,theresolutionandnoisepropertiesofthefinalPatlakimagesarenot
wellinvestigated.
The above two Patlak parametric estimation methods use sinogram frames, which
potentially lose information of the photon pairs arriving in the same sinogram acquisi-
tion time. To fully explore the benefits of temporal information of each arrival event,
wedevelopedanalternativepenalizedMLapproachthatuseslist-mode(List-modedata
storestheLORinformationtogetherwiththeexactarrivaltime. Thiswillbeintroduced
more in the next section) rather than binned sinogram data [61]. We use an incremen-
tal gradient algorithm originally developed for nonparametric dynamic imaging with a
B-spline temporal basis [62]. Here we investigate a modified version of this algorithm
that uses list-mode data from two acquisition frames. Preliminary results of this work
were reported in [63] and [61]. Simplifying the Patlak equation after the steady state
wouldyieldequation 1.6,whereisthePatlakslopeandqisthePatlakintercept. After
multiplying C
p
(t) at both hand sides, the activity rate at time t can be represented by:
ROI(t) =
R
t
0
C
p
()d +qC
p
(t). This indicates that the rate function after the steady
statecanberepresentedbytheweightedsumoftwobasisfunctionsthataredistincttem-
porally. As discussed before, this is a well-conditioned linear problem allowing using
partial data after the steady state to estimate the Patlak parameters. This idea finally
leads to our wholebody Patlak esitmation method, where we designed a dual-pass scan
protocolwith5bedpositionsperpass.
The direct estimation of Logan parameters is developed based on a modified ver-
sion [64] of the original Logan equation. The Logan equation is shown in ( 1.12),
where DV is the Logan slope and Int is the Logan intercept. As an approximation
to the Logan equation, Zhou [64] proposed the RE (relative equilibrium) equation
20
as in ( 1.13), which is derived based on the assumption that the activities in the 2
compartments gradually follow the plasma after sufficient time, i.e. C
1
(t)=C
p
(t) and
C
2
(t)=C
p
(t) approach constant. This formula allows direct estimation of the Logan
parameters to be possible. By multiplying C
p
(t) at both hand sides, the accumu-
lated activity from time=0min can be represented by the weighted sum of two terms:
R
t
0
ROI()d =DV
R
t
0
C
p
()d +qC
p
(t). Ifbloodinputfunctionisknown,thislinear
equation would allow direct estimation of DV and q possible from several accumulated
sinograms. However, by using accumulated sinogram to perform the estimation, this
method may potentially lose some temporal information of the data. Moreover, this
problem does not consider the correlation between the accumulated sinograms, which
mayincurbiasorextravariancetothefinalestimation.
R
t
0
ROI()d
ROI(t)
=DV
R
t
0
C
p
()d
ROI(t)
+Int (1.12)
R
t
0
ROI()d
C
p
(t)
=DV
R
t
0
C
p
()d
C
p
(t)
+q (1.13)
To obviate the data dependency problem and fully explore the temporal informa-
tion of each arrival event, an algorithm modeling the rate function of tissue activity and
performingevent-by-eventcalculationinsteadofusingsinogramsrepresentingaccumu-
latedactivitymaybeeffective. However,methodsofdirectestimationofLoganparame-
terswithratefunctionmodelingandevent-by-eventcalculationhasnotbeenestablished
before our work. This is because in the RE equation (1.13), the rate function of tissue
can be derived as ROI(t) = DVC
p
(t) +qC
0
p
(t) after the steady state. However after
the steady state,C
p
(t) usually approaches a single exponential function, i.e. The blood
input functionC
p
(t) = Ae
t
. Its derivativeC
0
p
(t) =Ae
t
, which is times the
bloodinputfunction. Thesimilarityof2basisfunctionsovertheentireacquisitiontime
21
afterthesteadystatemakestheestimationproblemquiteill-conditioned. Wedeveloped
a feasible solution to estimate the Logan parameters which allows rate function model-
ing and event-by-event computation after the steady state, by using an additional single
sinogram from time=0min to the steady state. This additional static sinogram turns out
togreatlyimprovetheill-conditioningoftheproblem. Ourmethoddoesnotsufferfrom
correlatednoisebecauseitdoesnottakemultipleaccumulatedsinograms.
1.8 PETattenuationcorrectionwithcontrastCT
AttenuationinPETimagingisthelossoftruecoincidenceeventsbecauseofabsorption
ofphotonsinthebodyortheComptonscatter. Thelossoftrueeventsduetoattenuation
in clinical PET can be up to 95%, especially for large patients. Without applying atten-
uation correction, image reconstruction may suffer from artifacts or distortion, such as
showingprominentactivityatbodysurfaceandlowactivityindeeperstructures. There-
fore, attenuation correction is critical for quantitative PET image reconstruction and
analysis.
Therearemainly2approachestoobtaintheattenuation. Beforethedevelopmentof
PET/CT scanners, the attenuation was obtained with a transmission source (e.g.
68
Ge)
directly measuring the photon attenuation at 511keV. The additional transmission scan
usually took 3-10 min or more for each bed position. This significantly increases the
time required for clinical whole-body procotols (typically 5-7 bed positions) and limits
the usage of the technique [65]. Another solution is to perform simultaneous PET and
CT(PET/CT)andthenmapthecoregisteredCTimagetothecorrespondingattenuation
coefficients with a pre-defined mapping table. Compared with the transmission source
method, the CT-based attenuation estimation is much faster because of the short acqui-
sition time of CT scans. In the image quantification aspect, transmission source scans
22
havinghighernoiseandlowerbiasbutwillsufferfromextrabiasfromthepost-injection
transmissionscan,andCTscanshavenegligiblenoisebutpotentialquantificationerrors
[66]andmaycauseartifactsduetothemismatchofPETandCTimagescausedbyrespi-
ratoryandcardiacmotion. Burger[67]andKinahan[68]havecomparedtheattenuation
obtained from both methods and found that the attenuation obtained by these 2 meth-
ods are generally equivalent. In clinical settings the CT based attenuation estimation is
widelyusedincombinedPET/CTscans[68].
TheCTHounsfieldunitreflectshoweasilytheX-raycanpenetratethroughthemat-
ter. The Hounsfield unit for water is defined to be zero in a conventional CT scan.
The bone has the highest Hounsfield unit (about 1000 HU) and the air has the lowest
Hounsfield unit (-1000HU). The conventional CT scans has poor soft tissue contrast,
because the Hounsfield unit of the soft tissues is usually around 30-70 HU such as
spleen, kidney, and liver. Recent PETCT scans are sometimes performed with con-
trastenhancedagentinCTsothatthestructureofdifferentorgansandtissuesespecially
the soft tissues will have a better contrast from the background and show more details.
These soft tissues will have an enhanced Hounsfield unit between 100 - 300 HU which
isquitedifferentfromtheirnormalCTvalues.
The contrast enhance CT is shown to be more useful than conventional CT in some
applications. Antoch[69]concludedthatthecontrastagentincontrastCTcanaddvalue
to lesion detection and characterization. Fletcher [70] reported that contrast enhanced
CT colonography is a promising method for detecting local recurrence, metachronous
disease, and distant metastases in patients with prior invasive colorectal carcinoma.
Oliver [71] found that in the application of multiple phase contrast CT, using arterial
phase and portal venous phase images is better than using the conventional and portal
venous phase images to reveal more hepatocellular carcinoma lesions. However, it is
23
necessary to point out that there are also some researchers claiming the opposite find-
ings. For example, Schaefer [72] reported that noncontrast CT is more sensitive and
specific than contrast-enhanced CT for evaluation of lymph node and organ involve-
ment, especially regarding exclusion of disease, in patients with Hodgkin disease and
high-gradenon-Hodgkinlymphomathan.
Thecontrastagentcausespotentialproblemforaccurateattenuationcorrection. The
contrastagentincreasestheHounsfieldunitoftheCTimageashighasseveralhundred
in some blood rich regions such as the spleen or aorta. However when the PET scan
is performed, the attenuation of the contrast agent is similar to the water in the
-ray
energy range. Therefore using the predefined table, also known as the CT-to-mumap
mappingtable(Theattenuationmapiscalledmumap. Itrepresentsthelinearattenuation
coefficient. For water in 511keV it is 0.096 cm
1
) for the contrast enhanced CT will
overestimate the attenuation coefficients for the contrast enhanced regions. Beside, the
attenuation coefficient for each LOR is not a linear function of the line integral of the
CT-converted mumap. Instead, it is an exponential relation, which further complicates
the problem. Many researchers reported overestimated SUV or attenuation coefficients
whencontrastagentisused[73],[74],althoughYau[75]reportedthatcontrastenhanced
CTdoesnotsignificantlyincreasetheSUV
max
.
There are generally two main directions in PET attenuation estimation when the
conventional CT and the transmission scan are not available. One main direction is
the segmentation-driven methods. These methods basically apply segmentation algo-
rithms to segment regions with constant attenuation and assign predefined coefficients
[76], [77]. Currently the segmentation based approach is widely used for attenuation
estimation in combined PETMR machine, by first segmenting the MR image and then
usingpriorinformationtoassignmumapcoefficientstodifferentregions. However,this
method does not reflect the individual character, and may incur large error for specific
24
patients, for example, when the patient has metal implant. The accuracy of segmen-
tation and registration also affect the quantitative results. The other direction is joint
estimation of activity and attenuation (also known as MLAA) from PET emission data
[78, 79, 80]. The theoretical foundation is the analytical consistency condition and
the discrete consistency condition. However, this joint estimation in the maximum-
likelihood or penalized maximum-likelihood framework is a nonconvex problem, and
results in ”cross-talk” between activity and attenuation, or a local maxima of the cost
function. This can be explained by the fact that the Fisher information matrix for non-
TOF PET data is ill-conditioned so there is no global maxima [81]. With some prior
information,however,theconditioningoftheproblemcanbeimprovedtoproduceuse-
ful local maxima that may be close to the global maxima [82, 83]. In particular, Nuyts
et al [83] proposed a jointly iterative algorithm for estimating emission and attenuation
simultaneously as shown in equation 1.14 and 1.15. In each iteration, equation 1.14
is updating the emission while fixing attenuation (normal MLEM) and equation 1.15
is updating the attenuation mumap while fixing the emission estimation. The prior that
the intensity of attenuation map should be centered at several narrow peaks indicating
normaltissueorairorboneismodeledintotheoptimizationprocedureprovidingalocal
optimalsolutionthatisusefulinpractice.

new
j
=

j
P
i
a
i
p
ij
X
i
p
ij
y
i
b
i
(1.14)
u
new
j
=u
j
+
P
i
p
ij
a
i
b
i

P
i
p
ij
y
i
+M
0
j
(u)
N
P
i
p
ij
a
i
b
i
M
00
j
(u)
(1.15)
Time-of-Fight (TOF) PET data provides more information than conventional
nonTOF PET data. Interestingly, Defrise showed that in time-of-fight PET, the atten-
uation sinogram is determined by the emission data except for a constant [84]. Rezaei
25
alsostudiedthejointestimationprobleminTime-of-FlightPETandtheoreticallyproved
theuniquenessofglobalmaximauptoaconstantfortheattenuation[81]. Thisproperty
of TOFPET makes MLAA expecially useful in many applications. Boellaard showed
that MLAA can improve quantitative accuracy of PET image in PET/MR studies [85].
HamillfoundthatTOF-MLAAreducedmotion-relatedPETimageartifactsinthelower
lungsinaphantomstudy[86].
Considering the large amount of nonTOF PET done with contrast enhanced CT, it
is quite useful to develop a method to estimate the attenuation accurately. We take the
approachofjointestimationforemissionandattenuation. Takingthecontrastenhanced
CT as the prior is a possible solution to improve the ill-conditioning of the joint esti-
mation. As the contrast agent is usually concentrated in the blood-rich soft tissues such
as liver, heart and spleen and rarely exists in low-CT intensity regions such as lungs or
high-CT intensity regions such as bone. The linear attenuation coefficients of a large
portion of the voxels can be accurately determined from the corresponding contrast CT
iamge. For any specific pixel with Hounsfield unit in the contrast-affected regions, the
contrastCTalsoprovidesusefulinformation. Itmaybelongtoboneregionswithpartial
volumeeffectthusshowslowerintensitythanpurebone(inwhichcaseweshouldmain-
tain the current CT value), or belong to the contrast affected soft tissue (in which case
we have to map it to the ordinary soft tissue attenuation coefficients instead of treating
it as noncontrast region and use the conventional mapping). Therefore although a strict
global optimal solution is difficult, we can use this prior information to obtain a useful
solutionthatiscloseenoughtotheglobaloptimalsolution.
26
1.9 Contributions
One major contribution of this dissertation is the development of a clinically feasible
wholebodyPatlakimageestimationmethod. Ourmethodisalsodifferentfromthetradi-
tionalapproachbyreconstructingawholebodyPatlakimagedirectlyfromlistmodePET
data. WefurtherappliedthealgorithmtoclinicaldataandcomparedPatlakimageswith
SUV images in a pilot study of 10 patients. Various practical issues such as dynamic
modeling of randoms and scatters, blood input function estimation, and motion correc-
tion are discussed. In both the simulated and real data analysis, we found superior per-
formance of the Patlak parameters compared with the SUV based quantities, including
smaller inter-subject variation and high tumor-to-background contrast. Compared with
otherPatlakestimationmethodsinliterature,ourworkdistinguishesitselfbyproposing
a clinically feasible and efficient method to obtain wholebody Patlak image allowing
voxelwiseanalysis.
Another contribution of this dissertation is the extension of our direct parameter
estimation method to the reversible tracers, where Logan analysis is usually used. We
showthatthereversiblenatureofthetracerfromthespecific-bindingcompartmenttothe
free compartment leads to a more ill-conditioned problem for direct estimation if data
is only collected after the steady state. We solve this ill-conditioned problem by adding
a static sinogram storing counts from injection to the steady state. Our method differs
from the sinogram-based estimation method proposed by other researchers in that we
develop a correct likelihood function reflecting the true characteristic of the data, while
the sinogram-based approach suffers from correlated noise between the accumulated
sinogram data. Simulation results show that the variance obtained by our method is
closer to the Cramer Rao bound than the sinogram-based approach. This is the first
feasiblelistmode-drivenestimationalgorithmfordirectLoganparameterestimation.
27
1.10 Organizationofthedissertation
Thisdissertationisorganizedasfollows: chapter2presentsthelistmodereconstruction
algorithm for our Patlak estimation problem. We introduced our model and algorithm,
the Cramer-Rao bound analysis and the techniques to achieve desired resolution prop-
erties in this chapter. Simulated ROC study is then performed to compare the Patlak
slopeandfractionalSUV.Chapter3concentratesontheapplicationtoclinical
18
F-FDG
studies. Wewilldescribetheestimationmethodfortheratefunctionwithdynamicscat-
terandrandommodeling,andanon-invasivemethodtoderivethebloodinputfunction.
Wewillalsodemonstratethatourmotioncorrectiontechniqueimprovestheimagequal-
ity. In Chapter 4 we introduce the direct Logan parameter estimation by modifying the
original listmode reconstruction algorithm to improve the ill-conditioning, and further
verify the effectiveness with Cramer-Rao analysis and Monte-Carlo simulation. Chap-
ter5describesthejointestimationalgorithmforemissionandattenuationwithcontrast
enhanced CT. The comparison with simple rescaling scheme shows that our method
with function designed using the contrast CT image can greatly improve the estima-
tion accuracy of emission and attenuation. We conclude our work with future research
directionsinChapter6.
28
Chapter2
DirectPatlakestimationfromlistmode
PETdata
2.1 Introduction
Positron emission tomography (PET) is a powerful technique for staging cancer and
assessing response to therapy. Whole body FDG (18F-Fludeoxyglucose) imaging is
among the most common clinical protocols in which the patient is stepped through the
scanner so that a 3D volumetric image of the whole body can be reconstructed. These
static images represent accumulated tracer uptake at the time of scanning (typically 1
hourpostinjection)andareofclinicalvaluefordetectionandsemiquantitativeanalysis
ofprimarytumorsandmetastaticlesions. However,staticPETimagingdoesnotexploit
the full potential of this inherently dynamic modality. An alternative approach is to
acquiredatacontinuouslyfromthetimethetracerisadministeredandtothenfitakinetic
model [29]. These models allow us to represent tracer uptake in terms of quantitative
rateparameterswhichhavethepotentialforimprovedtumorcharacterizationcompared
to more widely used semi-quantitative measures such as the SUV (standardized uptake
value)[30].
Standard whole body and dynamic imaging protocols are inconsistent with each
other since one requires imaging at multiple positions while the other requires acquisi-
tion of complete dynamic data at a single position. In this paper we develop an image
estimation approach that could be used for dynamic whole body imaging. To achieve
29
thisgoalweusethePatlakmethod[45],asimplifiedkineticmodelforirreversiblybound
tracerssuchasFDG,andweestimatetheparametersofthismodelfrompartialdynamic
data consisting of short frames of data acquired over two disjoint observation periods.
Wholebodydatacanthenbeacquiredattwodifferenttimepointsforeachbedposition
and the Patlak model fitted separately for each bed position. The resulting images can
thenbestitchedtogethertoformaquantitativewholebodydynamicimage.
Dynamic models are frequently fit to a sequence of reconstructed images with one
setofparameterscomputedfordataaveragedoverauserselectedregionofinterest. An
alternative approach, and that adopted here, is to directly fit the data to a kinetic model
at each voxel in the image. This idea was first proposed by Snyder [43] in 1984. Since
then several maximum likelihood [54], [55], [56] and penalized ML [57], [58] meth-
ods have been described for this purpose. The rate parameters in these kinetic models
are nonlinearly related to the data and result in nonconvex optimization problems. Fur-
thermore, the relatively high dimensionality of the search (four per voxel for the two
compartment model typically used for FDG) can result in instability of the estimates
andhighcomputationcost,limitingtheuseofthesedirectestimators.
The simplified Patlak model, applicable when tracers are irreversibly trapped in the
second compartment (k
4
= 0), uses an approximation that models the measured time
activity curve in steady state as a weighted sum of the input function and its integral.
The weights are referred to respectively as the Patlak slope and intercept. The slope in
particular is a physiologically meaningful quantity representing the net transfer rate or
influx constant and could be used in place of SUV values as a biomarker. For example,
Graham et al. [31] found that Patlak slope is a better predictor of outcome and a better
discriminator between normal tissue and tumor than SUV after comparing them in a
groupof40patientswithcoloncancermetastasizedtotheliver. WenotethatthePatlak
approximation is only valid for irreversible trapping, which is not always the case for
30
FDG. For example, Iozzo et al [87] show nonzero k
4
in liver. However, in practice the
Patlakmodelisfrequentlyandeffectivelyusedinstudiesoftumorviability.
Wang et al. [59] described a penalized ML method for voxel-wise reconstruction
of Patlak parameters using a preconditioned conjugate gradient method applied to a
sequence of frames. We developed an alternative penalized ML approach that uses list-
mode rather than binned data [61]. In that case, as in the current thesis, we used an
incremental gradient algorithm originally developed for nonparametric dynamic imag-
ing with a B-spline temporal basis [62]. Here we investigate a modified version of this
algorithmthatuseslist-modedatafromtwoacquisitionframes. Preliminaryversionsof
thisworkwerereportedin[63]and[61].
The concept of using two-pass whole body scans is not new, but has been exten-
sivelystudiedinthecontextofSUVvalues. ThesestudieslookatchangesinSUVvalue
betweentwoacquisitionframes,typicallyexpressedaspercentagefractionaldifference:
%DSUV = 100(SUV
2
SUV
1
)=SUV
1
%. By using fractional changes, these mea-
sures become invariant to the scaling by dose and body-weight used in standard SUV.
For this reason %DSUV reflects only dynamic changes in uptake and is robust to the
uncertainties governing the relationship between FDG uptake and patient weight that
make standard SUV values only semi-quantitative. Alkhawaldeh et al. [20] found that
dual time-point PET can improve diagnostic accuracy relative to standard single frame
SUV,increasingsensitivityandspecificity formalignant lung nodules diagnostic, espe-
cially for small lung lesions that have low SUVs. Prieto et al. [21] showed that dual
time-point PET can improve sensitivity for the identification and volume delineation of
high-grade brain tumors compared with standard PET studies. Hu et al. [88] reported
thatformediastinalnodalstaging,thespecificity,accuracy,andpositivepredictivevalue
of dual time-point scans were better than those of single-time-point by adding extra
%DSUV information. Zhuang et al. [22] and Kumar et al. [23] found that malignant
31
lesions showed a significant increase in SUV over time, while benign lesions showed a
decrease. As result, %DSUV can have advantages over SUV in differentiating malig-
nant from benign lesions. However, negative results have also been reported: Cloran et
al. [24]reportedthatdualtime-pointFDGPET maynotbe ofbenefit intheassessment
of pulmonary nodules with maximum SUV of less than 2.5 on initial imaging Chan et
al. [25] suggested that dual time-point FDG-PET does not improve the overall eval-
uation of pulmonary lesions. In the studies presented below we perform comparisons
between our Patlak based method and %DSUV since both can make use of identical
dualtime-pointdata.
Recently, Karakatsanis et al [60] [36] proposed a clinical protocol for wholebody
Patlak estimation composed of an initial 6min dynamic scan over the heart followed by
a sequence of PET scans (6-pass x 7 bed position each). The Patlak image estimate
is obtained by applying a hybrid linear regression model to the temporal static images.
Their work also aims at the quantification of tracer dynamics with Patlak model, with
results indicating the advantage of Patlak modeling over SUV. This differs from our
approach, both in the use of multiple bed frames per bed position and in the estimation
of Patlak images from reconstructed static images. The use of multiple frames per bed
positionmaylimitthepracticalutilityofthisapproach.
This chapter is organized as follows: we first introduce the list-model Patlak model
and image reconstruction methods. We then derive an expression for the Cramer Rao
boundsonthePatlakparameters. Theproblemofensuringcount-independentresolution
inthepenalizedMLimagesisaddressedafterwards. ResultspresentedincludeCramer-
Rao analysis and lesion detectability studies for a simplified single detector pair/single
voxel model. We then present a volumetric imaging simulation and comparison with
%DSUVanddemonstratetheperformanceofourcount-independentresolutionmethod.
32
We conclude with a discussion of the performance of the method and its extensions to
processingofclinicalwholebodydata.
2.2 Method
In this section we first formulate our Patlak approach using penalized maximum likeli-
hoodestimationfromlist-modePETdata. WethenderivetheCramer-Raoboundforan
unbiased estimator of the Patlak parameters, which provides a quick means to investi-
gatelowerboundsonthevarianceoftheestimatorasafunctionoftheparametersofthe
acquisitionprotocol. Usingthisresultwecancomparethebest-caseperformanceofour
list-mode algorithm with alternatives including sinogram-based Patlak estimation and
lesion detection using %DSUV. Finally we derive a spatially-variant weighting penalty
that enforces count independent resolution in the reconstructed Patlak images. Use of
this weighting ensures that the resolution of the Patlak slope and intercept are matched,
andalsothatresolutiondoesnotvaryacrosstheimageduetospatialvariationinactivity.
Usingthealgorithmandanalysisdevelopedinthissectionwearethenabletoinvestigate
theperformanceofdualtime-pointPatlakestimationunderarangeofconditionsandin
comparisontothealternativeapproachesofsinogram-basedestimationand%DSUV.
2.2.1 DirectPatlakEstimationfromtwoFrameList-ModeData
The Patlak model provides an approximate solution to the two compartment kinetic
model for irreversible tracers in the steady state period t T
0
in which changes are
effectively due to irreversible trapping in a single compartment. Let (t) be the tracer
time activity curve (TAC) with blood input function C(t). We can write the Patlak
equation[59],[63]as:
33
(t) =
Z
t
0
C()d +qC(t) (2.1)
where is the net influx rate or slope parameter, andq is the intercept of the Patlak
model. Our goal is to compute estimates of and q at each voxel in the source space
fromlistmodedata.
Using the Patlak model, we can represent the rate function at voxel j after steady
statetT
0
asalinearcombinationoftwobasisfunctionsB
1
(t)andB
2
(t):

j
(t) =
2
X
l=1
!
jl
B
l
(t) (2.2)
where!
j1
=
j
; !
j2
=q
j
; B
1
(t) =
R
t
0
C()d; B
2
(t) =C(t)
WenotethatsincethebasisfunctionB
1
(t)istheintegralofB
2
(t),thetwofunctions
will always be linearly independent. Consequently the rate function in the measured
sinogramspaceatlineofresponse(LOR) icanbewrittenas:

i
(t) =e
t=
nv
X
j=1
2
X
l=1
p
ij
!
jl
B
l
(t)+r
i
(t)+s
i
(t) (2.3)
wheretheexponentialtermaccountsforradioactivedecayofthetracerwithhalf-life
Ln(2);p
ij
istheprobabilityofaneventatvoxelj beingdetectedatdetectorpairi;and
n
v
isthenumberofvoxels. Thisformulationisthesameasthatusedinourearlierwork
on continuous time reconstruction from list-mode data [40] except that here we use the
input function and its integral rather than cubic b-splines as the temporal basis. In the
following we use a factored system model for p
ij
based on a product of attenuation
factors, normalization, geometric response, and a Monte Carlo model of the detector
response as described in [9] . The additional terms r
i
(t) and s
i
(t) in equation (2.3)
denote, respectively, the randoms and scattered coincidence rate functions. We assume
34
that random and scatter rate functions are estimated prior to image reconstruction and
that they are separable into spatial and temporal factors as described in [40], and will
describehowwemodeltheminthenextchapter.
Assuming we have continuous list mode data over time interval [T
0
;T] and the
arrival times in the list mode data follow an inhomogeneous Poisson model, the con-
tinuoustimelog-likelihoodfunctionofeventarrivaltimesisgivenby:
L(W) =
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)+
np
X
i=1
Z
T
T
0

i
(t)dt (2.4)
where a
ik
denotes the arrival time of the kth photon at detector pair i, x
i
is the
number of events detected in LOR i, and n
p
is the total number of LORs. In the case
whenwecollectdataovertwosubintervals [t
1
;t
2
]and [t
3
;t
4
],
L(W) =
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)+
np
X
i=1
(
Z
t
2
t
1

i
(t)dt+
Z
t
4
t
3

i
(t)dt) (2.5)
where the a
ik
are constrained to those events detected in the intervals [t
1
;t
2
] and
[t
3
;t
4
]. In principal we can use this formulation to investigate behavior for a single
frameorforanarbitrarilylargenumber. However,ourgoalinthisworkistofocusona
practical approach to whole body dynamic imaging so we consider only the two frame
case. Two frame data could be acquired using similar protocols now used for %DSUV
methodsinwhichthepatientisscannedtwiceineachbedposition.
To compute images of the Patlak slope and intercept parameters we regularize the
log-likelihoodfunctionwithaquadraticspatialpenaltyfunction:
F(W) =
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)+
np
X
i=1
(
Z
t
2
t
1

i
(t)dt+
Z
t
4
t
3

i
(t)dt)
+1=2!
T
R
S
!
(2.6)
35
wherew =w
jl
isa2n
v
lengthvectorcontainingthePatlakslopeandinterceptimage
values for each voxel;w
jl
is indexed by voxel indexj = 1;:::n
v
, and basisl = 1;2 (for
slopeandinterceptrespectively). R
S
computesthesquaredpairwisedifferencebetween
eachvoxelandits26nearestneighborsforboththeslopeandinterceptparameters.
Thespatialpenaltyforourproblemcanbewrittenas
1=2!
T
R
S
! =
2
X
l=1
nv
X
j=1
X
j
0
>j;j
0
2N
j

jl;j
0
l
(!
jl
!
j
0
l
)
2
(2.7)
where the weights
jl;j
0
l
can either (i) define a shift invariant penalty by setting
them equal to the inverse of the Euclidean distance between the two voxels, or (ii) be
precomputed from the data to ensure count independent resolution in the Patlak slope
andinterceptimageaswedescribebelow.
This formulation is very similar to that described in our work on list-mode based
reconstruction of continuous time dynamic images, with the main difference being the
use of the Patlak basis functions in place of a cubic b-spline basis. Consequently, we
can compute solutions using the same numerical optimization procedure. Here we use
the 4D incremental gradient method we describe in [62] . This algorithm iterates over
subsetsofdatainsuchawayastoensureconvergencetoamaximizeofF(W)inequation
(2.6). Theupdateequationforthisalgorithmhastheform:
!
n;m
=P

[!
n;m1
+
n
D(!
n;m1
)rf
m
(!
n;m1
)] (2.8)
where n denotes the iteration number, m denotes the subiteration number, f
m
denotes the mth sub-objective function, and D(!
n;m1
) is a diagonal preconditioning
matrix, and P

is a projection operator used to ensure convergence. Unlike for the full
parametric model, the Patlak formulation results in a concave cost function so that any
36
local maximum is a global maximum. We refer readers to the details in [62] for a com-
pletedescriptionofthealgorithmandupdateequations.
2.2.2 TheCramer-Raoboundofthealgorithm
To gain insight into the potential performance of this approach from the perspective
of estimating Patlak parameters from limited data, we computed Cramer Rao lower
bounds (CRLBs) on these parameters. While we give the general expression for the
CRLBs below, in result section we explore the bounds only for the simplified case of a
single LOR and a single voxel. In that case, a spatial penalty is not required to obtain
a stable estimate so that we might expect that bounds computed from the un-penalized
log likelihood will reflect the actual variances seen in our estimators. We will examine
the bias and variance of the estimator further in the result section using Monte Carlo
simulation.
TheCRLBprovidesalowerboundonthevarianceofanyunbiasedestimatorofthe
Patlak parameters. The bounds are found from the inverse of the Fisher information
matrix(FIM)definedas
F
mj;nl
=E[
@L
2
(W)
@w
mj
@w
nl
] (2.9)
where m and n are voxel indices, and j and l are the indices of the basis functions
(j,l=1,2). Ratherthandirectlyderivetheclosed-formCRLBforthelist-modelikelihood
(2.5),aswedoin[61],weinsteadusethebin-modelikelihoodapproximationdescribed
in [89]. In this approximation we assume data are binned into multiple sinograms so
thatinthelimitastheframedurationgoestozero,thesinogramlikelihoodconvergesto
that of the list-mode data. In the closely related problem of direct estimation of Patlak
parametersfromsinogramPETdata[59],thelog-likelihoodfunctioncanbewrittenas:
37
L(W) =
X
np
i=1
y
i
log y
i
+
X
np
i=1
y
i
(2.10)
where y
i
=
R
T
T
0

i
(t)dt is the mean number of events over the sinogram,
i
(t) is the
ratefunctionasin(2.3),andy
i
isthenumberofdetectedeventsbetweentheithdetector
pair. ThepartialderivativesofL(W)are
@L
2
(W)
@w
mj
@w
nl
=
X
np
i=1
y
i
p
im
p
in
R
T
T
0
~
B
j
(t)dt
R
T
T
0
~
B
l
(t)dt
(
R
T
T
0

i
(t)dt)
2
(2.11)
where we write
~
B
j
(t) and
~
B
l
(t) as decay corrected basis functions, i.e.
~
B
j
(t) =
e
t

B
j
(t) and
~
B
l
(t) = e
t

B
l
(t). Taking the expectation we obtain the single frame
sinogramFIM:
E[
@L
2
(W)
@w
mj
@w
nl
] =
X
np
i=1
p
im
p
in
R
T
T
0
~
B
j
(t)dt
R
T
T
0
~
B
l
(t)dt
R
T
T
0

i
(t)dt
(2.12)
TofindtheFIMforsingleframelist-modedataweusethebin-modeapproximation
ofthelog-likelihoodfunction(2.4)byassumingthedataarebinnedintoalargenumber
ofsinogramsn
t
eachofdurationt:
L(W) =
X
np
i=1
X
nt
n=1
y
i
(n)
log(
Z
tn
t
n1

i
(t)dt)+
X
np
i=1
(
Z
T
T
0

i
(t)dt) (2.13)
where
i
(t) is the rate function as in (2.3) and n
t
is the number of time bins with
t
1
=T
0
andt
nt
=T. TakingpartialderivativesofL(W)leadsto
@L
2
(W)
@w
mj
@w
nl
=
X
np
i=1
X
nt
n=1
y
i
(n)
p
im
p
in
R
tn
t
n1
~
B
j
(t)dt
R
tn
t
n1
~
B
l
(t)dt
(
R
tn
t
n1

i
(t)dt)
2
(2.14)
38
Takingtheexpectationgives:
E[
@L
2
(W)
@w
mj
@w
nl
] =
X
np
i=1
X
nt
n=1
p
im
p
in
R
tn
t
n1
~
B
j
(t)dt
R
tn
t
n1
~
B
l
(t)dt
R
tn
t
n1

i
(t)dt
(2.15)
Finally,inthelimitast! 0,wehave
E[
@L
2
(W)
@w
mj
@w
nl
] =
X
np
i=1
p
im
p
in
Z
T
T
0
~
B
j
(t)
~
B
l
(t)

i
(t)
dt (2.16)
Thisresultextendsdirectlytothetwoframecaseas:
E[
@L
2
(W)
@w
mj
@w
nl
] =
X
np
i=1
p
im
p
in
(
Z
t
2
t
1
~
B
j
(t)
~
B
l
(t)

i
(t)
dt+
Z
t
4
t
3
~
B
j
(t)
~
B
l
(t)

i
(t)
dt) (2.17)
We have compared numerical values computed using the approximate finite sum
over n = 1;:::n
t
in (2.15) to those found using the more complex closed-form expres-
sion in [61] and found negligible differences for values ofn
t
> 100 We use the expres-
sion(2.15)tocomputetheCRLBsforthestudiesdescribedlater.
2.2.3 EnforcingCountIndependentResolutioninPatlakImages
It is well known that using a spatially invariant penalty function with a Poisson like-
lihood leads to a count-dependent image resolution [90]. This is because noise vari-
ance scales with mean number of counts causing the smoothing effect of the penalty
function to also be count dependent. This issue has previously been addressed using a
spatially variant weighting function for static reconstruction [90], [33], [34]. Asma and
Leahy [89] extended this approach to ensure count and temporally invariant resolution
in reconstruction of dynamic images. Here we describe a modification of the approach
39
in [91] to develop a spatially variant weighting function to ensure that resolution in the
Patlakimagesiscountindependent.
Aswedescribebelow,resolutionisafunctionoftheFIM.Weagainusethebin-mode
approximation of the list-mode likelihood from which we then derive an approximate
expression for the local impulse response of the penalized maximum likelihood esti-
matesofthePatlakimages. Thisapproximateexpressionisafunctionoftheproductof
the Fisher information and the weights
jl;j
0
l
in (2.7) and we can therefore locally vary
theseweightstoensurecountindependentresolution.
The local impulse response (LIR) at voxel j corresponds to the mean change in
the reconstructed image in response to a perturbation of the true source distribution
at that pixel. In the case of linear estimators this corresponds to the standard source-
independent impulse response of the system. However, for nonlinear estimators the
LIR is dependent on the underlying source distribution. As shown in [90] and further
investigated in [33], the LIR at voxel j and the covariance matrix for penalized ML
estimationoftheform(2.6)canbeapproximatedrespectivelyas:
LIR
j
(~ !) = [F +R
s
]
1
Fe
j
(2.18)
cov(~ !) = [F +R
s
]
1
F[F +R
s
]
1
(2.19)
Where ~ ! is the estimated Patlak image,F is the FIM, e
j
is the unit vector at voxel ,
andR
s
isthequadraticpenaltyfunctionwith aglobal(scalar)smoothingparameter.
Following [91] , the bin-mode approximation of the FIM can be written in matrix
formas:
F = (P
B)
T
Df1= y
(n)
i
g(P
B) (2.20)
40
where
denotes the Kronecker product, Df1= y
(n)
i
g is a diagonal matrix with ele-
ment (i
n
;i
n
) equal to the reciprocal of the mean of the data at LOR i in time bin n,
P is the system matrix introduced in (2.3), and B is a decay corrected temporal basis
function matrix of dimension n
t
by 2 whose (n;l)th element is given by
R
t
n
t
n1
~
B
l
(t)dt.
Noteherethematrixrepresentationassumesafinitenumberoftimebins.
To compute the LIR based on (2.18) we need to find a tractable approximation to
theinverse of the FIM. Webeginwith the approximation originally due to Fessler etal.
[90]andusedfordynamicreconstructionin[91]:
F D

D
v
1
[(P
T
P)
(B
T
B)]D
v
1
D

(2.21)
whereD

andD
v
arediagonalmatriceswithelements:

jl
=
s
X
n
d
i=1
p
ij
2
Z
T
T
0
~
B
2
l
(t)
i(t)
dt (2.22)
v
jl
=
s
X
n
d
i=1
p
ij
2
Z
T
T
0
~
B
2
l
(t)dt (2.23)
Note that in writing the integrals in (2.22) and (2.23) we are taking the limit as
t! 0,inthebinmodelikelihood,aswedoin(2.16). Similarly,inthelimit,thematrix
B
T
B hastheform:
0
@
(
R
T
T
0
~
B
1
(t)dt)
2
R
T
T
0
~
B
1
(t)dt
R
T
T
0
~
B
2
(t)dt
R
T
T
0
~
B
1
(t)dt
R
T
T
0
~
B
2
(t)dt (
R
T
T
0
~
B
2
(t)dt)
2
1
A
(2.24)
ThisapproximationisbasedontheobservationthattheFisherinformationmatrixis
concentratedalongitsdiagonalsothatareasonableapproximationshouldatleastmatch
the diagonal elements. Substituting (2.22) and (2.23) in (2.21), it is straightforward to
showthatthediagonalvaluesof(2.21)matchthoseof(2.20).
41
Denoting D

= D

D
v
1
, we can plug the FIM representation (2.21) directly into
(2.18)and(2.19)toobtaintheexpressionfortheLIRateachvoxel:
LIR
jl
(~ !) = [D

(P
T
P)
(B
T
B)D

+R
s
]
1
[D

(P
T
P)
(B
T
B)D

e
j
(2.25)
Wecanfurtherrewritethisas
LIR
jl
(~ !) =
jl
D

1
[(P
T
P)
(B
T
B)+D

1
R
s
D

1
]
1
(P
T
P)
(B
T
B)e
j
(2.26)
Similarly,thecovarianceis:
cov(~ !) =D

1
[(P
T
P)
(B
T
B)+D

1
R
s
D

1
]
1
(P
T
P)
(B
T
B)[(P
T
P)
(B
T
B)+D

1
R
s
D

1
]
1
D

1
(2.27)
While(2.26)canbeusedtoinvestigatethepropertiesofimagesreconstructedusing
the method described above, here we are specifically concerned with modifying the
penaltyfunction(2.7)soastoensureapproximatelycountinvariantresolution.
TheLIRisbyconstructionameasureoflocalblurringaroundvoxelj. Thefactorsin
(2.22)and(2.23)tendtovaryslowlyinthevicinityofanyvoxelj sothatwecanlocally
approximate
jl
D

1
= I. The matrix P
T
P associated with combined forward and
backprojection is approximately locally spatially shift invariant, andB
T
B is dependent
on the temporal bases only and therefore strictly spatially shift invariant. Consequently
the matrix P
T
P
B
T
B can be regarded as approximately spatially shift invariant. It
follows that to achieve approximate count and shift invariance in the LIR it is sufficient
42
for the term D

1
R
s
D

1
in (2.26) to be independent of the true source distribution
andspatiallyshiftinvariant.
Supposeweselecttheweights
jl;j
0
l
in(2.7)suchthatwecanwrite

R
s
=D

R
s
D

(2.28)
whereR
s
isaspatiallyinvariantoperator. Thenwewouldhave:
LIR
jl
(~ !) = [(P
T
P)
(B
T
B)+

R
s
]
1
(P
T
P)
(B
T
B)e
jl
(2.29)
cov(~ !) =D

1
[(P
T
P)
(B
T
B)+

R
s
]
1
(P
T
P)
(B
T
B)[(P
T
P)
(B
T
B)+

R
s
]
1
D

1
(2.30)
which gives us the desired count and shift invariant LIR for the Patlak slope and
intercept images. We would also like to achieve matched resolution between the Patlak
slope and intercept. With reference to (2.29) it follows that this can be achieved with a
shiftinvariantmatrixB
T
B. Wecanachievethisbyscalingthebasisfunction
~
B
2
(t)by
=
R
T
T
0
~
B
1
(t)dt
R
T
T
0
~
B
2
(t)dt
(2.31)
Withthisscalingthediagonalelements in(2.24)are equal so thatthe LIR will have
the same height for both slope and intercept images. After reconstruction the intercept
imageisthenmultipliedby 1= tofindthecorrectlyscaledvalues.
Itnowremainstodeterminetheweights
jl;j
0
l
suchthat

R
s
hasthedesiredinvariant
property. Beforeproceedingwenotethatwiththespatialinvariancein

R
s
theLIRdoes
still depend on , which is desirable as it gives us a single parameter with which to
control the trade-off between noise and resolution. We note also that while this method
43
should produce invariance to count rates, there will still be some variation in spatial
resolution. The reason for this is that the intrinsic resolution of PET scanners, and
therefore the matrixP
T
P, is slowly spatially variant, with decreasing radial resolution
inthetransaxialplaneasonemovesawayfromthecentralaxisofthescanner. Sincewe
donotcompensateforthisintheaboveanalysis,wewillseeasimilarradialdependence
inresolutionaswemoveawayfromthecenterofthefieldofview. Itispossibletoalso
compensateforthiseffect,asdescribedin[9],butthiscausesnoisevariancetoincrease
away from the center of the field of view, and in practice it is preferable to allow this
smalldecreaseinresolutionratherthanamplifyingnoise.
Toconclude,theappropriateweightstouseinthepenaltyfunctiontoassureapprox-
imatecountindependentresolutioninthereconstructedimagesare

jl
=
q
P
n
d
i=1
p
ij
2
R
T
T
0
~
B
2
l
(t)
i(t)
dt
q
P
n
d
i=1
p
ij
2
R
T
T
0
~
B
2
l
(t)dt
(2.32)
Onmodifyingthebin-modelikelihoodfunctionin(2.17)toincludethetwoseparate
acquisition frames: [T
1
;T
2
] and [T
3
;T
4
], it is straightforward to show that
jl;j
0
l
has the
followingform:

jl
=
s
n
d
P
i=1
p
2
ij
R
t
2
t
1
~
B
2
l
(t)dt

i
(t)
+
R
t
4
t
3
~
B
2
l
(t)dt

i
(t)
s
n
d
P
i=1
p
2
ij
R
t
2
t
1
~
B
2
l
(t)dt+
R
t
4
t
3
~
B
2
l
(t)dt
(2.33)
Similarly the correct scaling of B
2
(t) to match resolution between the slope and
interceptimagesfortwoframedatais:
=
v
u
u
t
R
t
2
t
1
~
B
2
1
(t)dt+
R
t
4
t
3
~
B
2
1
(t)dt
R
t
2
t
1
~
B
2
2
(t)dt+
R
t
4
t
3
~
B
2
2
(t)dt
(2.34)
44
WehaveshownthatbyimposingadiagonalweightingmatrixD

asshownin(2.33
) to the standard penalty, theLIR will be spatially invariant. The consequent formulas
ofLIR andpenaltyisgivenin(2.29)and(2.28).
2.3 Results
2.3.1 CramerRaoAnalysis
Toexplorethepotentialperformanceofdualtime-pointPatlakimagingweusetheFIM
expressions derived above to compute CR bounds on the estimated Patlak parameters.
The CR bound provides a lower bound on the variance of any unbiased estimator of
a parameter. In practice these bounds are frequently tight and reflect the achievable
variance of maximum likelihood estimators. In the following we first compare CR
bounds with estimator variances computed using Monte Carlo simulation to show that
the bounds are tight. We then use these results as the basis for a comparison of lesion
detection performance based on dual time-point imaging using our Patlak approach in
comparison to fractional SUV. Since the difference between fractional SUV and Patlak
analysis lies largely in the handling of the dynamic rather than the spatial information,
we can obtain insight into the problem through consideration of a simplified problem
corresponding to a single image voxel from which we directly observe Poisson counts
generatedaccordingtothedynamicTACassociatedwithaparticularmodel.
We used a two compartment model with irreversible trapping in the second com-
partment, k
4
= 0. The studies reported in the following subsections used the range of
kinetic parameters listed in Table 2.1 for tumor and background, which were found in
theliteratureonPatlakmodeling.
WeusedtheFIMexpressionstocomputeCramer-Raoboundsforthekineticparam-
etersandcomparedthesetomaximumlikelihood(ML)parametersestimatedfromdata
45
Tumorregion normaltissue
SetID K1 k2 k3 K1 k2 k3
1 0.198 0.228 0.035 0.533 0.980 0.012
2 0.076 0.134 0.115 0.194 0.557 0.054
3 0.088 0.091 0.037 0.117 0.176 0.008
4 0.062 0.055 0.074 0.719 0.788 0.003
5 0.057 0.045 0.077 0.889 1.062 0.004
6 0.069 0.029 0.121 0.845 1.019 0.004
7 0.075 0.078 0.080 0.708 0.865 0.004
8 0.076 0.042 0.082 0.786 0.750 0.004
9 0.064 0.061 0.078 1.151 1.348 0.007
10 0.110 0.084 0.062 0.452 0.589 0.006
11 0.149 0.041 0.098 1.251 1.325 0.009
Table 2.1: Pairs of kinetic model parameters for background and tumors identified in
theliterature
simulated for each model. In this single voxel/detector model no smoothing penalty is
required and parameters were estimated by maximizing the likelihood function (2.5).
Thedataweresimulatedbygeneratingarrivaltimesforphotonscorrespondingtoatime
inhomogeneousPoissonprocessdefinedbytheTACforeachparametersetforarealistic
inputfunction. Weincludedscatterproportionaltothetruerateandrandomsproportion
to the square of the true rate, each scaled to contribute at 30% of the mean true counts
whenaveragedovertheentireTAC.Tocomputethebounds,themeanrandomandscat-
terrateswereassumedknown. Wethenretainedthosearrivaltimescorrespondingtothe
twoframes: 40-45minsand80-85mins. Differentcountlevelswereachievedbyscaling
the input function. We performed 100,000 Monte Carlo trials for each parameter set
andcomputedthemeanandstandarddeviation(s.d.) forthesetrialsforthePatlakslope
parameter.
Table 2.2 compares the CR bound and standard deviation of Patlak slope for the
tumor kinetic parameters in the first row of Table 2.1. The true value of Patlak slope
for these parameters is .0262/min. For all count levels the error in the mean from the
ML estimator is less than 0.5%. Therefore we can regard the estimator as effectively
46
Countlevel 1 2 3 4 5
Meantumorcounts 59.8 69.8 79.7 89.7 99.7
MonteCarlomean 0.0262 0.0263 0.0262 0.0262 0.0263
CramerRaobound(s.d.) 0.0143 0.0133 0.0124 0.0117 0.0111
MonteCarlo(s.d.) 0.0145 0.0134 0.0125 0.0119 0.0112
Table2.2: Comparisonofthes.d. ofPatlakslopeobtainedfromtheCRboundwithmean
and s.d. computed by Monte Carlo simulation with maximum likelihood parameter
estimation
unbiased. This is consistent with the asymptotically unbiased property of the maxi-
mumlikelihoodestimator. TheCRvarianceboundsarewithin2%ofthecorresponding
MonteCarlovalues. ConsequentlyitisreasonabletousetheCRboundasaclosedform
estimateofthestandarddeviationofMLestimatedparametersaswedointhefollowing
subsection. In the following we also use the fact that the ML estimators approximately
follow a Gaussian distribution. This property is illustrated in Fig. 2.1 which shows the
sampledistributionofthePatlakslopeparametersforoneparametersetfor1,000Monte
CarlotrialsoverlaidwithaGaussianwithmatchedmeanandvariance.
Figure2.1: HistogramoftheestimatedPatlakslopefittingaGaussiancurve
47
We also used Cramer-Rao analysis to compare Patlak estimation from two frame
list-mode vs. two frames of binned sinogram data. Again, we used the single voxel
modelwithCramer-RaovaluescomputedusingtheFIMbasedonthebin-modeapprox-
imation(2.15)forlist-modedata,andthetwo-frameextensionof(2.12)forthesinogram
data. Figure2.2showsstandarddeviationsforthetwoestimatorsasafunctionofframe
duration. The starting time for the 2 frames was fixed at 40 and 80min using parame-
ters from the first row in Table 2.1. As one might expect, the methods are matched for
shorter frame duration, but the list-mode method improves relative to sinogram data as
durationincreases.
Figure 2.2: Relative standard deviation of Patlak slope and intercept for list-mode and
sinogrambasedestimationmethodwithdifferentdurationofacquisitionframes.
2.3.2 ROCAnalysis: DualTime-PointPatlakvs. FractionalSUV
Using the results of the preceding analysis we next performed a Receiver Operating
Characteristic (ROC) study to determine the relative performance of dual time-point
48
PatlakandfractionalSUVindifferentiatingtumorfrombackground. Weagainusedthe
singlevoxel/detectormodeldescribedabove.
To compute ROC curves we find the distribution of the Patlak slope parameter esti-
mates for tumor and background. The ROC curve is then formed by plotting the rate
of true positive vs. false positive tumor detection as a function of the threshold applied
to the Patlak parameter. For comparison we also computed ROC curves for fractional
SUV:
%DSUV = (C
2
C
1
)=C
1
=C
2
=C
1
1 (2.35)
where C
1
and C
2
are the ML estimates of activity in each of the two frames (note
thatthenormalizationfactorsusedtocomputeSUVvaluesfromcalibratedPETimages
cancel in the fractional SUV model so arent considered here). With the dual time-point
imaging protocol, many researchers have used %DSUV as a common approach while
the absolute difference of SUV values is not generally used [92], [88], [22] since these
valuesaresensitivetonormalization. Similarly,inourexperimentswerestrictouratten-
tion to %DSUV in comparison with the Patlak method. For the single voxel/detector
model the ML estimator for a static frame is equal to the number of counts observed
over the frame. To determine the distribution of %DSUV we used the Monte Carlo
simulation method outlined above with 10,000 repetitions and computed the sample
distribution.
Examplesofthedistributionswithkineticparameterset1fromTable2.1areshown
inFig. 2.3withROCareas0.7565and0.5166forPatlakand%DSUVrespectively.
We repeated the same procedure for a range of different count rates, frame dura-
tions and acquisition start times, in each case computing the area under the ROC curve
for the Patlak slope and %DSUV measures. In each case, randoms and scatters were
added to the data. Area under the ROC curves was computed using the CR analysis
49
Figure 2.3: An example of the distribution of Patlak slope parameter (upper) and
%DSUV(lower)fortumorandbackground(BG)
described above for both list-mode and sinogram based Patlak estimation. These are
shown relative to %DSUV analysis in Fig. 2.4. Figure 2.4 upper left shows ROC areas
for 5 different count levels for parameters equal to the first row of Table 2.1 with two
5min frames starting at 40min and 80min respectively. Figure 2.4 upper right shows
ROCareasforthelowestcountlevelinTable2.2withacquisitionagainbeginningat40
and 80mins. In this case we kept the total acquisition to 10min but varied the relative
duration of the first and second frames. Figure 2.4 lower left keeps the same count rate
and total duration and 40min inter-frame interval, but varies the starting point for the
first frame. Figure 2.4 lower right shows the ROC areas for Patlak and %DSUV for all
kinetic parameter sets except the first row in Table 2.1 for two 10min frames beginning
at 40 and 80mins and a count rate approximately equal to the lowest value reported in
Table2.2.
From Fig 2.4 we see that list-mode based estimation slightly outperforms the sino-
gram based method in terms of the AUC, although the differences are small. These
50
Figure2.4: DetectionoftumorsrelativetobackgroundusingrateparametersfromTable
2.1. We compare area under ROC curves (AUC) for Patlak slope (sinogram and list-
modebasedestimation)and%DSUVasafunctionof(a)countrate,(b)durationoffirst
frame, and (c) first frame start time. Panel (d) shows the area under the ROC curves
averagedoversimulationsfordifferentparametersetslistedinTable2.1
resultsindicatethereislittledifferenceinperformancebetweenlist-modeandsinogram
based estimation for the range of parameters tested here. However the fact that the
list-mode approach always performs as well or better than sinogram based estumation
supportsitsuseinsubsequentstudies.
We also observe consistently better performance in distinguishing tumor and back-
groundwithPatlakcomparedto%DSUVeventhoughbothmethodsusethesamedata.
Since %DSUV has already shown advantages over SUV in cancer staging as we have
51
reviewed in the introduction, these results are promising indications of practical util-
ity for our dual time-point strategy in whole body cancer imaging. These simulations
also may provide guidance for optimizing scan parameters for two frame imaging. By
extending these simulations to a wider range of parameters we can optimize start times
for the two frames and their relative duration give a fixed total scan time so as to maxi-
mizeareaundertheROCcurve.
2.3.3 SimulatedImagingStudies
Forthesestudieswesimulatedasmallscale3DPETsystem(diameter: 148.4mm,detec-
tor size: 2.423mm2.423mm; number of rings: 4) with a total of 13 2D sinograms each
with 84 angles of view by 96 radial LORs. A uniform cylinder phantom of diameter
31.4mm is centered in the scanner and contains 5 cylinders (tumors) of diameter 1.0,
1.8, 2.6, 3.4, and 4.2mm as shown in Fig. 2.5(a). We simulated time activity curves
for tumors and background using the first set of kinetic parameters in Table 2.1. The
blood input function, tumor and background TACs are shown in Fig. 2.6. As in the
Cramer Rao analysis, we added scatters and randoms to the data. Both processes were
assumed to be separable, with the temporal variation in scatter and randoms propor-
tionalrespectivelytotheinstantaneoustrueandsquared-trueactivityintegratedoverall
LORs. Total activity contributed by the scatters and randoms were both set at 30% of
thetotaltrueactivity. Themeannumberoftotalcountsforthefollowingstudiesoverthe
entire scanning time was 17M. The TACs for each sinogram element are computed by
forward projection through a system matrix generated for the scanner using the model
described in [9]. Each of these TACs represents the time-inhomogeneous Poisson pro-
cesses that determine emission rates for each LOR. List-mode data are then generated
by sampling from these TACs. In the results presented here, we window the data to
52
retain only those events corresponding to two five minute frames starting at 40min and
80min. Wegeneratedatotalof100trialsforourMonteCarlosimulations.
Figure 2.5: (a) Uniform cylinder (background) with five cylinders (tumor) of different
sizes. (b) The mean (over 100 Monte Carlo trials) of the Patlak slope image recon-
structedfromtwo5minframelist-modedatastartingat40and80minwithTACsasin
Fig. 2.6.
Weusedthemodified4Dincrementalgradient(4DIG)algorithm[63]toreconstruct
the Patlak slope and intercept images. The mean of the Patlak slope images over the
100 trials is shown in Fig. 2.5b. Wealso computed static penalized ML reconstructions
of each frame using the method described in [9] with spatially variant smoothing to
ensure count independent resolution as described in [34]. To ensure a fair comparison
wematchedresolutionofthestaticandPatlakimagesasdescribedbelow.
For each of the 100 trials we computed the average Patlak slope parameter for each
of 5 tumor regions of interest (ROIs) and 5 corresponding background ROIs with sizes
matching those of the tumor ROIs. Similarly we computed the mean intensity in each
ROI for each static reconstruction and then computed %DSUV values. The mean and
s.d. resultsforPatlakand%DSUVforeachofthefiveROIsareshowninFig. 2.7.
53
Figure 2.6: Simulated time-activity-curves for tumor, background, and the blood input
functionwithkineticparametersfromrow1inTable2.1.
Asacomparisonbetweenlist-modeandsinogrambasedPatlakestimation,wecom-
pute the standard deviation of the mean Patlak slope within the biggest tumor. The
results are displayed in Fig. 2.8. The relative sd of the list-mode based estimation is
uniformly smaller than that of the sinogram based method, which is consistent with
whatweobserveintheCramerRaoanalysis.
AssumingGaussiandistributionsfortheROImeansineachcase,wecomputedareas
under the ROC curve for detection of each of the tumors relative to a size-matched
background ROI. These results are shown in Fig. 2.9. In this figure we include ROC
resultsforMonteCarlosimulationsforparameterscorrespondingtothefirstthreerows
inTable2.1andincluderesultsfor10minaswellas5minframes.
WeseethatasthesizeoftheROIincreases,s.d. in%DSUVdropsinFig. 2.7leading
to improving ROC values in Fig. 2.9 with increasing ROI size. In the Patlak case, the
54
Figure 2.7: Mean and s.d. (error bars) for each of five ROIs shown in Fig. 2.5(a) for
list-mode Patlak and %DSUV. %DSUV values were computed by applying (2.39) to
activityineachframeafteraveragingovertheROI.
performance also improves as ROI size increases, but overall the s.d.s are far lower
than in the %DSUV images so that ROC results are considerably better for the smaller
ROIs. Theseresultsareconsistentwithourearliersinglevoxel/detectorstudyandagain
indicate potential benefits of using Patlak over %DSUV for detection, especially for
smallertumors.
2.3.4 Count-InvariantResolution
To investigate the effectiveness of our technique to ensure count-invariant resolution
wesimulatedacoolcylindricalbackgroundwithonehotandonewarmsmallercylinder
insert,eachwithapointsourceattheircenter. Theentirephantomhasthesamedynamic
activitycurveateachvoxelbutwiththerelativescalepointsource: hotcylinder: warm
cylinder: backgroundof50: 10: 1: 0.4. Thephantomandaprofilethroughthe2point
sourcesareshowninFig. 2.10.
55
Figure2.8: Relatives.d. ofPatlakslopeestimatecomputedfrom100Monte-Carlotrials
forthebiggesttumor.
We simulated data for this phantom using the same method as described earlier.
Randoms and scatters were added to the true events as described above. The data for
the entire acquisition time was 23M counts. We then reconstructed images of Patlak
slopeandinterceptusingthecount-invariantspatialweightingofthesmoothingpenalty
defined in (2.37) and, for comparison, a spatial weighting equal to the inverse of the
Euclidean distance between each pair of pixels. The reconstructed Patlak slope image
and a profile for one noisy reconstruction using these two different penalties are shown
in Fig. 2.11. The figure shows that the reconstructed height of the two point sources
is quite different using the spatially invariant weighting but approximately equal when
usingthecount-dependentspatiallyvariantweighting.
To explore this performance quantitatively we computed the peak value of the local
impulse response at the location of the point source in the hot and warm cylinders as
56
Figure2.9: Areaunder theROCcurveforkinetic parametersets1, 2, 3inTable2.1for
Patlakand%DSUVbasedtumordetection. Seetextfordetails.
a measure of contrast recovery coefficient (CRC). These values are tabulated for Pat-
lak slope and intercept in Table 2.3 for three different values of the global smoothing
parameter when using a uniform spatial weight for the penalty term. These results
clearly show that the hot and warm cylinders have differing CRCs for each value of.
Furthermore, the slope and intercept also have distinct CRCs. In contrast, for the same
data using a shift-variant weighting we obtain CRCs that are approximately equal in
the hot and warm cylinders for both the slope and intercept images for each value of
as shown in Table 2.4. Furthermore, we see that after modifying the count rate for
57
Figure 2.10: (left) Digital phantom for validating count-independent resolution and (b)
profilethroughthe2pointsourcesinthephantom
Beta 0.1 0.1 0.5 0.5 1 1
Warm Hot Warm Hot Warm Hot
slope 0.430 0.418 0.311 0.258 0.236 0.191
inter 0.460 0.532 0.450 0.465 0.429 0.433
Table2.3: CRCat2pointsourcesforthelowcountlevelandtheuniformspatialpenalty
Beta 0.1 0.1 0.5 0.5 1 1
Warm Hot Warm Hot Warm Hot
LowCountslope 0.471 0.469 0.428 0.434 0.350 0.352
LowCountinter 0.467 0.472 0.437 0.437 0.363 0.367
HighCountslope 0.478 0.512 0.448 0.449 0.374 0.372
HighCountinter 0.475 0.514 0.470 0.483 0.421 0.440
Table2.4: CRCat2pointsourcesforlowandhighcountlevelsandourmodifiedspatial
penalty
this study from low to high, we again obtain CRCs that are invariant for fixed and
approximatelyequalforboththehighandlowcountrates.
58
Figure 2.11: Profiles through one reconstructed Patlak slope image with (left) spatially
invariantsmoothingand(right)spatiallyvariantsmoothingforcount-independentreso-
lutionusing(2.37).
2.3.5 MatchingResolutioninPatlakand%DSUVimages
TomakethecomparisonsmeaningfulitisimportantthattheresolutionofthePatlakand
staticimagesfromwhichwecomputed%DSUVvaluesarematched. ThepenalizedML
staticimageswerereconstructedusingthespatiallyvariantsmoothingmethoddescribed
in [34]. This ensures that resolution is matched in the two frames in each dual frame
studyandfurtherthatitcanbecontrolledusingasingleparameterwhichplaysthesame
role as the parameter in the Patlak image smoothing function in (2.6) and (2.21). We
computed a calibration curve of vs. full width at half maximum (FWHM) resolution
59
as shown in Fig. 2.12. The curve is based on the approximation that the mean of the
reconstructedimageequalsthereconstructionofthemeanofthedata[90]. Wetherefore
reconstruct noiseless data for a point source at the center of the field of view and then
computetheFWHMresolutionformultipledifferentvaluesofthesmoothingparameter
toobtainthecurveinFig. 2.12.
Figure 2.12: Resolution versus log of smoothing parameter in static MAP reconstruc-
tionscomputedatthecenterofhottumorforthephantominFig. 2.5.
The resolution of the Patlak parametric image vs cannot be obtained from the
mean of the data because we are using a list-mode likelihood function that depends on
the individual arrival times. Rather than modify the algorithm to compute resolution
based on bin-mode likelihood we used a Monte Carlo method to empirically determine
resolution. We computed the FWHM resolution for the mean image averaged over the
100MonteCarlotrialsdescribedaboveforthepointsourceatthecenterofthehottumor
forthephantominFig. 2.5. Theresolutionwasmeasuredtobe2.05mm. Wethenfound
60
the value that matched this resolution in the static images and use this for all of the
studiesreportedintheresultsection. Amorecompletestudywouldalsoinvestigatethe
relativeperformanceof%DSUVandPatlakmethodsasafunctionresolutionparameter
,howevertheprimarydifferencebetweenthetwomethodsisthemannerinwhichthe
dynamic data is handled. Both use the same spatial model so we would expect relative
performancetobesomewhatrobusttochangesinresolution. Limitedsimulationstudies
wehaveperformedbutnotincludedhereappeartosupportthisconjecture.
2.4 Discussionandconclusion
TheROCresultsshowthatforasingledata/singlevoxelCramer-Raostudyoverawide
rangeofparametervalues,Patlakslopeestimatesfromdualtime-pointdataconsistently
out-perform fractional SUV measures computed from the same data. The advantage of
theCramer-Raoapproachisthatwecanrapidlyevaluateperformanceformanydifferent
acquisition parameters. Consequently in addition to comparisons to fraction SUV, we
canalsousethisapproachtooptimizeacquisitionintermsofstartandendtimesforthe
twoframeswhendesigningwholebodyprotocols.
Ouranalysisofresolutionprovidesthetheoreticalbasisformodifyingthecostfunc-
tion to ensure that reconstructed Patlak images have a resolution that does not depend
on activity or count rate and furthermore that the resolution of the slope and intercept
parameters are matched. This is important for the purposes of qualitative and quan-
titative interpretation of the Patlak values. For example, matched resolutions allow
meaningfull comparisons of Patlak slope values, which are inevitably affected by par-
tial volume effects for small lesions, across subjects and over time for a single subject.
Removingdependenceonactivitywillmakethesevaluesmorereliablebiomarkersthat
canpotentiallybeusedforstagingandassessingresponsetotherapy.
61
Thesimulatedimagingstudiesmakeuseofthecount-independentresolutionmethod
to compare performance of Patlak and fractional SUV methods at matched resolutions.
Through Monte Carlo simulation over a range of lesion sizes, we see in Figs. 2.7 and
2.9 that, as with the Cramer Rao analysis, Patlak slope parameters are more accurate in
differentiating tumor and background than is fraction SUV. This is particularly the case
forthesmallerlesionswherethemuchsmallerrelativevarianceofPatlakmeasureslead
tobetterdifferentiationoftumorandbackground.
In summary, the results presented support the conjecture that Patlak values esti-
matedfromdualtime-pointPETdatacanproducemoremeaningfulindicationsoftracer
dynamics than are extracted using the fractional SUV approach. Comparisons of list-
modeandsinogrambasedPatlakimageestimationshowlittledifferenceforshort(5min)
duration frames but with improvements in list-mode performance relative to sinograms
as the duration increases. In practice there may be little difference in performance
between the two approaches, but the list-mode approach always performs at least as
wellasthesinogrammethod.
Our reason for pursing the two time-point method, rather than the more standard
method of acquiring a full dynamic scan, is that we plan to extend this approach to
whole body studies where the patient must be stepped through the scanner to produce
the whole body image. The method described in this paper is directly applicable to two
time-pointwholebodydata.
Inprinciple,dataacquisitionforwholebodyPatlakimagingrequiresasimplemodi-
ficationofcurrentwholebodyprotocolsinthatthepatientissimplysteppedthroughthe
scannertwiceratherthanonce. However,anumberofpracticalissuesariseinacquiring
and processing this data. First there needs to be a delay of 40-60mins between frames.
Assuming 5mins per bed position and 5 positions to cover the patient, the patient must
remain in the scanner for 65-85mins. It is difficult to remain stationary for that long, so
62
thatcoregistrationbetweenframeswillbeneededwhenprocessingpracticallyacquired
data. Asecondcomplicationisthatabloodinputfunctionisrequired. Arterialsampling
isnotpracticalfortheenvisionedapplication,butapopulationormodelbasedapproach
can be applied for this purpose. Finally, the usual calibrations for attenuation, scatter,
randomsandsystemsensitivityneedtobeincludedwhenimplementingthismethodfor
experimental whole body data. Methods for addressing these practical issues, and the
resultingwholebodyreconstructions,willbedescribedinthenextchapter.
63
Chapter3
Clinicalapplicationofdirect
wholebodyPatlakestimation
3.1 Abstract
We describe a wholebody Patlak parametric image estimation method for
18
F-FDG
positron emission tomography (PET) scans. A practical method for processing clinical
PET data has been developed, including estimation of blood input function, modeling
of random and scattered events, and motion correction. We present preliminary clinical
results comparing Patlak images with the standardized uptake value (SUV) measure-
ments. Methods: We have developed a dual time-point PET scan protocol. Whole-
body Patlak parametric image are then estimated from the acquired list mode data. An
input function is estimated using a hybrid method combining population information
and imaging data. Dynamic models are developed for random and scattered events. A
non-rigidregistrationmethodisusedtocompensateforthemotionbetweentwoframes.
We selected regions of interest (ROIs) in the liver, thigh muscle and lumbar spine (L2)
for a set of eleven patients and computed the mean and standard deviation (s.d.) of the
Patlakslopeandinterceptoverthepopulation. WealsocomputedSUVvaluesforthese
ROIs for each frame, from which the percentage change in SUV (%DSUV) was also
calculated. Results: Comparisons of ROI statistics show that Patlak slope and inter-
cept have consistently lower fractional standard deviation across the population than
64
SUV values in the same ROIs. In comparison to %DSUV the fractional standard devi-
ations are considerably lower (by a factor of three or more). Population-wise s.d. of
the ROI vlaues in normal livers, thigh muscles, and lumbar spines for Patlak slope are
15.4%,12.2%,and15.1%,comparedto17.1%,22.3%,and20.3%fortheSUVofearly
frame and 15.6%, 21.0%, and 24.4% for the SUV of late frame. Corresponding s.d. of
%DSUV are 49.9%,100%, and 59.5% respectively. Conclusion: Wholebody Patlak
parametric imaging is feasible for clinical
18
F-FDG PET studies. Patlak images show
consistently lower inter-subject variability than the semiquantitative SUV-based mea-
sures. PatlakestimationhaspotentiallybettertumorcontrastthanSUV-basedmeasures.
3.2 Introduction
Positron Emission Tomography (PET) using
18
F-FDG is widely used in clinical oncol-
ogyfortumordetection,stagingandassessmentofresponsetotherapy. Currentclinical
practiceofwhole-body
18
F-FDGPETinvolvestaticscansthattake2-5minutesperbed
position(15-20cm)startingapproximately50-70minutesafterinjectionassuggestedby
European Organization for Research and Treatment of Cancer (EORTC) and National
Cancer Institute (NCI) guidelines ([93], [1]), when the tracer distribution is assumed to
be stable. Standardized uptake value (SUV), a semi-quantitative measurement, is com-
putedandusuallyinterpretedvisuallybyradiologists([94]),althoughrecentlyquantita-
tivecriteriaforevaluatingresponsetotherapyinsolidtumorhasbeenproposed([13]).
OnelimitationtotheapplicationofSUVisthehighvariability. Manytechnological
andbiologicalfactorsaffecttheaccuracyofSUV([95],[96]). Themostimportantfactor
is probably the uptake time (time between tracer injection and PET scan). Standardiz-
ing the uptake time improves the reproducibility of SUV ([97]) but does not solve the
65
problemcompletelybecausethepeaktimeofFDGindifferenttumorsmayvarysignif-
icantlydependingonthetumortissuekinetics([1]). AnotherlimitationisSUVmaynot
beabletodifferentiatemalignanttumorsfromphysiologicalorbenignprocessessuchas
inflammation([98],[99],[100]),becauseSUVdoesnotprovideareliablemeasurement
ofthekineticsofFDGinthetissue([95]).
To measure the underlying dynamics of
18
F-FDG uptake, a dual-time-point SUV
method has been proposed ([101]). This method measures the percentage change in
SUV (%DSUV), which is the change of SUV values between two time points, for-
mulated as %DSUV=(SUV2-SUV1)/SUV1. Alkhawaldeh et al ([20]) have shown that
dual time-point PET improves diagnostic accuracy relative to standard SUV, increasing
sensitivity and specificity for malignant lung nodules diagnostic, especially for small
lung lesions that have low SUVs. Similar results have been shown for high-grade brain
tumors ([21]) and breast cancer ([102]). However, contradictory results have also been
reported. For example, Cloran et al. reported that dual-time-point
18
F-FDG PET may
not be of benefit in the assessment of pulmonary nodules with maximum SUV of less
than 2.5 on initial imaging ([24]) . Chan et al. suggested that dual-time-point
18
F-FDG
PET does not improve the overall evaluation of pulmonary lesions ([25]). This is not
surprising because %DSUV is based on SUV measurement thus suffers from the same
limitations.
Kinetic analysis has the ability to determine metabolic rate and greater robustness
to variations that affect SUV measurements. It is also possible to distinguish tumor
frombenignprocessesandnormaltissuesfromthekineticinformation,whichhasbeen
demonstratedfromanimal([22])andhuman([103])studies.
A full kinetic analysis requires measurement of the time course of the radioactivity
in tissue and blood from the moment tracer is injected, which is impractical in a clinic
environment. As a result, graphical methods such as Patlak analysis have been used
66
frequentlyfor
18
F-FDGPETduetoitssimplicityandrobustness([44]). Aninvestigation
of patients with colon cancer metastatic to liver found that the rate constant determined
from Patlak analysis was far less variable compared to SUV, although they are strongly
correlated ([31]). Studies have shown that Patlak slope is a more accurate index of
glucosemetabolicratethanSUVbecauseoftheun-metabolizedFDGmeasuredbySUV
and the inability of SUV to account for the available dose ([104]). It has also been
suggested that the tumor volume can be better delineated from Patlak slope image than
SUV image ([105]). In addition, Studies have shown the Patlak method has improved
specificityatasimilarsensitivityforbreastcancer([106])andskeletalcancer([103]).
WholebodyPETimagingisusedextensivelyintheclinicforstagingandmonitoring
the treatment response of primary and metastatic tumors. Due to the limited axial field
of view (FOV), the patient is usually stepped through the scanner gantry such that the
body between skull and mid thigh is covered. The clinical setting poses challenges
to dynamic imaging, which usually requires imaging at different time points. This is
apparently in conflict with multi-bed position wholebody imaging, which only allows
partial frame data for each bed position. We solved this by developing an estimation
method that uses partial data for each bed position, so as to make dynamic wholebody
imagingpossible.
A method to estimate Patlak parametric images from list-mode PET data has been
proposedinourearlierpublication([42]). UsuallyinPatlakanalysis,multipleframesof
PET data are acquired and reconstructed, then ROIs are drawn and a linear regression
modelisappliedtoestimatethePatlakslopeandintercept. Directestimationfromsino-
gramhasalsobeenproposed([107]). Ourmethodisdifferentfromtraditionalapproach
or sinogram-based approach by directly using the list mode PET to estimate the Pat-
lak parameters, which fully utilizes the temporal information of the events. We note
67
that recently Karakatsanis ([35], [36]) proposed a clinical protocol for wholebody Pat-
lak estimation composed of an initial 6 min dynamic scan over the heart followed by
a sequence of PET scans (6 pass x 7 bed position each). A hybrid linear regression
approach was used to obtain the Patlak image. In contrast, we propose to use the list
modedataandworkinaBayesianframework.
In this paper, we describe the practical considerations of Patlak parametric imaging
for clinical whole body
18
F-FDG PET. In the following chapters we will first illustrate
thealgorithmofourmethod. Wethendescribehowwemodelthescattersandrandoms,
correct patient motion, and estimate the blood input function. Finally we show patient
resultsandcompareROImeasurementsofPatlakparameterswithSUVand%DSUV.
3.3 Materialsandmethods
3.3.1 WholebodyPatlakAnalysis
The Patlak model provides an approximate solution to the two-tissue compartment
model for irreversible tracers in the steady state period. The rate function in sinogram
canbewrittenas([42]):
(t) =e
t=
(P
B(t))W +r(t)+s(t) (3.1)
where in this vector form representation W = [w
T
1
w
T
2
:::w
T
N
]
T
is the Patlak image
and each element w
n
is constructed by the Patlak slope K
i
and the Patlak intercept q.
B(t) = [B1(t)B2(t)] where B1 is the integral of blood input function and B2 is the
blood input function.
is the Kronecker product. The exponential term accounts for
radioactive decay of the tracer (e.g., half-life of
18
F is 110 min); P is the PET system
68
matrix ([9]). The additional terms r(t) and s(t) are the randoms and scattered coinci-
dence rate functions. A dual-pass whole-body listmode scan protocol for performing
Patlak analysis is proposed. For each bed position with acquisition time [t
1
;t
2
] and
[t
3
;t
4
], the penalized log-likelihood function of event arrival times is derived with an
inhomogeneousPoissonmodelforthearrivingphotons:
F(W) =
X
i
x
i
X
k
log
i
(a
ik
)

X
i
(
Z
t
2
t
1

i
(t)dt+
Z
t
4
t
3

i
(t)dt)
W
T
R
S
W
(3.2)
where a
ik
denotes the arrival time of the kth photon at LOR (line-of-response) i in the
intervals [t
1
;t
2
] and [t
3
;t
4
],x
i
is the number of events detected in LORi. R
S
computes
the weighted sum of squared pairwise difference between each voxel and its 26 nearest
neighborsforboththeslopeandinterceptparameters. Maximizingthispenalizedlikeli-
hood function provides a direct estimation of the Patlak parametric images. We apply a
spatially variant smoothing parameter to achieve count-independent and approximately
spatially invariant resolution of the Patlak image ([42]), based on our previous work on
static PET image reconstruction ([34]). The weights in R
s
are functions of the Fisher
information matrix chosen to compensate for the count-dependent effect. The resolu-
tionoftheslopeandinterceptimagesarealsomatchedbyrescalingthetemporalmatrix
B. Unlike OSEM reconstructions, the resolution of Patlak images does not vary with
activitylevel.
69
3.3.2 DataProcessing
BloodInputFunctionEstimation: Routine measurement of the input function using
arterialsamplingisimpracticalinclinicalprocedure. Ahybridapproachthatcombinesa
populationbasedestimateforearlyscantime(t1h)withanimage-basedexponential
modelfortherestofthescantime(t1h)isdevelopedforthiswork.
Thehybridinputfunctionestimationwasevaluatedusingdynamic
18
F-FDGdatasets
of a previous liver cancer patient study. Seven patients were scanned for 1 hour imme-
diately after injection, using a listmode scan protocol which only covers 1 bed position
where the liver was located. ROIs in the abdominal aorta were identified in each sub-
ject from coregistered CT images. The list mode data was histogrammed into multiple
frames and each frame reconstructed separately using an iterative method ([9]). The
activity in the aorta of each reconstructed frame was used to fit a continuous model.
Each input function was then scaled to have unit integral under the curve. The stan-
dard input function was then generated by averaging over these scaled input functions.
The average RMSE (relative mean square error) of these scaled input functions to the
standard input function is 0.03), strongly supporting the validity of population based
approach.
The multi-exponential function is commonly used to model the input function.
Ferl ([108]) modeled the blood input function as a weighted sum of 4 exponentials
B(t) = A
1
exp(b
1
t) + A
2
exp(b
2
t) + A
3
exp(b
3
t) (A
1
+ A
2
+ A
3
)exp(b
4
t).
Wong et al ([109]) also modeled the input function as a weighted sum of 4 exponen-
tials,withfollow-upPatlakanalysisshowingthatthePatlakslopeK
i
estimatedwiththe
approximationmethodandthebloodsamplebasedmethodmatchwell. Usuallybecause
of the different decaying rates for the exponential functions, only one function will be
70
dominant after sufficient time elapses (Table 3 in ([108])). For our purpose the expo-
nential model is used 1 hour after injection, and 2 additional blood samples computed
fromtheaortaintensityafter1hourwillbeusedforthemodelfitting.
The bed position containing the abdominal aorta in the dual-pass was chosen for
this purpose. Usually this bed position starts 65-min in the first pass and 105-min in
thesecondpassafterinjection. StaticMAPimagereconstructionisperformedforthe2
framesofdata. AnROIwithintheaortaismanuallydelineatedandtheaverageactivity
in the ROI computed for both frames, followed by decay correction, to obtain the 2
bloodsamples.
Thefollowinghybridinputfunctionestimationmodelisthenusedforeachincoming
patient,wheret
1
isapproximately1hourafterinjection. Theconstants;;
aredeter-
minedbyfittingthemodeltothe2sampleswiththeconstraintthattheinputfunctionis
continuousatt
1
.
C(t),f
C
standard
(t) tt
1
exp(
t) tt
1
(3.3)
Thisapproachcouldbemodifiedinastraightforwardmannertoincludedatafromaddi-
tionalframesthatalsoincludelargearterialvessels.
ModelingofRandomsandscatters:Acontinuous-timeestimateofbothrandoms
and scatters is needed for direct dynamic image reconstruction. The whole-body proto-
col is composed of 2 listmode acquisition frames each with 5 bed positions, in prepro-
cessing10pairsofstaticrandomsandscattersestimateswereobtainedfrom10listmode
files.
Thetemporalandspatialdistributionofscattersandrandomsareassumedtobesep-
arable([40])foreachofthe2bedpositions. Sothattheestimatedrandomsratefunction
can be modeled as the product of the integrated random sinogram and a temporal vari-
ation function reflecting the dynamic activity change of randoms over time. Similarly,
71
thescatterratefunctioncanbemodeledastheproductoftheintegratedscattersinogram
andatemporalvariationfunctionreflectingthedynamicactivitychangeofscattersover
time.
TheintegratedrandomsarecomputedfromtheproductofsingleratesforeachLOR
in each bed position. The integrated scatter is computed from the sinograms generated
by histogramming the list-mode data, in combination with a co-registered x-ray CT
scan, using sofware supplied by the manufacturer (Siemens PET image reconstruction
e7 tools). The temporal rate function of scatters for each bed position is assumed to
be proportional to the rate function of prompts, and then normalized to unity over the
frame. The rate function of all prompts for each bed is obtained by first dividing the 5
min time into ten 30 sec subsets and counting the number of prompts for each subset,
andthenperforminginterpolationbetweenthe10datapoints. Asimilarprocessisused
forestimatingtherandomsratefunction,butallowsittovaryproportionaltothesquare
ofactivity.
PatientMotionCompensation: The patient motion compensation is integrated
into the Patlak image reconstruction. We assume that the intra-frame motion is neg-
ligible and only correct for the inter-frame motion between the two passes. The static
PETimagesforeachframeateachbedpositionwerefirstreconstructedwithOSEMfor
9 patients with no detected tumor in report, and with MAP (maximum a posteriori) for
the other 2 patients with detected tumors, using a spatially variant smoothing function
selected to achieve a constant count independent resolution ([34]). The deformation
fields obtained with OSEM and MAP images are quite consistent, given that qualita-
tively the OSEM and MAP images are similar to each other. Therefore, either can be
used for normal size ROI analysis. The reason we use MAP images for the 2 patients
withtumoristhattheresolutioninMAPiscountindependentandmoreuniform,which
mayleadstobetterregistrationforthesmallstructureliketumor.
72
We investigated two approaches to perform registration. For 9 studies the nonrigid
registration Matlab toolbox ([110]) was used, and the registration was performed in an
iterative and multi-resolution way by deforming the moving image with B-splines and
align the images in each grid size. For the 2 patient studies with tumor the registra-
tion function of SVREG was used ([111], [112]), based on the extension of a surface
registrationmethodtovolumetricregistrationintheapplicationofbrainimageregistra-
tion. Both methods are set to obtain the deformation fields between 2 listmode frames
by maximizing mutual information of the two static images reconstructed with OSEM
or MAP. The registration quality was improved by subtracting the 2 images after reg-
istration to visually identify any mismatch and adjust the smoothing parameter of the
penalty and the initial grid size (for MATLAB code) if necessary. The MATLAB reg-
istration code is slower but captures large motions more rapidly by initializing with a
coarsegridsize. TheSVREGmethodisfasterbutpotentiallysuffersfromartifactwhen
the smoothing parameter is set small to allow larger deformation, possibly because the
surface-to-volume approach is more ill-conditioned for small penalty in mutual infor-
mationbasedregistration. However,withcarefultunning,itcanyieldbetterregistration
qualityindetailswhentheoverallmismatchisnotsignificant.
ThePatlakimagesarethencomputedusingtheobtaineddeformationfield. Sincewe
cannot warp the data for the 2
nd
list-mode data, the reconstruction method is modified
so that the spatial mapping between the 1
st
and 2
nd
frame is propagated through the
forwardandbackprojectionoperatorswhencomputingthegradientofthecostfunction
withrespecttothe 2
nd
frame([42]).
3.3.3 Patients
Eleven patients referred to USC PET Center for whole-body
18
F-FDG PET scan were
imaged. The study has been approved by the Institutional Review Board at USC, and a
73
written informed consent form was obtained from each patient. The patients include
6 males and 4 females, with weight 81.119.2 kg. The
18
F-FDG dose injected is
532.835.6MBq.
3.3.4 ImagingProtocol
The patients were scanned using a modified protocol with two listmode frames in addi-
tion to the routine clinical whole-body PET/CT scan. The two listmode frame start 60
minand100minafterinjectionrespectively,eachcomposedwith5bedpositions. Each
bedpositionis5-minlong. Thenon-contrastCT(10patients)andthecontrastenhanced
CT (one liver cancer patient) in the clinical whole-body scan was used for attenuation
correction. The Patlak parameters for the contrast CT patient is quantitatively close to
the statistics for the normal CT patients 3.1. The routine clinical whole-body PET/CT
scanwasusedforstandardpatientcare. Thelist-modefilesoftheadditionaltwoframes
wereusedforPatlakanalysis.
The whole-body patient studies were performed on a Siemens Biograph PET/CT
scanner (Siemens Molecular Imaging, Knoxville, TN). The 3D PET system is 855.20-
mm in diameter with 55 detector rings. The reconstructed image has 109 planes with
4.05-mmslicethicknessand4-mmvoxelsizeinthetransaxialplane. Staticreconstruc-
tions at matched resolutions were reconstructed using manufacturer-supplied recon-
struction software, from each frame with OSEM (3 iterations and 14 subsets) with
PSF(point-spread-function)and5-mmGuassianpostsmoothing. Wealsoreconstructed
MAPwholebodyimagesusingthe methoddecribedin([9])for eachofthetwoframes,
with smoothing penalty chosen to approximately match the profile of MAP and OSEM
reconstructedimagesintheabdomen. ThePatlakslopeK
i
andinterceptq imageswere
reconstructed with algorithm described in ([42]). The images were then stitched to
form the wholebody images. Finally for resolution matching purpose, the wholebody
74
Palakimagesarepost-smoothedwithGaussianfiltersothattheprofilesacrossabdomen
approximatelymatchtheprofilesoftheOSEMimages.
3.3.5 StatisticalAnalysis
The SUV was computed in each listmode frame as SUV =
C(kBq=g)
dose(MBq)=weight(kg)
. The
percentage change in SUV was then computed as %DSUV = (SUV2-SUV1)/SUV1,
where SUV1 and SUV2 are the SUV images from frame 1 and frame 2 respectively.
For%DSUVonlyROIstudieswereperformedbecauseoftheexcessivenoisyduetothe
division operation. The Patlak images were reconstructed with smoothing parameters
chosentomatchtheresolutionsofSUVimages.
ROIs in the liver, thigh muscle, and lumbar spine (L2) were identified for the SUV,
%DSUV, and Patlak images, from the coregistered CT image. The mean and variance
of the Patlak parameters as well as SUV and %DSUV values are computed. Note that
the ROIs were propagated from the first to second wholebody SUV image after coreg-
istraitonusingdeformationsmentionedinthemotioncompensationsubsection.
For one of the patients with breast tumor and inflammation, the contrast of tumor
ROI against background is computed. ROIs in the SUV images are identified to strictly
lie in the tumor and background to minimize partial volume effect. Considering the
potentialerrorinmisregistrationofCTandPET,theROIsareidentifieddirectlyinPET
imageinsteadofCT.ThemeanvaluesofROIsobtainedfromtheSUVimagesandPat-
lakimagesarethencomputedtocalculatethecontrastastumor
mean
=background
mean
.
75
3.4 Results
3.4.1 InputFunctionEstimation
For each of the 7 liver cancer patients in the 1 hour
18
F-FDG dynamic study, data was
sorted into a total of 20 frames (6x30sec, 7x1min, 4x5min, 3x10min). Each frame was
reconstructedwithOSEM(3iterationsand14subsetsand5-mmGuassianpostsmooth-
ing) and the activity in a hand-drawn ROI in the aorta is computed for each frame to
obtain the blood samples. A continuous time function was fitted using spline inter-
polation and the curve was normalized to integrate to unity from 0 to 60 min. The
7 normalized input functions are shown in Figure 3.1A. They are clearly very similar
(average RMSE (relative mean square error) to their mean is 0.03). The standard
input function is thencomputed as the average of the 7 normalized input functions. For
each patient 2 samples are obtained from the aorta ROI in MAP reconstructed images
to fit the exponential function, and the standard input function is scaled to meet the
value of exponential function at 1h so that the whole input function is continuous. Fig-
ure 3.1B shows the estimated blood input function for one patient, which is derived by
concatenating the scaled standard input function (0-60 min) and exponential approxi-
mation (60-140 min). The method takes the values in MAP images and does not need
to represent the blood input function with specific unit, because the Patlak model cali-
brates itself in cancelling out the scale in regressor and regressant. In Figure 3.1B we
extendtheexponentialapproximationtoasearlyas35minafterinjection,whichyields
a25minoverlapregionwiththescaledstandardinputfunction. Itisclearthattheexpo-
nential approximation matches well with the weighted standard input function in this
overlapped region (RMSE 0.02). This indicates that after 35 min post-injection, the
bloodinputfunctioncanbeaccuratelyapproximatedbyasingleexponentialfunction.
76
0 500 1000 1500 2000 2500 3000 3500 4000
0
0.5
1
1.5
2
2.5
x 10
−3 7 normalized blood input function
time (seconds)
normalized unit
A
0 2000 4000 6000 8000 10000
0
0.005
0.01
0.015
0.02
Wholetime blood input function with hybrid estimation
time (seconds)
activity

training part
exponential part
B
Figure 3.1: (A) the input function from 7 liver scan patients and (B) wholetime estima-
tionoftheinputfunctionforoneofthewholebodypatients.
3.4.2 RandomandScatteredEvents
The accuracy of randoms r(t) and scatters s(t) modeling is assessed by ensuring
dynamic tail matching of dynamic prompt rate function p(t) and r(t) + s(t). Each
5-min bed position is divided into 5 1-min segments. The integral of the r(t) + s(t)
in each segment is computed. The tail of random plus scatter matches quite well with
that of the prompts for all 5 segments (relative error 0.05 on average tail mismatch-
ing), demonstrating that the model accurately captures the dynamic change of scatters
and randoms. Figure 3.2 shows the 2 integrated scatter sinogram corresponding to the
second bed position in the 50th direct plane. It is clear that the scatters are of different
activityrate,althoughthespatialdistributionsareclose.
3.4.3 PatientMotionCompensation
Figure 3.3 shows the difference image of 2 MAP images reconstructed from the 2 list-
mode frames of one patient before and after registration with MATLAB registration
code. In this example the patient left the bed between the two scans presenting a more
difficult registration challenge than in the typical case. A coarse grid size was initial-
ized so that the large motion can be captured more rapidly. The bright and dark regions
77

0
0.05
0.1
0.15
0.2

0
0.05
0.1
0.15
Figure 3.2: The 2 integrated scatter sinogram corresponding to the second bed position
inthe50thdirectplane.
in Figure 3.3A, B shows misregistration in the heart, brain, liver and vertebrae before
registration. The post-registration results in Figure 3.3C, D clearly shows significantly
better alignment between the 2 frames. The same deformation field was then used in
computingthePatlakimages. Meanwhile,forsomepatientswithsmallmotionbetween
2 frames, we havefound that SVREG provides slightly better registration by tuning the
penaltytermcarefully.
Figure3.4showsPatlakinterceptimagewithdynamicrandomsandscatterscorrec-
tionandmotioncompensationagainsttheonewithoutdynamicrandomsscatterscorrec-
tion and motion compensation. Overall, these corrections produce noticeable improve-
ments in contrast and noise suppression as well as improved boundary definition as a
resultofcoregistration. Thispatienthadametalimplantintheleftlegwhichproducesa
clear artifact when no correction is applied. Also noticeable is that the bed shows up in
theupperimagewhenrandomsandscattersarenotcalibrated. Theseartifactsdisappear
inthefullycalibratedimages.
78

−3000
−2000
−1000
0
1000
2000
3000
4000
5000
A

−3000
−2000
−1000
0
1000
2000
3000
4000
5000
B

−4000
−2000
0
2000
4000
6000
8000
C

−4000
−2000
0
2000
4000
6000
8000
D
Figure 3.3: Difference of 2 MAP reconstructed images of frame 1 and frame 2 before
(A,B)andafterregistration(C,D)
B C D A
Figure 3.4: Patlak intercept images: (A, B) scatter correction only for the first frame
data,nomotioncorrection. (C,D)scatter,randomcorrection,motioncorrected.
79
3.4.4 Patlakv.s. SUVand%DSUV:
Examples of the reconstructed Patlak slope and intercept images with 2 corresponding
SUV images from a patient with detected tumor are shown in Figure 3.5. The Patlak
slopeK
i
imagehasthehighestcontrast with slightly increased noisecompared tothe 2
SUVimages.
Figure3.5: ExampleofreconstructedSUV1,SUV2,Patlakslopeandinterceptimages
We draw ROIs manually in the liver, thigh muscle and the Lumbar Spine (L2), and
themeanvaluesarecalculated. Table 3.1liststhemeanvalueofPatlakslope,intercept,
SUV1, SUV2, and the %DSUV values in the liver, thigh muscle, and Lumbar Spine
(L2)forallpatients. NotumorwasdetectedneartheROIsforallthepatients.
Wealsocalculatedthemeanandrelatives.d. oftheROImeansfromthe11patients
andlistedtheminTABLE3.1. Forall3ROIs,thes.d. acrosspatientsissmallestforthe
Patlakparameters,followedbystaticSUVs,andhighestfor%DSUV.Onlyfortheliver
ROI,s.d. ofPatlakinterceptissmallerthanK
i
.
Forliverandspine,theSUV
mean
attwoacquisitiontimediffersmorethan10%. The
%DSUVhashighlyvariabilitybecauseofthedivisionoperation.
80
Table 3.1: Patlak parameters, SUV and %DSUV values from liver, thigh muscle, and
lumbarspine(L2)ROIs
ID K
i
(mL/min/100g) V
B
(mL/g) SUV1
mean
SUV2
mean
%DSUV
mean
liver
1 0.43 50.82 2.220 2.114 -4.77%
2 0.38 50.16 2.781 2.508 -9.82%
3 0.46 55.02 2.496 2.157 -13.58%
4 0.49 47.16 2.900 2.686 -7.38%
5 0.31 43.92 3.512 2.618 -25.46%
6 0.50 43.38 3.623 3.137 -13.41%
7 0.41 42.36 2.665 2.458 -7.77%
8 0.38 45.18 2.537 2.212 -12.81%
9 0.46 41.70 2.681 2.311 -13.8%
10 0.32 50.34 2.174 1.694 -22.08%
11 0.47 27.40 2.466 2.266 -8.11%
Meansd 0.420.064 45.27.24 2.730.47 2.380.37 -12.6%6.31%
thighmuscle
1 0.14 14.28 0.702 0.684 -2.56%
2 0.13 12.36 0.816 0.807 -1.10%
3 0.18 11.46 0.672 0.663 -1.34%
4 0.16 10.68 0.788 0.812 3.05%
5 0.17 10.02 0.718 0.637 -11.28%
6 0.16 9.48 0.895 0.904 1.01%
7 0.14 14.58 0.803 0.780 -2.86%
8 0.18 11.16 0.541 0.627 15.90%
9 0.18 10.62 0.954 0.963 0.94%
10 0.15 10.08 0.467 0.462 -1.07%
11 0.136 9.42 0.497 0.554 11.5%
Meansd 0.160.019 11.291.78 0.710.16 0.720.151 1.10%7.3%
lumbarspine(L2)
1 0.59 19.5 1.343 1.386 3.20%
2 0.76 16.8 2.876 3.122 8.55%
3 0.73 12.1 1.856 2.002 7.87%
4 0.80 17.6 2.413 2.940 21.84%
5 0.88 18.1 2.216 2.434 9.84%
6 0.78 13.0 2.735 3.273 19.67%
7 0.60 18.1 2.069 2.536 22.57%
8 0.66 13.6 1.762 1.912 8.51%
9 0.76 13.7 2.289 2.615 14.24%
10 0.88 11.2 2.577 3.407 32.21%
11 0.34 10.86 2.125 2.346 10.4%
Meansd 0.730.11 15.612.92 2.200.45 2.540.62 14.45%8.61%
81
Figure 3.6 shows result from a 55-year-old patient diagnosed with colorectal liver
metastasis.. An invasive tumor was confirmed using ultrasound guided biopsy. As
shown quantitatively in Table 3.2, both PatlakK
i
and the intercept has higher contrast
(9.25,6.00)than2SUVimages(4.04,4.31).
Table 3.2: Contrast of breast tumor and liver tumor against background for SUV and
Patlakmethods
ROImean K
i
(mL/min/100g) V
d
(mL/g) SUV1
mean
SUV2
mean
breasttumorpatient
Tumor 0.35 16.75 2.16 2.24
BG 0.038 2.79 0.53 0.52
Contrast 9.25 6.00 4.04 4.31
livercancerpatient
Tumor 0.73 32.64 2.86 3.12
BG 0.48 27.78 2.50 2.35
Contrast 1.521 1.175 1.144 1.328
The 2 SUV images and Patlak images for a liver cancer patient is shown in Fig-
ure 3.7. Visually the Patlak slope image has the highest tumor-to-background contrast.
The quantitative analysis is shown in Table 3.2, where Patlak K
i
again has the highest
contrast(1.521)than2SUVimages(1.144and1.328respectively). However,thePatlak
interceptiswithcontrastlowerthanthesecondSUVimage.
3.5 Discussion
SUV is known to have large variation across different patients due to technological
factors such as scanner calibration, reconstruction algorithms, size of ROI, methods of
analysis(e.g.,maximumvs. mean)andbiologicalfactorsincludingbodysizemeasure-
ment,bloodglucoselevelanduptaketime. Inourstudy,thesefactorsarenotaccounted
for when computing the s.d. of SUV related measurement. It is fair to conclude that
the actual s.d. of SUV will be considerably larger than the one provided in the table.
82
B
C D
A
Figure3.6: CoronalviewofSUV1(A),SUV2(B),Patlakslope(C),andPatlakintercept
(D)images. Thearrowsindicatesthetumorintheleftbreast.
In addition, SUV is, at best, a semi-quantitative measurement for the metabolic rate
in the tissue. Accuracy of SUV depends critically on the assumption that the integral
of plasma FDG TAC is proportional to the ratio of injected dose and body weight (or
body surface area), which may be violated for example when patient is going through
chemotherapy and has impaired renal function ([95]). On the other hand, our prelimi-
naryresultssuggestthatPatlakslopeprovidesamorereliablesurrogatemeasurementof
the FDG metabolic rate in the tissue. Thus our method has the potential to fully utilize
thequantitativenatureof18F-FDGPET.
83
Figure3.7: CoronalviewofSUV1(A),SUV2(B),Patlakslope(C),andPatlakintercept
(D)images. Thearrowsindicatethetumorintheliver.
Patlak analysis has not been applied to routine clinical whole-body
18
F-FDG PET.
Onereasonisthecurrentwhole-bodyscanprotocolsdonotallowtheacquisitionofmul-
tiple time courses of the tissues as required by the traditional Patlak analysis. Another
difficulty is the blood input function. And finally the complication in the data correc-
tions including patient motion correction. In this paper, we developed techniques to
addresseachproblem,whichallowsPatlakanalysisinroutineclinicalPET/CTscans.
The whole-body protocol can be further optimized to reduce scanning time and
maintain image quality, and may also allow different data acquisition time for differ-
entbedpositionsinsomespecificapplications. The2imageregistrationapproacheswe
84
have can successfully compensate for the motion effect between 2 frames. Idealy for
patient with significant motion, the MATLAB toolbox can be used in big-ROI studies.
Otherwise the SVREG can be used to match more local structures with careful tuning.
Itiscriticalthattheregistrationresultbeinvestigatedwelltoproduceaccuratedeforma-
tion fields. In our work, we model the rate functions for randoms and scatters based on
thepromptrate. However,thismightnotbenecessaryforshortframes,becausetherates
of randoms and scatters will not change significantly in several minutes. In our study
wefoundthatthetemporalchangewithineach5-minbedpositionisverytrivialsothat
modeling the temporal function of randoms and scatters as constant in each 5-min bed
positionalsoworkswell.
The other advantage of the proposed wholebody Patlak estimation method is that
it is quantitative voxelwise, which is different from the ROI studies in conventional
approaches. The count independent and spatially invariant resolution property of the
proposedapproachalsoallowsmoreaccurateimageanalysisandlesiondetectiontasks.
Our estimated liver Patlak slope values are consistent with published values. Munk
([113]) reported 0.0037/min for the mean Patlak slope values for the pig liver. Iozzo
([87] ) obtained a mean value of Patlak slope for the fasting group and clamp group
0.0022/min and 0.0042/min respectively. The Patlak slope value is 0.003/min in the
study([114])and0.00564/minand0.00389/minfor2groupsofmicein([115]).
Ourpreliminaryresultsincludingabreastcancerpatientandapatient withcolorec-
tal liver metastasis suggest that Patlak method has less variability and higher tumor to
background contrast, which may result in higher sensitivity and specificity in lesion
detection than SUV. In future we would acquire more data and further investigate the
performance difference in lesion detection. specificity than SUV method. In future we
would acquire more data and further investigate the performance difference in lesion
detection.
85
3.6 Conclusion
We addressed several practical issues in wholebody Patlak image estimation. We intro-
duced a hybrid population-based method for estimating the blood input function. We
also demonstrated the effectiveness of our random and scatter correction and of motion
compensation. Wholebody patient studies show that we can produce Patlak images
of equivalent visual quality to static PET scans. This is in contrast to dual-time point
%DSUV values which do not produce good quality images. The data analysis sug-
gests that Patlak as a quantitative method yields less variation among population than
the semi-quantitative methods such as SUV or %DSUV. The contrast recovery results
in Table 3.2 also indicates that Patlak is more suitable in the task of lesion detection.
Clearlyfurtherstudiesarerequiredtodeterminewhetherthisapproachleadstoimprove
diagnosticpowerrelativetothealternatives.
86
Chapter4
DirectLoganparameterestimation
fromlistmodedataforreversibletracer
4.1 Abstract
WeproposedaneffectivemethodtoestimateLoganparametersfromlistmodePETdata
based on our previous Patlak estimation framework. A static sinogram from time of
injection to the steady state is incorporated to improve the ill-conditioning of the prob-
lem. Cramer-Rao analysis and Monte-Carlo simulation demonstrate that our method
significantly improves the conditioning of estimation and provides more robust result.
Our method is also shown to be advantageous over the sinogram based approach in
providing higher contrast recovery and smaller standard deviation in Monte-Carlo sim-
ulation. In future more quantitative analysis will be carried out to better compare the
performanceofourmethodwiththesinogrambasedapproach.
4.2 Introduction
In dynamic imaging, the data is usually acquired continuously from the time the tracer
is administered and then fit to a quantitative kinetic model [29]. These models allow
us to represent tracer uptake in terms of quantitative rate parameters which have the
87
potential for improved tumor characterization compared to the more widely used semi-
quantitativemeasuressuchastheSUV[30]. IthasbeenshownthatdynamicPETimag-
ingimproveslesiondetectionsensitivityandspecificityforcertaintumors[31][22].
The traditional approach in dynamic PET imaging requires multiple static image
reconstructionsforeachframesfollowedbykineticanalysisappliedtothereconstructed
images [116]. This approach is not ideal for voxelwise parametric images and suffers
from high noise, as the frames are reconstructed independently. Besides, the recon-
struction for multiple frames is time consuming and matching image resolution for the
reconstructed images for quantitative analysis is difficult. An alternative approach is to
directly fit the data to a kinetic model at each voxel in the image. This idea was first
proposed by Snyder [43]. Since then several maximum likelihood [54], [55], [56] and
penalized ML [57], [58] methods have been described for this purpose. However, the
rate parameters in these kinetic models are nonlinearly related to the data and result in
nonconvex optimization problems. Furthermore, the relatively high dimensionality of
the search (four per voxel for the two compartment model typically used for FDG) can
resultininstabilityoftheestimatesandhighcomputationcost,limitingtheuseofthese
directestimators.
Severalsimplifiedgraphicalmethodshavebeenproposedasmorerobustalternatives
for the kinetic analysis, including Patlak graphical analysis [44] for irreversible tracers
andLogangraphicalanalysis[49]forreversibletracers. Previously,severaldirectPatlak
parametric reconstruction algorithms have been developed for irreversible tracers such
as
18
F-FDG (sinogram based approach [59] and listmode based approach [42]). The
focus in this paper is on direct estimation of Logan parameters. To allow the activity
represented with weighted basis functions for the purpose of performing direct Logan
estimation from data, Zhou proposed the relative equilibrium (RE) method [117] as a
surrogate of original Logan equation. Wong [118] has compared the Logan plot and
88
the RE plot and concluded that the RE plot provides less noisy parametric images and
gives consistent and reliable regional DVR (distribution volume ratio) estimates. In
direct Logan parameter estimation involving reversible binding tracers, only sinogram
based approach has been described [119]. This method uses the relative equilibrium
(RE) graphical analysis formulation [117] with a generalized AB-EM algorithm [120].
However, the correlation between the dependent variables (cumulated activity for each
time point from time 0min) is not considered. Therefore the algorithm does not reflect
the true likelihood of the data. Furthermore, the timing information for each count is
lost in the histogramming operation. Wang et al. [121] tried to model the activity rate
function by taking the derivative of the RE equation and use listmode data to perform
estimation. However,he did not implement the method or study the conditioning of the
approach. As we will show in the method section below, simply using listmode based
estimationmethodafterthesteadystateisaveryill-conditionedproblem. Therefore,the
directestimationforLoganparameterswithlistmodedatahasnotbeenwellestablished.
OurmethodisdifferentfrompreviousworkinthatweusethelistmodePETdatato
exploit more from the timing information of counts, rather than using sinogram based
method to bin counts arriving in each frame to a single sinogram. Besides, the method
we derived strictly follows the true likelihood of the data model, which is distinct from
thesinogrambasedmethodwhichdoesnotconsiderthecorrelationbetweendependent
variables. More importantly, we distinguish our method from the method proposed by
Wang [121] by using an extra sinogram from beginning time to the steady state, which
effectivelyimprovestheill-conditioning.
89
4.3 Materialsandmethods
4.3.1 DirectLoganparameterEstimation
The Logan model provides an approximate solution to the two compartment kinetic
model for reversible tracers in the steady state. The simplified Logan graphical model
afterthesteadystatecanberepresentedbythefollowingequation:
R
t
0
()d
(t)
=DV
R
t
0
C()d
(t)
+q (4.1)
where DV is the Logan slope (distribution volume) and q is the Logan intercept
(bindingpotential). (t)istheactivityrateandC(t)isthebloodinputfunction. Follow-
ing the RE equation in [117], we can approximate the Logan equation after the relative
equilibrium state as in equation( 4.2). In the following we use the RE equation in our
directestimationalgorithm.
R
t
0
()d
C(t)
=DV
R
t
0
C()d
C(t)
+q (4.2)
Multipling equation( 4.2) byC(t) at both hand side and taking derivative yields the
representationoftheratefunctionafterthesteadystateinthemeasuredsinogramspace
atlineofresponse(LOR) i:
(t) =C(t)+qC
0
(t) (4.3)
Thisindicatesthatwecanrepresenttheratefunction(t)atvoxeljaftersteadystate
tT
0
asalinearcombinationoftwobasisfunctionsB
1
(t)andB
2
(t):

j
(t) =
2
X
l=1
!
jl
B
l
(t) (4.4)
90
where!
j1
istheLoganslope, !
j2
istheLoganintercept, B
1
(t) =C(t)istheblood
inputfunction,andthat B
2
(t) =C
0
(t)isthederivativeoftheinputfunction.
Theratefunctionateachdetectorpairafterthesteadystatecanthenbewrittenas:

i
(t) =e
t=
nv
X
j=1
2
X
l=1
p
ij
!
jl
B
l
(t)+r
i
(t)+s
i
(t) (4.5)
where the exponential term accounts for radioactive decay of the tracer with half-
life Ln(2); p
ij
is the probability of an event at voxel j being detected at detector pair
i; and n
v
is the number of voxels. In the following we use a factored system model
for p
ij
based on a product of attenuation factors, normalization, geometric response,
and a Monte Carlo model of the detector response as described in [9] . The additional
terms r
i
(t) and s
i
(t) denote, respectively, the randoms and scattered coincidence rate
functions.
Assuming we have continuous list mode data over 2 time interval [t
1
;t
2
], and that
the arrival times in the list mode data follow an inhomogeneous Poisson model, the
continuoustimepenalizedlog-likelihoodfunctionofeventarrivaltimesisgivenby:
F(W) =
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)+
np
X
i=1
(
Z
t
2
t
1

i
(t)dt)
+1=2!
T
R
S
!
(4.6)
wherea
ik
denotesthearrivaltimeofthekthphotonatdetectorpairiintheintervals
[t
1
;t
2
], x
i
is the number of events detected in LOR i, and n
p
is the total number of
LORs. R
S
computes the squared pairwise difference between each voxel and its 26
nearest neighbors for both the slope and intercept parameters. We compute solutions
usingthe4Dincrementalgradientmethodasdescribein[62].
91
TheaboveLoganestimationframeworkissimilarwiththePatlakestimationframe-
work we previously proposed [42]. However, the estimation for Logan parameters is
a more ill-conditioned problem than Patlak parameter estimation, because of the simi-
larity of basis functions after the steady state. We will investigate in this problem and
provideaneffectivesolutioninthenextsection.
4.3.2 Overcomingill-conditioning
The blood input function can usually be well modeled as the combination of multiple
exponentialfunctions[108]. Aftersufficienttimeelapses(usuallyafterthesteadystate),
it is sufficiently accurate to approximate it with a single exponential function, because
of the different decay rates for exponential functions. Assuming that approximately we
have C(t) = e
t
for t > T

, where T

is the steady state start time, we can obtain
C
0
(t)as
C
0
(t) =e
t
(4.7)
ClearlyC
0
(t) =C(t) is a scaled version of C(t) aftert > T

under the assump-
tionofsingleexponentialfunction. Thereforeequation(4.3)canberewrittenas:
(t) = (q)C(t) (4.8)
whichisanill-conditioned problemand that it will be verydifficultto estimateDV
andq separately.
We try to improve the ill-condition of this estimation problem by utilizing extra
information before the steady state. We have shown that the log-likelihood function for
thelistmodedatain [T

;T]canbewrittenas[42]
92
F
l
(W) =
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)+
np
X
i=1
(
Z
T
T

i
(t)dt) (4.9)
For the events collected from injection to the steady state start time T

, we use a
Poisson model for the cumulated sinogram data. The log-likelihood for sinogram data
[0;T

]canbewrittenas:
F
s
(W) =
np
X
i=1
N
i
log
Z
T

0

i
(t)dt+
np
X
i=1
(
Z
T

0

i
(t)dt) (4.10)
So the likelihood function for dynamic PET data acquired between injection and
someT >T

istheadditionofF
l
(W)andF
s
(W):
F(W) =
np
X
i=1
N
i
log
Z
T

0

i
(t)dt
np
X
i=1
x
i
X
k=1
log
i
(a
ik
)
+
np
X
i=1
(
Z
T
0

i
(t)dt)
(4.11)
Note here that the listmode data x
i
is constrained to listmode data after the steady
stateonly,excludingthecumulatedsinogramfromtimet = 0tillthesteadystateT

.
We will show the Fisher information matrix of equation( 4.9) and equation( 4.11)
to show the improvement of conditioning by adding the sinogram from injection to the
steady state. From the log-likelihood function ( 4.9), we derived the following expres-
sion for the second derivative with a simplified case which only has one voxel and we
directobservethePoissoncounts:
93
E[
@L
2
(W)
@
2
] =
x
i
X
k=1
C(a
k
)
2
(a
k
)
2
(4.12)
E[
@L
2
(W)
@q
2
] =
x
i
X
k=1
C
0
(a
k
)
2
(a
k
)
2
(4.13)
E[
@L
2
(W)
@@q
] =
x
i
X
k=1
C(a
k
)C
0
(a
k
)
(a
k
)
2
(4.14)
wherei = 1istheonlydetectorpairsdirectlyobservercountsfromthesinglevoxel.
With the approximated binmode cost function [42], we can approximate equation( 4.12
to 4.14)asfollows:
E[
@L
2
(W)
@
2
] =
(
R
T
T

C(t)dt)
2
R
T
T

(t)dt
(4.15)
E[
@L
2
(W)
@q
2
] =
(
R
T
T

C
0
(t)dt)
2
R
T
T

(t)dt
(4.16)
E[
@L
2
(W)
@@q
] =
R
T
T

C(t)dt
R
T
T

C
0
(t)dt
R
T
T

(t)dt
(4.17)
Clearly, the product of the diagonal of this 2 by 2 FIM equals the product of the
off-diagonal elements, which yields a highly ill-conditioned Fisher information matrix.
However,wefindthatusingequation(4.11)thederivativesarechangedas:
E[
@L
2
(W)
@
2
] =
[
R
T

0
C(t)dt]
2
R
T

0
(t)dt

[
R
T
T

C(t)dt]
2
R
T
T

(t)dt
(4.18)
E[
@L
2
(W)
@q
2
] =
[C(T

)]
2
R
T

0
(t)dt

[
R
T
T

C
0
(t)dt]
2
R
T
T

(t)dt
(4.19)
94
E[
@L
2
(W)
@@q
] =
C(T

)
R
T

0
C(t)dt
R
T

0
(t)dt

R
T
T

C(t)dt
R
T
T

C
0
(t)dt
R
T
T

(t)dt
(4.20)
Because of the interaction of the term introduced by data from [0;T

], the Fisher
information matrix can be easily shown to be not singular. This can also be explained
by that fact that the blood input function before the steady state does not follow an
exponentialfunctionsothatitbreaksdownthelineardependencyofcolumnsintheold
FIM ( 4.12 to 4.14). Therefore ( 4.11) can be used for more stable Logan parameter
estimation than listmod-only and steady-state-only algorithm. We will show this in the
resultsection.
4.4 Results
In this section, we first investigate in the validity of Logan plot and RE plot after the
steadystateforvariouskineticparametersinsimulation. Wethenshowtheimprovement
in overcoming ill-conditioning of the problem using our method as in equation( 4.11)
instead of the listmode-only approach as in equation( 4.9) with Cramer-Rao analysis
and Monte-Carlo simulation. More importantly, we compare the performance of our
methodwiththesinogrambasedapproachasin[119],andshowthatourmethodyields
higher contrast recovery and lower variance for small tumor ROI against the back-
ground. Finally, the clinical example reconstructed with our method and the sinogram
basedmethodareshown.
The studies reported in the following subsections used the kinetic parameters listed
inTable 4.1,whichwerefoundintheliterature.
95
SetID K1 k2 k3 k4
1 .0918 0.4484 0.141 0.1363
2 0.303 0.391 2.27 0.572
3 0.325 0.437 1.82 0.277
4 0.105 0.305 0.07 0.1
5 0.347 0.065 0.007 0.077
Table 4.1: reversible kinetic parameters for several brain or body regions reported in
literature
4.4.1 LoganandREplotvalidation
We first investigate the validation of the Logan and RE plot for the kinetic parameters
in our simulation. The Logan equation is valid when C
2
(t)=(C
1
(t) +C
2
(t)) becomes
constant. The RE equation is valid after C
1
(t)=C
p
(t) and C
2
(t)=C
p
(t) are both con-
stant. Figure4.1showsthedynamicchangeofC
2
(t)=(C
1
(t)+C
2
(t))andC
1
(t)=C
p
(t),
C
2
(t)=C
p
(t) simulated with kinetics as in the first row of Table 4.1. From the plot, the
termC
2
(t)=(C
1
(t)+C
2
(t)) approaches a constant after 1 hour, whileC
1
(t)=C
p
(t) and
C
2
(t)=C
p
(t)becomeaconstantroughlyatthesametime. Weobservethesamebehavior
for other kinetic parameters in Table 4.1. Therefore it is reasonable to assume that the
steady state is reached after 1 hour and the RE equation can be used to perform Logan
parameterestimationasanapproximatelyequivalentalternative.
4.4.2 Cramer-Raoanalysis
To gain insight into the potential performance of our approach in estimating Logan
parameters, we computed Cramer Rao lower bounds (CRLBs) as a measurement of
the estimation variance. Same as that for the Patlak analysis, we explore the bounds
only for the simplified case of a single LOR and a single voxel. In that case, a spatial
penalty is not required to obtain a stable estimate so that we might expect that bounds
computed from the un-penalized log likelihood will reflect the actual variances seen in
96
0 1000 2000 3000 4000 5000 6000 7000 8000
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time (secs)

Logan C2(t)/(C1(t)+C2(t))
RE C1(t)/bld(t)
RE C2(t)/bld(t)
Figure4.1: ThedynamicchangeofC
2
(t)=(C
1
(t)+C
2
(t))andC
1
(t)=C
p
(t),C
2
(t)=C
p
(t)
asafunctionoftime.
our estimators. We aim at comparing the variance obtained with listmode data after
the steady state against our modified algorithm with additional data before the steady
state to demonstrate that the modified algorithm greatly improves the conditioning of
theestimationproblem.
WeusedtheFIMexpressionstocomputeCramer-Raoboundsforthekineticparam-
etersandcomparedthesetomaximumlikelihood(ML)parametersestimatedfromdata
simulatedforeachmodel. Thedataweresimulatedbygeneratingarrivaltimesforpho-
tons corresponding to a time inhomogeneous Poisson process defined by the TAC for
eachparametersetforarealisticinputfunction. Figure4.2showsthebloodinputfunc-
tionandthetime-activity-curvefortheregionwithkineticparametersasinthefirstrow
inTable 4.1.
To compute the bounds, the mean random and scatter rates were assumed known.
We then retained those arrival times corresponding to 60-140min allowing estimation
with 80-min listmode data after the steady state. Then for comparison, we only use
97
0 1000 2000 3000 4000 5000 6000 7000 8000
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
time (sec)
activity

blood
Brain region
Figure 4.2: The blood input function and the time-activity-curve for the tumor region
forthefirstrowinTable 4.1.
60-80minslistmodedata,butmeanwhileincludethesinogramdatafrom0-60min. This
makes a fair comparison between the traditional listmode based method and our modi-
fiedmethod,asthetotalscantimeisequalto1hourforbothcases.
Table 4.2 shows the standard deviation computed from Cramer Rao bound for the
listmode-onlymethodandourmodifiedmethodforall5kineticparametersinTable4.1,
where we use the same blood input function. The listmode-only estimation method
gives huge variation because the extremely high correlation of 2 basis functions after
the steady state prevents the separation of the 2 Logan parameters. On the contrary,
our modified method yields a much smaller relative standard deviation in all cases. We
alsotunethedurationoflistmodeacquisitionandthelistmodestartingtimeandfindthat
uniformlythestandarddeviationofourmethodismuchsmaller,whilethelistmode-only
basedalgorithmconsistentlydoesnotproducereliableresults.
Meanwhile,toshowtheadvantageofourestimationmethodoverthesinogrambased
approach, we investigate in a one-pixel study and compare the estimation accuracy in
98
Kinetics 1 2 3 4 5
slope(withsinogram) 3.4% 1.1% 0.9% 2.8% 0.9%
slope(nosinogram) 1e8% 3.9e7% 7e7% 1.1e8% 5e7%
inter(withsinogram) 14.2% 15.2% 13.2% 13.7% 14.0%
inter(nosinogram) 5e8% 6.7e8% 1.3e9% 6.8e8% 9.6e8%
Table 4.2: Comparison of the relative s.d. (s.d. divided by the mean) of Logan param-
eters obtained from the CR bound with and without cumulated sinogram data from
time=0-mintothesteadystate60-min
termsofvarianceinMonte-Carlosimulation. Thisisjustifiedbythefactthatourmethod
differs from the sinogram based approach mainly in the handling of timing information
ofdatainsteadofspatialinformation. WecomputethestandarddeviationofLoganslope
and intercept for our algorithm and the sinogram based approach. We allow the two
methodstousethewholetimeseriesofdatafromt=0mintot=120min,andsett=60min
asthestartingpointofusinglistmodedataforourmethod. Weuseawindowof4minto
histogramsinogramframesforperformingsinogrambasedestimation. Figure4.3shows
thestandarddeviationforthe2methodscomparedwiththeCramerRaoresults. Clearly
our method yields smaller standard deviation across the Monte-Carlo simulation. This
suggests that our method is capable of predicting Logan parameters with less variation.
Wethinkthebenefitcomesfromthecorrectmodelingofdatastatisticsandtheobviation
ofcorrelatednoisefromcorrelateddata.
4.4.3 Monte-Carlosimulation
We simulated a small scale 3D PET system (diameter: 148.4mm, detector size:
2.423mm x 2.423mm x 2.025mm; number of rings: 4) with a total of 13 2D sino-
grams each with 84 angles of view by 96 radial LORs. A uniform cylinder phan-
tom of diameter 31.4mm is centered in the scanner and contains 5 cylinders (tumors)
of diameter 1.0, 1.8, 2.6, 3.4, and 4.2mm as shown in Fig. 4.4(a). We simulated
99
A B C D E
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
s.d.
s.d. of estimated Logan slope

CR s.d.
our s.d.
sinogram s.d.
A B C D E
0
200
400
600
800
1000
1200
s.d.
s.d. of estimated Logan intercept

CR s.d.
our s.d.
sinogram s.d.
Figure4.3: standarddeviationcomputedfromCramerRaoanalysisandsimulation,both
withourmethodandthesinogrammethod.
time activity curves for 2 different regions in the brain using 2 sets of kinetic param-
eters: K
1
= 0:303;k
2
= 0:396;k
3
= 2:270;k
4
= 0:572 for simulating tumor and
K
1
= 0:092;k
2
= 0:448;k
3
= 0:141;k
4
= 0:136 for simulating the background. The
blood input function, tumor and background TACs are shown in Fig. 4.5. The TACs
foreachsinogramelementarecomputedbyforwardprojectionthroughasystemmatrix
100
generated for the scanner using the model described in [9]. Each of these TACs repre-
sentsthetime-inhomogeneousPoissonprocessesthatdetermineemissionratesforeach
LOR. List-mode data are then generated by sampling from these TACs. The total list-
mode events in the 140-min acquisition time is 6.6M. The dynamics of randoms and
scatters are also modeled. We include scatter proportional to the true rate and randoms
proportion to the square of the true rate, each scaled to contribute at 20% of the mean
truecountswhenaveragedovertheentireTAC.
In the results presented here, we window the data to retain only those events corre-
spondingto60-160minforthelistmode-onlyalgorithm. Incomparisonforourmodified
methodweusethe60-100minlistmodedataandaddthe0-60minsinogramdata. Forthe
sinogrambasedapproach,weuses4cumulatedsinograms,withstartingtime0-minand
acuquisition time ending at 60, 80, 100, 120-min respectively. We generated a total of
100trialsforourMonteCarlosimulations. Inthefollowingweshowthattheestimation
from the modified method has much better accuracy and much less variation than the
listmode-only method, and also outperforms the sinogram based method with higher
contrast recovery and lower variance. We used the modified 4D incremental gradient
(4DIG) algorithm [63] to reconstruct the Logan slope and intercept images. The mean
of the Logan slope images with our method and the sinogram based approach over 100
trialsisshowninFig.4.4(b)Fig.4.4(c).
Note here that we do not include the reconstruction from the listmode-only method
results,butwehavedonethesamesimulationwiththeonlydifferenceinthatwedidnot
includesinogramdatafrominjectiontothesteadystate. Thenumericalresultisthatthe
algorithm did not even converge, which turns out to be not useful quantitatively. This
is consistent with our Cramer Rao analysis in the previous subsection where we find
that the relative variance of Logan estimates are on the scale of e8 for listmode-only
algorithm.
101

0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5

0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2

0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 4.4: (a) Uniform cylinder (background) with five cylinders (tumor) of different
sizes. (b) The mean (over 100 Monte Carlo trials) of the Logan slope image recon-
structed with our method with TACs as in Fig. 4.5. (c) The mean (over 100 Monte
Carlo trials) of the Logan slope image reconstructed with the sinogram based approach
withTACsasinFig. 4.5.
102
0 2000 4000 6000 8000
0
1
2
3
4
5
6
x 10
−3
time (secs)
activity

tumor
background
input function
Figure 4.5: Simulated time-activity-curves for tumor, background, and the blood input
function.
Thedifferenceofour method and the sinogram based approach differmost in small
ROIs. Figure 4.6 shows the profile through the smallest tumor rod for the mean Logan
slope images computed with our method and the sinogram based method. Clearly the
contrastrecoveryishigherforourmethodforthesmallesttumorrod. Wealsocompute
the variance of the mean values in the 5 tumor rod ROIs from 100 Monte-Carlo simu-
lations and find that the variance is smaller for our method for all ROIs and the bias is
lower, as in Table 4.3. This demonstrates that our method is more advantageous than
the sinogram based approach in producing more accurate and robust Logan parametric
images.
ROI 1 2 3 4 5
ourmethod 2.990.12 3.380.039 3.480.036 3.630.02 3.670.02
sinobased 2.770.13 3.250.043 3.370.039 3.450.03 3.510.03
Table 4.3: mean and s.d. of the 5 tumor rod ROIs in 100 Monte-Carlo simulation for
ourmethodandthesinogrambasedmethod
103
60 80 100 120 140 160 180 200
0
0.5
1
1.5
2
2.5
3
3.5
4
voxel index
absolute Logan slope value

our method
sino based
Figure 4.6: profile through the smallest tumor rod for the mean Logan slope images
computedwithourmethodandthesinogrambasedmethod.
4.5 Discussion
We proposed a new approach to directly estimate the Logan parameters from listmode
PET data. The advantage of our method over the sinogram based approach is that (i)
the true likelihood of the data is modeled; (ii) the rate function modeling may improve
on image estimation. On the contrary, the sinogram based approach uses cumulated
sinogram data in all the dependent variables so that the model will suffer from corre-
lated noise. The advantage of our method over the listmode-only approach is that we
significantlyimprovetheconditioningoftheproblembyaddingastaticsinogram.
Estimation of the blood input function is important in obtaining an accurate esti-
mateoftheLoganparameters. Inoursimulationweusedaknownbloodinputfunction.
Usually the Logan plot in clinic settings is used for reversible tracers in the brain scan.
Therefore one of the possibilities is to allow the scan cover the carotid artery and seg-
ment it to estimate the blood input function. The other way is to use reference region
104
in the brain with no specific binding to the tracer and with a known k
2
. This will be a
futureexploration.
From the theoretical derivation and the computational results, we have shown that
our modified method can greatly improve the ill-conditioning of the Logan estimation
problem when listmode data is used. Meanwhile, we haven not yet thoroughly studied
the relationship between the performance with the protocol parameters. Given a fixed
subsetoflistmodedataafterthesteadystate,thechoiceofdurationofthesinogramdata
at the beginning is not yet optimal in maximizing the estimation performance. Push-
ing the sinogram collection time to before the steady state may result in invalidity of
the Logan or RE equation, and extending it long may result in a less well-conditioned
inverse problem if the listmode data is short, because the sinogram data will dominate
inthedatacollection. Wewillalsoinvestigatemoreinthisissue.
The disadvantage, and perhaps all other current methods in performing Logan anal-
ysis, is that the data collection has to start from the time of injection. This is different
from the Patlak approach, which does not need continuous data acquisition from the
time of injection. Practically, however, the Logan equation for performing graphical
analysis is usually approximately valid even before the steady state. Therefore we may
reducethesinogramdatacollectiontimeatthebeginningofthescan,althoughitisstill
requiredthatthepatientstaysinthescannerfromtheinjectiontime. Ourfinalgoalisto
develop a method which only needs short-time partial data after the steady state, which
willbemorefavorableinclinicalusage.
105
Chapter5
Jointestimationofactivityand
attenuationforcontrastenhanced
PET/CT
5.1 Introduction
AttenuationinPETimagingisthelossoftruecoincidenceeventsbecauseofabsorption
inthebodyortheComptonscatteroutofthefield-of-view. Thelossoftrueeventsdueto
attenuationinPETcanrangeupto5095%,especiallyforlargepatients. Withoutapply-
ing attenuation correction, image reconstruction may suffer from artifacts or distortion,
suchasshowingprominentactivityatbodysurfaceandlowactivityindeeperstructures.
Therefore, attenuation correction is critical for quantitative PET image reconstruction
andanalysis.
Before the development of PET/CT scanners, the attenuation was obtained with a
transmission source (e.g.
68
Ge) directly measuring the photon attenuation. The addi-
tional transmission scan usually took 3-10 min or more for each bed position. This
significantly increases the time required for oncologic applications with whole-body
procotols(typically5-7bedpositions)andlimitstheusageofthetechnique[65].
AnothersolutionforattenuationestimationistocombineaCTscannerwiththePET
scanner,andusethecoregisteredCTimagetoestimatetheattenuationcoefficientsinthe
PETenergyrange(511keV)withaprecomputedmappingtable[122]. Severalmapping
106
methods have been proposed for this purpose. Kinahan performed tissue segmentation
first and then assigned predefined values to different tissues [68]. GE Healthcare uses
a bilinear transformation combining the segmentation and scaling methods [67] [123].
Compared with the rotating source approach, the CT-based attenuation estimation is
muchfasterbecauseoftheshortacquisitiontimeofCTscans. Italsoimprovesthepre-
cisionoftheattenuationcorrection[124]. Intheimagequantificationaspect,ithasbeen
shownthatthetransmissionsourcescanshavinghighernoise,lowerbiasbutwillsuffer
from extra bias from the post-injection transmission scan, and CT scans have negligi-
ble noise [66]. Burger [67] compared the attenuation obtained from both methods and
foundthatforhumanbeingstheattenuationobtainedbythese2methodsareessentially
equivalent.
Recent PETCT scans are sometimes performed with contrast enhanced agent in CT
(oral or inject depends on the specific application) so that some blood-rich organs and
tissues (e.g. liver, heart, spleen) will show more detailed structures and have a better
contrast from the background. This enhancement has been shown to benefit in oncol-
ogy. Antoch [69] concluded that the contrast agent in contrast CT can add value to
lesion detection and characterization. Fletcher [70] reported that contrast enhanced CT
colonography is a promising method for detecting local recurrence, metachronous dis-
ease, and distant metastases in patients with prior invasivecolorectal carcinoma. Oliver
[71]foundthatintheapplicationofmultiplephasecontrastCT,usingarterialphaseand
portalvenousphaseimagesisbetterthanusingtheconventionalandportalvenousphase
images to reveal more hepatocellular carcinoma lesions. However, negative reports of
using contrast-enhanced CT also exist. Schaefer [72] has shown that the noncontrast
CT is more sensitive and specific than contrast-enhanced CT for evaluation of lymph
node and organ involvement, especially regarding exclusion of disease, in patients with
Hodgkindiseaseandhigh-gradenon-Hodgkinlymphomathan.
107
Thecontrastagentcausespotentialproblemforaccurateattenuationestimation. The
contrastagentincreasestheHounsfieldunitoftheCTimageashighasseveralhundred
in some blood rich regions such as the spleen or aorta. If the conventional CT-to-umap
(umapistheattenuationcoefficientsmap)mappingisused,theseenhancedregionswill
beassignedattenuationcoefficientsmuchhigherthanthevaluescorrespondingtowater
and normal soft tissue. However, photons penetrating through these regions containing
the contrast agent in the PET scan is usually attenuated similar as penetrating through
water,andconsequentlyshouldbegivenamuchlowerumapvalue. Thisapparentlywill
cause overestimation of the umap. Although there are researchers reporting the negli-
gible error caused by contrast in PET quantification, (for example, Yau [75] reported
that contrast enhanced CT does not significantly increase the SUV
max
), many other
researchersreportedoverestimatedSUVorattenuationcoefficientswhencontrastagent
is used [73], [74]. Therefore, it is necessary to process contrast enhanced PET/CT data
differently than conventional PET/CT to compute accurate attenuation coefficients for
quantitativePETimageanalysis.
There are generally two main directions in PET attenuation estimation when the
conventional CT scan and the transmission scan are not available. Several researchers
have developed segmentation based approaches to estimate attenuation. Their methods
basically apply segmentation algorithms to segment regions with constant attenuation
coefficients and assign predefined values [76] [77]. However, this method does not
reflect the individual character, and may incur large error for specific patients. For
example,thepatientswithmetalimplantshouldhavehighattenuationcoefficientsinthe
implantedmetalregion,yetthisisnotconsideredwiththesegmentationbasedapproach.
Theaccuracyofsegmentationandregistrationalsoaffectthequantitativeresults.
108
Severalotherresearchershaveworkedonjointestimationofactivityandattenuation
(also known as MLAA) from PET emission data [78, 79, 80]. The theoretical founda-
tionistheanalyticalconsistencyconditionandthediscreteconsistencycondition. How-
ever,thisjointestimationinthemaximum-likelihoodorpenalizedmaximum-likelihood
framework is a nonconvex problem, and results in ”cross-talk” between activity and
attenuation, or a local maxima of the cost function. This can be explained by the fact
thattheFisherinformationmatrixfornon-TOFPETdataisill-conditionedsothereisno
global maxima [81]. With some information as prior, however, the conditioning of the
problemcanbeimprovedtoproduceusefullocalmaximathatmaybeclosetotheglobal
maxima [82, 83]. In particular, Nuyts et al [83] proposed a jointly iterative algorithm
for estimating emission and attenuation simultaneously. The prior that the intensity of
attenuation map should be centered at several narrow peaks indicating normal tissue
or air or bone is modeled into the optimization procedure providing a local optimiza-
tion solution that is useful in practice. Time-of-Fight (TOF) PET data provides more
information than conventional nonTOF PET data. Interestingly, Defrise showed that in
time-of-fight PET, the attenuation sinogram is determined by the emission data except
for a constant [84]. Rezaei also studied the joint estimation problem in Time-of-Flight
PET and theoretically proved the uniqueness of global maxima up to a constant for the
attenuation [81]. This property of TOFPET makes MLAA expecially useful in many
applications. Boellaard showed that MLAA can improve quantitative accuracy of PET
imageinPET/MRstudies[85]. Hamillfound thatTOF-MLAAreducedmotion-related
PETimageartifactsinthelowerlungsinaphantomstudy[86].
In this chapter, we follow the joint estimation approach proposed by Nuyts et al
[83],andsubstitutetheiterativeOSEM-basedupdateforemissionsourcewiththeMAP-
based update. This ensures that the noise does not propagate after multiple iterations.
109
Our application of this joint estimation is for clinical nonTOF PETCT scans with con-
trastagent. ThereforewealsoincorporatethecontrastenhancedCTimageintoourjoint
estimation problem as prior information. We show in simulation that with the contrast
enhanced CT and further prior information of the range of CT values, the attenuation
coefficients could be estimated within a certain tolerable range and the error caused to
thefinalPETimageiskeptwithinanarrowrange.
5.2 Materialsandmethods
5.2.1 BilinearmappingfromCTtoattenuationcoefficient
Thecontrastagentonlyaffectslimitedregionsinsidethepatientbody,primarilyinblood
richregionssuchasaorta,spleen,andliver. Theseorgansandtissueshaveintermediate
CT values (30-70HU) in the CT histogram. Other regions such as bone or lung are
almost unaffected by the contrast agent. In the result section, we will show that the
range of CT values affected by the contrast agent in the contrast CT is concentrated
in [70 350]HU, by carefully investigating a 2D histogram for noncontrast and contrast
CT. This is a very important prior leading to useful solution for the non-convex joint
estimationproblem.
People have used the same approach for conventional CT-to-umap mapping in con-
verting contrast CT image to attenuation map, with the difference in that the bilinear
mapping curve may be unchanged or altered. As a simple implementation, we have
generatedamodifiedcurveformappingCTvaluetoattenuationcoefficientforcontrast
CT. Figure 5.1 shows the standard mapping curve from Siemens for noncontrast CT
(blue)andoursimpleimplementedcurveforcontrastenhancedCT(red). Noteherewe
assume that the Hounsfield unit range [70 350] HU are all soft tissues enhanced by the
110
contrast agent, and therefore mapping this range to the same umap value. This is the
mainreasonwhyusingmodifiedsimplecurvecausesbiasinattenuationestimation.
Figure 5.1: the standard converting curve from Siemens (blue) and our implemented
curve for contrast enhanced CT (red color indicating the change made upon original
mappingcurve).
5.2.2 Jointestimationofemissionandattenuation
Wecanrepresentthelog-likelihoodfunctionofPETdataasfollows:
L(;u) =
X
i
r
i
+y
i
logr
i
(5.1)
where is the emission activity in image domain, u is the attenuation coefficient,
i is the sinogram index, y
i
is the observed sinogram, and r
i
is the expected number of
counts. Wecanalsowritethefollowingequationstodescribethedataattenuation:
111
b
i
=
X
j
p
ij

j
(5.2)
a
i
=exp(
X
j
l
ij
u
j
) (5.3)
r
i
=a
i
b
i
(5.4)
Hereb
i
denotes the activity in LOR i before attenuation, a
i
denotes the attenuation
factorforLORi. l
ij
strictlyfollowthelineintegralmodelbutinpracticeitcanbetreated
asequalwithp
ij
.
Whentheattenuationcoefficientisfixed,itcanbereadilyshownthattheupdatefor
emissioncanbecarriedoutwiththefollowingMLEMalgorithm,asintroducedin[83].

new
j
=

j
P
i
a
i
p
ij
X
i
p
ij
y
i
b
i
(5.5)
However considering the propagated noise when the iteration number increases, we
use the conjugate gradient of MAP cost function as the updates when fixing the attenu-
ation.
For updates of the attenuation when fixing the emission source, the following for-
mulaisadopted:
u
new
j
=u
j
+
P
i
p
ij
a
i
b
i

P
i
p
ij
y
i
+M
0
j
(u)
N
P
i
p
ij
a
i
b
i
M
00
j
(u)
(5.6)
whereM(u) is the prior function on the attenuation map and can vary dramatically
over different applications. In our study where the contrast CT is provided, the penalty
M(u)isdesignedasfollows:
112
M(u) =
8
<
:
jjuu
0
jj
2
jju0:096jj
2
70HUCT350HU
1 otherwise
(5.7)
This design aims at keeping the regions that are not affected by the contrast agent
to maintain original attenuation coefficient, but meanwhile imposing a penalty for the
contrast-affected regions in CT image between 70 HU and 350 HU. The product of the
2 squares aims at pushing the estimates of u towards one of the values. The value u
0
corresponds to primarily the bone region where the updates for u should keep its high
attenuationcoefficient,where0.096isencouragingtheconvergencetowardsthenormal
softtissueattenuationcoefficient.
5.3 Results
5.3.1 2DhistogramofnoncontrastedandcontrastCT
WefirststudythechangeofCTHounsfieldunitafterconsumingcontrastenhanceagent
with the 4-phase CT data. The 4-phase CT data refers to the contrast enhanced CT
collected at 4 different time points: pre-phase, arterial phase, venous phase, delayed
phase. Weselectthepre-phaseandthevenousphase(usuallycorrespondstothecontrast
enhance CT in clinic) for comparison. Figure 5.2 shows the 2D histogram of the pre-
phase and venous phase CT image after image registration. Clearly there are 2 peaks,
correspondingtotheenhancementofCTvaluewiththeadditionalinfluenceofcontrast
agent. Note here the CT value is normalized by adding 1000 so that the intensity of air
have 0 HU. By carefully studying the 2D histogram, we determine that the soft tissue
regions that suffer from the influence of the contrast agent will have a CT Hounsfield
unitintherange[70350]HU.
113
1200
1400
1600
1000
1200
1400
1600
0
5
10
15
x 10
4

venous phase CT value (HU)
pre phase CT value (HU)
0
5
10
15
x 10
4
Figure 5.2: the 2D histogram of pre-phase with no contrast agent and the venous phase
afterthepatientconsumescontrastagent. The2peakscorrespondtotheincreaseofCT
valuesinsofttissuessuchasliverandheart.
5.3.2 3Dsimulationforestimationofemissionandattenuationmap
We simulate a small scale Siemens mCT system. The number of views is 192 and the
number of rays per angle is 400. To make computationally efficient simulation, the
number of rings is decreased to 9 and the maximum ring difference is 8. The image we
simulatehasdimension256x256x19.
We previous have reconstructed a PET image from clinical dataset with no contrast
enhanced CT with MAP [9]. This reconstructed emission image could be regarded
as ground truth for the emission activity. The reconstructed emission image and the
attenuation umap are shown in Fig. 5.3, after truncating the axial direction so that the
image’sdimensionis256x256x19. Thistruncatedemission-attenuationpairisthenused
inoursimulation.
We then artificially change the umap value in the blood-rich regions such as heart
andliverwithreferencetotheaveragevalueafterapplyingtheconventionalCT-to-umap
114
Figure5.3: Thereconsructedemissionimagewithcorrectattenuationmap(top)andthe
originalnoncontrastedattenuationmap(bottom).
mapping to the contrast CT. By doing this we generate a new umap map with consider-
able error in blood-rich area such as heart and liver. This is a fairly good simulation of
the umap image we will get after performing the CT-to-umap mapping for the noncon-
trastedCTtothecontrastCTimage. Ofcourse,ifweusethisartificialcontrast-affected
attenuationmapinperformingimagereconstruction,wewouldexpecttoseelargeerror
in the final PET image. In the following sections, we will show this biased result in
comparisonwith ourjointestimationresults, which demonstratesthat thejointiterative
115
relativeerror emission attenuation
simple(spine) 11% 6.1%
simple(heart) 0.7% 0.5%
joint(spine) 1.9% 0.9%
joint(heart) 0.4% 0.2%
Table5.1: Comparisonoftherelativeerrorofthesimplerescalingmethodandourjoint
estimation method in estimating the attenuation coefficient and the emission activity.
The numbers are relative difference of the mean of ROIs for estimated values and the
groundtruth.
updates could accurately estimate emission and attenuation simultaneously. To speed
up convergence, we run a fast OSEM algorithm with this artificial attenuation map to
obtain a PET image as initialization emission data. The artificial attenuation map is
the initialization attenuation map. Our method takes these two as input and iteratively
updatesbothwithpriorinformationfromthecontrastenhancedCTimage.
Figure 5.4 shows the estimated emission and attenuation image computed from
simplerescalingschemeandourjointestimationalgorithm. ROIsfortheheart,spineis
chosen and the relative error with respect to the true value and the result is summarized
as in Table 5.1. Clearly for both the emission and attenuation, estimation from our
methodyieldssmallererrorthanthesimplerescalingmethod.
5.4 Discussion
Our study currently depends on a fairly good initialization of the umap and emission
image. Meanwhile,thepenaltytermwedesignallowsthecostfunctiontoapproachthe
global maxima in its neighborhood. It is possible that if the initialization of attenuation
umapispoor,i.e. farfromthetrueattenuationmap,thepriorimposedontheattenuation
mapmayleadthealgorithmtoconvergetothewrongpeak,classifyingboneattenuation
116
Figure5.4: Thereconstructedemissionandattenuationimagewiththesimplerescaling
method(top)andtheoneswithourjointestimationalgorithm(bottom).
coefficient as normal tissue or vice versa. Therefore a good initialization of the atten-
uation map and a preliminary reconstructed emission image is necessary to make the
algorithmworkwell.
As shown the prior information can lead the algorithm to converge towards a useful
localsolution. Howeveritisalsotheimpactfactorthatcaninfluencethefinalestimation.
Inourmodelweassumethattheundeterminedattenuationcoefficientsareconcentrated
near2peaks,correspondingtothenormaltissueandbone. Howevergenerallyspeaking
117
the algorithm will reconstruct more accurate attenuation image by having a more com-
plex model for the distribution of the attenuation coefficients for the contrast affected
regions.
A finer joint estimation algorithm may incorporate the segmentation of different
organs and tissues in the range of the ”fuzzy” CT values. In each update for the atten-
uation map, different regions will be treated differently. For example, the region in the
spine will keep the attenuation value unchanged while the regions in the heart or liver
willallowattenuationvaluemovetowardsthenormalvalueforsofttissues.
The application of this work can be extended to PETMR, where the MR field of
viewisusuallysmallerthanthatofthePETandcannotprovideattenuationinformation
of the limbs. The MR based attenuation estimation will provide accurate quantitative
estimation for the in-field of view region and the joint estimation will provide the extra
informationfortheoutside-fieldofviewregions.
5.5 Conclusion
We adopted the joint estimation method for emission and attenuation from [83] and
further modify it with the MAP-derivative for updating the emission source. The prior
from contrast enhanced CT is designed to encourage convergence towards 2 peaks cor-
responding to normal soft tissue and bone respectively. The result shows that the joint
estimation can provide more accurate estimates of both emission and attenuation map
thanthesimplerescalingmethod.
118
Chapter6
ConclusionandFuturework
6.1 Conclusion
We developed a wholebody Patlak parametric image estimation method for
18
F-FDG
PETscans. ThedynamicdataismodeledtofollowaninhomogeneousPoissonprocess.
We use a single data/single voxel Cramer-Rao study over a wide range of parameter
valuesandfoundthatPatlakslopeestimatesfromdualtime-pointdataconsistentlyout-
perform fractional SUV measures computed from the same data. We also derived the
formula for calculating the spatially variant penalty to ensure that reconstructed Patlak
imageshavea resolution thatdoes not depend on activityor count rate and furthermore
that the resolution of the slope and intercept parameters are matched. Through Monte
Carlo simulation over a range of lesion sizes, as with the Cramer Rao analysis, we find
that Patlak slope parameters are more accurate in differentiating tumor and background
thanisfractionSUV.Thisisparticularlythecaseforthesmallerlesionswherethemuch
smaller relative variance of Patlak measures lead to better differentiation of tumor and
background. Our listmode based estimation method also provides a more accurate esti-
mationmethodthanthesinogrambasedmethod,becauseofthebetterusageoftemporal
informationinthelistmodedata.
We then applied this wholebody Patlak estimation algorithm to clinical PET data
and developed practical solutions for several challenges, including estimation of blood
input function, modeling of random and scattered events, as well as motion correction.
119
We present preliminary clinical results comparing Patlak images with the SUV mea-
surements. Input function is estimated using a hybrid method incorporating population
informationandimagingdata. Dynamicmodelsaredevelopedforrandomandscattered
events. A non-regid registration method is used to compensate the motion between two
frames. ThewholebodypatientstudiesshowthatwecanproducePatlakimagesofsim-
ilarqualitytostaticPETscans. Thisisincontrasttodual-timepointfractionalSUVval-
ueswhichdonotproducegoodqualityimages. ThedataanalysissuggeststhatPatlakas
a quantitativemethod yields lessvariation among population than the semi-quantitative
methods such as SUV or %DSUV. The contrast recovery results also indicates that Pat-
lak is more suitable in the task of lesion detection. To conclude, our wholebody Patlak
estimation method is feasible for clinical
18
F-FDG PET studies, providing more stable
and reliable means for lesion detection, compared to the semiquantitative SUV mea-
surements.
As an extension of our estimation method to the reversible tracer, we developed a
modifiedalgorithmforestimatingLoganparametersfromlistmodedata. Thetheoryand
quantitativesimulationresultsstronglysupportthatourmethodcanachievesatisfactory
performance, while the direct usage of the original Patlak framework for data after the
steadystateyieldshugevariancebecauseoftheill-conditioning. Ourmethodalsoshow
its advantage over the sinogram based approach in that it uses an accurate statistical
noise model, with Monte-Carlo simulation demonstrating that the variance of Logan
estimators in our method is uniformly lower than the ones provided by the sinogram
basedapproach.
Finally we studied the influence of the contrast enhanced CT to the quantitative
analysisinPET.SimplyapplyingtheCT-to-umapmappingtothecontrastCTwillincur
large error to the attenuation thus influence the estimation of emission too. We used
MLAA (maximum likelihood estimation of activity and attenuation) to perform a joint
120
estimationofactivityandattenuationmapusingpriorinformation. Itisshownthatwith
the prior and a reasonable model for the distribution of attenuation coefficients. The
estimation of attenuation and emission images can be estimated for practical use and
hasmuchsmallestimationerrorthanthesimplescalingscheme.
6.2 Futurework
Extending the ROI study to multiple organs and tissues in the wholebody Patlak image
will be one of our primary interests in our future research. There are several issues that
affect the performance of our estimation method which need to be further investigated.
Thefirstchallengeistoenhancetheconvergenceinoptimizingthecostfunction. Right
now we are using 4D incremental gradient method in optimization. However, the step-
sizeinourframeworkisempirical. Asweknowforgradientbasedmethod,thestepsize
canensureconvergenceifitsatisfiesthefollowingcondition:
1
X
i=1
a
i
>1
1
X
i=1
a
i
2
<1 (6.1)
Rightnowweuseasimplestepsizechoice
a
i
=
1
i=30+1:5
(6.2)
However, this might not be optimized as we have seen in phantom simulation that
the convergence in some regions will be quite slow (usually 100-200 iterations). For
example for the region that is assigned the activity curve identical to the blood input
function, we should expect to have estimated Patlak slope and intercept equal to 0 and
1 respectively. However, we did not observe convergence until 150 iterations. This
might be impractical to real data, as it will take more than several hours to obtain the
121
results. Considering different data size and the complexity of real data, it may take
evenlongertoachieveconvergence. Thereforeoptimizingthestepsizeisveryimportant
to achieve fast convergence rate. Intuitively we can enhance the convergence speed
by increasing the stepsize, yet this may lead to divergence of the algorithm. What we
recent implemented is an Armijo-based stepsize choosing technique. We initialize the
stepsizeateachiterationtoberelativelybig. Thecostfunctionaftereachupdatewillbe
comparedtothepreviouscost. Ifthecurrentcostishigher,thestepsizewillbedecreased
by half. This comparison of cost function will continue until the new cost function is
smaller than the previous one. This method yields faster convergence as well as ensure
convergence. Another method is to analyze the convergence rate for different regions
andoptimizethestepsizeforeachregionindividually.
The second issue is to optimize the wholebody dynamic PET scan protocol for rou-
tine clinical studies. Right now we start our first frame acquisition around 1 hour after
injection. WethenstarttheclinicalscanandobtainCTandwholebodystaticPETimage
after the first frame. After the clinical scan we start the second frame acquisition. The
patientstaysinthePETscannerforalongtime(about90minutes). Wewillinvestigate
whether decreasing the scanning time for each bed position and reduce the wholebody
scanto4bedpositionswillretaintheimagequalityandstatisticalaccuracy. Wewillalso
investigate whether decreasing the waiting time between 2 listmode frames or shorten
the time from injection to the first frames starting time will impair the performance of
ourmethod. Shortentheframedurationandtimebetweentwoframescanreducemotion
artifactsthusimproveimagequality.
The scan duration issue is even more important for the Logan estimation. In logan
analysis, data startingfrom the injectiontime till the steady state is required to improve
the ill-conditioning of the estimation problem. What is really desirable is an algorithm
with comparible conditioning to our method yet only requires short collection time
122
after the steady state. This is difficult based on the current Logan model. In future
we will investigate in the 2 tissue compartment model hopefully in a totally different
way and come up with a new physiologically meaningful biomarker that allows well-
conditioninglistmodebasedestimationwithonlypartiallistmodedata.
The joint estimation of emission activity and the attenuation is very challenging for
nonTOF PET. The study of [81] has shown that the Fisher information matrix of the
joint unknowns is ill-conditioned. As we have shown, the contrast enhanced CT can be
used as useful prior to the joint estimation algorithm. However, it is difficult to adjust
the weights of the CT intensity prior. A possible improvement to the current method
is to combine the joint estimation with segmentation. The segmentation will add more
confidencetotheunknownattenuationofeachROIandcanavoidthealgorithmsticking
tothelocalsub-optimalsolutions.
The final goal and most important application of this wholebody approach is to
improve sensitivity and specificity in clinical
18
F-FDG studies. In this thesis we have
performedquantitativeanalysisonasmallnumberofpatientdatasetsandcomputedthe
absolute value of Patlak parameters and SUV values as well as tumor-to-tissue contrast
foronepatient. PreliminaryresultsindicatethatPatlakslopemayhavebettersensitivity
and specificity in oncologic applications such as staging and evaluation of response to
therapy. Therefore we aim at interpreting the results of future cancer patients and com-
pare for example the contrast of tumor against background for the Patlak slope image
andSUVimages. Wealsoplantostudythebehaviorofthe2methodsinregionscorre-
spondingtoinflammationaswell. Weknowthatinflammationalsoshowhighintensity
against normal tissue background, however the activity will gradually decrease over
time. Therefore it is possible that with Patlak method, we will not observer quite high
intensity while with SUV method we will see the high intensity cause by inflamma-
tion. We are interested in applying the method to clinical data with reversible tracer
123
as well, aiming at improving the quantitative results with our direct Logan parameter
estimation.
124
ReferenceList
[1] L.K.Shankar,J.M.Hoffman,S.Bacharach,M.M.Graham,J.Karp,A.A.Lam-
mertsma, S. Larson, D. A. Mankoff, B. A. Siegel, A. Van den Abbeele, et al.,
“Consensus recommendations for the use of 18f-fdg pet as an indicator of thera-
peutic response in patients in national cancer institute trials,” Journal of Nuclear
Medicine,vol.47,no.6,pp.1059–1066,2006.
[2] E. M. Rohren, T. G. Turkington, and R. E. Coleman, “Clinical Applications of
PETinOncology,” Radiology,vol.231,no.2,pp.305–332,2004.
[3] G. El Fakhri, S. Surti, C. M. Trott, J. Scheuermann, and J. S. Karp, “Improve-
ment in lesion detection with whole-body oncologic time-of-flight pet,” Journal
of Nuclear Medicine,vol.52,no.3,pp.347–353,2011.
[4] L. A. Shepp and B. F. Logan, “The fourier reconstruction of a head section,”
Nuclear Science, IEEE Transactions on,vol.21,no.3,pp.21–43,1974.
[5] J. G. Colsher, “Fully-three-dimensional positron emission tomography,” Physics
in medicine and biology,vol.25,no.1,p.103,1980.
[6] M.E.Daube-WitherspoonandG.Muehllehner,“Treatmentofaxialdatainthree-
dimensionalpet,” J Nucl Med,vol.28,no.11,pp.1717–1724,1987.
[7] R. M. Leahy and J. Qi, “Statistical approaches in quantitative positron emission
tomography,” Statistics and Computing,vol.10,no.2,pp.147–165,2000.
[8] M.YavuzandJ.A.Fessler,“Maximumlikelihoodemissionimagereconstruction
forrandoms-precorrectedPETscans,”inNuclearScienceSymposiumConference
Record, 2000 IEEE,vol.2,pp.15–229,IEEE,2000.
[9] J. Qi, R. M. Leahy, S. R. Cherry, A. Chatziioannou, and T. H. Farquhar, “High-
resolution 3d bayesian image reconstruction using the micropet small-animal
scanner,” Physics in medicine and biology,vol.43,no.4,p.1001,1998.
125
[10] H.M.HudsonandR.S.Larkin,“Acceleratedimagereconstructionusingordered
subsets of projection data,” Medical Imaging, IEEE Transactions on, vol. 13,
no.4,pp.601–609,1994.
[11] S.AhnandJ.A.Fessler,“Globallyconvergentimagereconstructionforemission
tomography using relaxed ordered subsets algorithms,” Medical Imaging, IEEE
Transactions on,vol.22,no.5,pp.613–626,2003.
[12] L. G. Strauss and P. S. Conti, “The applications of pet in clinical oncology.,”
Journal of nuclear medicine: official publication, Society of Nuclear Medicine,
vol.32,no.4,pp.623–48,1991.
[13] R. L. Wahl, H. Jacene, Y. Kasamon, and M. A. Lodge, “From recist to per-
cist: evolving considerations for pet response criteria in solid tumors,” Journal
of nuclear medicine,vol.50,no.Suppl1,pp.122S–150S,2009.
[14] R. Boellaard, “Need for standardization of 18f-fdg pet/ct for treatment response
assessments,” Journal of Nuclear Medicine, vol. 52, no. Supplement 2, pp. 93S–
100S,2011.
[15] P.Masa-Ah,M.Tuntawiroon,andS.Soongsathitanon,“Anovelschemeforstan-
dardized uptake value (suv) calculation in pet scans,” International Journal of
Mathematical Models and Methods in Applied Sciences, Naun, vol. 4, pp. 291–
299,2010.
[16] K. Wiyaporn, C. Tocharoenchai, P. Pusuwan, T. Ekjeen, S. Leaungwutiwong,
and S. Thanyarak, “Factors affecting standardized uptake value (suv) of positron
emission tomography (pet) imaging with 18f-fdg,” J Med Assoc Thai, vol. 93,
no.1,pp.108–114,2010.
[17] C.NahmiasandL.M.Wahl,“Reproducibilityofstandardizeduptakevaluemea-
surements determined by 18f-fdg pet in malignant tumors,” Journal of Nuclear
Medicine,vol.49,no.11,pp.1804–1808,2008.
[18] R. Boellaard, “Optimisation and harmonisation: two sides of the same coin?,”
European journal of nuclear medicine and molecular imaging, vol. 40, no. 7,
pp.982–984,2013.
[19] M. A. Lodge, M. A. Chaudhry, and R. L. Wahl, “Noise considerations for pet
quantification using maximum and peak standardized uptake value,” Journal of
Nuclear Medicine,vol.53,no.7,pp.1041–1047,2012.
[20] K. Alkhawaldeh, G. Bural, R. Kumar, and A. Alavi, “Impact of dual-time-point
18f-fdg pet imaging and partial volume correction in the assessment of solitary
pulmonarynodules,”Europeanjournalofnuclearmedicineandmolecularimag-
ing,vol.35,no.2,pp.246–252,2008.
126
[21] E. Prieto, J. M. Mart´ ı-Climent, I. Dom´ ınguez-Prado, P. Garrastachu, R. D´ ıez-
Valle, S. Tejada, J. J. Aristu, I. Pe˜ nuelas, and J. Arbizu, “Voxel-based analysis
of dual-time-point 18f-fdg pet images for brain tumor identification and delin-
eation,” Journal of Nuclear Medicine,vol.52,no.6,pp.865–872,2011.
[22] H. Zhuang, M. Pourdehnad, E. S. Lambright, A. J. Yamamoto, M. Lanuti, P. Li,
P. D. Mozley, M. D. Rossman, S. M. Albelda, and A. Alavi, “Dual time point
18f-fdg pet imaging for differentiating malignant from inflammatory processes,”
Journal of Nuclear Medicine,vol.42,no.9,pp.1412–1417,2001.
[23] R. Kumar, V. A. Loving, A. Chauhan, H. Zhuang, S. Mitchell, and A. Alavi,
“Potential of dual-time-point imaging to improve breast cancer diagnosis with
18f-fdgpet,”JournalofNuclearMedicine,vol.46,no.11,pp.1819–1824,2005.
[24] F.J.Cloran,K.P.Banks,W.S.Song,Y.Kim,andY.C.Bradley,“Limitationsof
dualtimepointpetintheassessmentoflungnoduleswithlowfdgavidity,”Lung
Cancer,vol.68,no.1,pp.66–71,2010.
[25] W.-L. Chan, S. C. Ramsay, E. R. Szeto, J. Freund, J. M. Pohlen, L. C. Tarlinton,
A. Young, A. Hickey, and R. Dura, “Dual-time-point 18f-fdg-pet/ct imaging in
theassessmentofsuspectedmalignancy,”Journalofmedicalimagingandradia-
tion oncology,vol.55,no.4,pp.379–390,2011.
[26] G. Antoch, F. M. Vogt, L. S. Freudenberg, F. Nazaradeh, S. C. Goehde,
J. Barkhausen, G. Dahmen, A. Bockisch, J. F. Debatin, and S. G. Ruehm,
“Whole-body dual-modality PET/CT and whole-body MRI for tumor staging in
oncology,” Jama,vol.290,no.24,pp.3199–3206,2003.
[27] C. K. Hoh, R. A. Hawkins, J. A. Glaspy, M. Dahlbom, Y. T. Nielson, E. J. Hoff-
man, C. Schiepers, Y. Choi, S. Rege, E. Nitzsche, et al., “Cancer detection with
whole-body pet using 2-[18f] fluoro-2-deoxy-d-glucose,” Journal of computer
assisted tomography,vol.17,no.4,pp.582–589,1993.
[28] G. Jerusalem, Y. Beguin, M.-F. Fassotte, F. Najjar, P. Paulus, P. Rigo, and G. Fil-
let, “Whole-body positron emission tomography using18f-fluorodeoxyglucose
for posttreatment evaluation in hodgkins disease and non-hodgkins lymphoma
has higher diagnostic and prognostic value than classical computed tomography
scanimaging,” Blood,vol.94,no.2,pp.429–433,1999.
[29] P.E.Valk, Positron emission tomography: basic sciences. Springer,2003.
[30] K.R.ZasadnyandR.Wahl,“Standardizeduptakevaluesofnormaltissuesatpet
with 2-[fluorine-18]-fluoro-2-deoxy-d-glucose: variations with body weight and
amethodforcorrection.,” Radiology,vol.189,no.3,pp.847–850,1993.
127
[31] M.Graham,L.Peterson,andR.Hayward,“Comparisonofsimplifiedquantitative
analyses of fdg uptake,” Nuclear medicine and biology, vol. 27, no. 7, pp. 647–
655,2000.
[32] M. D. Walker, M. Asselin, P. J. Julyan, M. Feldmann, P. Talbot, T. Jones, and
J. Matthews, “Bias in iterative reconstruction of low-statistics pet data: benefits
of a resolution model,” Physics in medicine and biology, vol. 56, no. 4, p. 931,
2011.
[33] J.QiandR.M.Leahy,“Resolutionandnoisepropertiesofmapreconstructionfor
fully 3-d pet,” Medical Imaging, IEEE Transactions on, vol. 19, no. 5, pp. 493–
506,2000.
[34] Q.Li,B.Bai,S.Cho,A.Smith,andR.M.Leahy,“Countindependentresolution
anditscalibration,” Proc. 10th Fully3D Image Recon,2009.
[35] N.A.Karakatsanis,M.A.Lodge,A.K.Tahari,Y.Zhou,R.L.Wahl,andA.Rah-
mim, “Dynamic whole-body PET parametric imaging: I. Concept, acquisition
protocol optimization and clinical application,” Phys Med Biol, vol. 58, no. 20,
p.7391,2013.
[36] N. A. Karakatsanis, M. A. Lodge, Y. Zhou, R. L. Wahl, and A. Rahmim,
“Dynamic whole-body pet parametric imaging: Ii. task-oriented statistical esti-
mation,” Physics in medicine and biology,vol.58,no.20,p.7419,2013.
[37] R.E.Carson,“Tracerkineticmodelinginpet,”inPositronEmissionTomography,
pp.127–159,Springer,2005.
[38] M. Takesh, “The potential benefit by application of kinetic analysis of pet in the
clinicaloncology,” ISRN oncology,vol.2012,2012.
[39] Y. Sugawara, K. R. Zasadny, H. B. Grossman, I. R. Francis, M. F. Clarke, and
R. L. Wahl, “Germ cell tumor: Differentiation of viable tumor, mature teratoma,
and necrotic tissue with fdg pet and kinetic modeling 1,” Radiology, vol. 211,
no.1,pp.249–256,1999.
[40] T. E. Nichols, J. Qi, E. Asma, and R. M. Leahy, “Spatiotemporal reconstruction
of list-mode pet data,” Medical Imaging, IEEE Transactions on, vol. 21, no. 4,
pp.396–404,2002.
[41] E. Asma, D. W. Shattuck, and R. M. Leahy, “Lossless compression of dynamic
petdata,”NuclearScience,IEEETransactionson,vol.50,no.1,pp.9–16,2003.
[42] W.Zhu,Q.Li,B.Bai,P.Conti,andR.Leahy,“Patlakimageestimationfromdual
time-pointlist-modePETData,”vol.33,no.4,pp.913–924,2014.
128
[43] D.L.Snyder,“Parameterestimationfordynamicstudiesinemission-tomography
systemshavinglist-modedata,” Nuclear Science, IEEE Transactions on,vol.31,
no.2,pp.925–931,1984.
[44] C.S.PatlakandR.G.Blasberg,“Graphical evaluationof blood-to-brain transfer
constants from multiple-time uptake data. generalizations,” J Cereb Blood flow
metab,vol.5,no.4,pp.584–590,1985.
[45] C.S.Patlak,R.G.Blasberg,J.D.Fenstermacher,etal.,“Graphicalevaluationof
blood-to-braintransferconstantsfrommultiple-timeuptakedata,”JCerebBlood
Flow Metab,vol.3,no.1,pp.1–7,1983.
[46] F. O’Sullivan, J. Kirrane, M. Muzi, J. N. O’Sullivan, A. M. Spence, D. A.
Mankoff,andK.A.Krohn,“Kineticquantitationofcerebralpet-fdgstudieswith-
outconcurrentbloodsampling: statisticalrecoveryofthearterialinputfunction,”
Medical Imaging, IEEE Transactions on,vol.29,no.3,pp.610–624,2010.
[47] Z. J. Wang, P. Qiu, K. R. Liu, and Z. Szabo, “Model-based receptor quantization
analysis for pet parametric imaging,” in Engineering in Medicine and Biology
Society, 2005. IEEE-EMBS 2005. 27th Annual International Conference of the,
pp.5908–5911,IEEE,2006.
[48] M.Naganawa,Y.Kimura,J.Yano,M.Mishina,M.Yanagisawa,K.Ishii,K.Oda,
and K. Ishiwata, “Robust estimation of the arterial input function for logan plots
using an intersectional searching algorithm and clustering in positron emission
tomography for neuroreceptor imaging,” Neuroimage, vol. 40, no. 1, pp. 26–34,
2008.
[49] J. Logan, “Graphical analysis of pet data applied to reversible and irreversible
tracers,” Nuclear medicine and biology,vol.27,no.7,pp.661–670,2000.
[50] R. N. Gunn, S. R. Gunn, F. E. Turkheimer, J. A. Aston, and V. J. Cunningham,
“Positron emission tomography compartmental models: a basis pursuit
strategy for kinetic modeling,” Journal of Cerebral Blood Flow & Metabolism,
vol.22,no.12,pp.1425–1439,2002.
[51] M. Slifstein and M. Laruelle, “Effects of statistical noise on graphic analysis
of pet neuroreceptor studies,” Journal of Nuclear Medicine, vol. 41, no. 12,
pp.2083–2088,2000.
[52] D.Feng,Z.Wang,andS.-C.Huang,“Astudyonstatisticallyreliableandcompu-
tationallyefficientalgorithmsforgeneratinglocalcerebralbloodflowparametric
images with positron emission tomography,” Medical Imaging, IEEE Transac-
tions on,vol.12,no.2,pp.182–188,1993.
129
[53] D. Feng, S.-C. Huang, Z. Wang, and D. Ho, “An unbiased parametric imaging
algorithm for nonuniformly sampled biomedical system parameter estimation,”
Medical Imaging, IEEE Transactions on,vol.15,no.4,pp.512–518,1996.
[54] R.E.CarsonandK.Lange,“Comment: Theemparametricimagereconstruction
algorithm,” Journal of the American Statistical Association, vol. 80, no. 389,
pp.20–22,1985.
[55] A. J. Reader, F. C. Sureau, C. Comtat, R. Tr´ ebossen, and I. Buvat, “Joint estima-
tion of dynamic pet images and temporal basis functions using fully 4d ml-em,”
Physics in Medicine and Biology,vol.51,no.21,p.5455,2006.
[56] J. Matthews, D. Bailey, P. Price, and V. Cunningham, “The direct calculation of
parametric images from dynamic pet data using maximum-likelihood iterative
reconstruction,” Physics in medicine and biology,vol.42,no.6,p.1155,1997.
[57] M. E. Kamasak, C. A. Bouman, E. D. Morris, and K. Sauer, “Direct reconstruc-
tionofkineticparameterimagesfromdynamicpetdata,”MedicalImaging,IEEE
Transactions on,vol.24,no.5,pp.636–650,2005.
[58] G. Wang and J. Qi, “Generalized algorithms for direct reconstruction of para-
metric images from dynamic pet data,” Medical Imaging, IEEE Transactions on,
vol.28,no.11,pp.1717–1726,2009.
[59] G. Wang, L. Fu, and J. Qi, “Maximum a posteriori reconstruction of the patlak
parametricimagefromsinogramsindynamicpet,”Physicsinmedicineandbiol-
ogy,vol.53,no.3,p.593,2008.
[60] N.A.Karakatsanis,M.A.Lodge,A.K.Tahari,Y.Zhou,R.L.Wahl,andA.Rah-
mim,“Dynamicwhole-bodypetparametricimaging: I.concept,acquisitionpro-
tocol optimization and clinical application,” Physics in medicine and biology,
vol.58,no.20,p.7391,2013.
[61] W. Zhu, Q. Li, and R. M. Leahy, “Dual-time-point patlak estimation from list
modepetdata,”inBiomedicalImaging(ISBI),20129thIEEEInternationalSym-
posium on,pp.486–489,IEEE,2012.
[62] Q. Li, E. Asma, S. Ahn, and R. M. Leahy, “A fast fully 4-d incremental gradient
reconstruction algorithm for list mode pet data,” Medical Imaging, IEEE Trans-
actions on,vol.26,no.1,pp.58–67,2007.
[63] Q. Li and R. M. Leahy, “Direct estimation of patlak parameters from list mode
pet data,” in Biomedical Imaging: From Nano to Macro, 2009. ISBI’09. IEEE
International Symposium on,pp.390–393,IEEE,2009.
130
[64] Y. Zhou, J. Tang, and D. Wong, “Direct 4d parametric image reconstruction
with plasma input and reference tissue models in reversible binding imaging,”
in Nuclear Science Symposium Conference Record (NSS/MIC), 2009 IEEE,
pp.2516–2522,IEEE,2009.
[65] S. L. Bacharach, “Pet/ct attenuation correction: breathing lessons,” Journal of
Nuclear Medicine,vol.48,no.5,pp.677–679,2007.
[66] P.E.Kinahan,B.H.Hasegawa,andT.Beyer,“X-ray-basedattenuationcorrection
for positron emission tomography/computed tomography scanners,” in Seminars
in nuclear medicine,vol.33,pp.166–179,Elsevier,2003.
[67] C. Burger, G. Goerres, S. Schoenes, A. Buck, A. Lonn, and G. Von Schulthess,
“Petattenuationcoefficientsfromctimages: experimentalevaluationofthetrans-
formation of ct into pet 511-kev attenuation coefficients,” European journal of
nuclear medicine and molecular imaging,vol.29,no.7,pp.922–927,2002.
[68] P. Kinahan, D. Townsend, T. Beyer, and D. Sashin, “Attenuation correction for
a combined 3d pet/ct scanner,” Medical physics, vol. 25, no. 10, pp. 2046–2053,
1998.
[69] G. Antoch, L. S. Freudenberg, T. Beyer, A. Bockisch, and J. F. Debatin, “To
enhance or not to enhanced 18f-fdg and ct contrast agents in dual-modality 18f-
fdgpet/ct,”JournalofNuclearMedicine,vol.45,no.1suppl,pp.56S–65S,2004.
[70] J. Fletcher, C. D. Johnson, W. R. Krueger, D. A. Ahlquist, H. Nelson, D. Ilstrup,
W. S. Harmsen, and K. E. Corcoran, “Contrast-enhanced ct colonography
in recurrent colorectal carcinoma: feasibility of simultaneous evaluation for
metastatic disease, local recurrence, and metachronous neoplasia in colorectal
carcinoma,” American Journal of Roentgenology, vol. 178, no. 2, pp. 283–290,
2002.
[71] J.Oliver3rd,R.Baron,M.Federle,andH.RocketteJr,“Detectinghepatocellular
carcinoma: valueofunenhancedorarterialphasectimagingorbothusedincon-
junction with conventional portal venous phase contrast-enhanced ct imaging.,”
AJR. American journal of roentgenology,vol.167,no.1,pp.71–77,1996.
[72] N. G. Schaefer, T. F. Hany, C. Taverna, B. Seifert, K. D. Stumpe, G. K. von
Schulthess, and G. W. Goerres, “Non-hodgkin lymphoma and hodgkin disease:
coregisteredfdgpetandctatstagingandrestaging-doweneedcontrast-enhanced
ct?,” RADIOLOGY-OAK BROOK IL-,vol.232,no.3,pp.823–829,2004.
[73] Y. Nakamoto, B. B. Chin, D. L. Kraitchman, L. P. Lawler, L. T. Marshall, and
R. L. Wahl, “Effects of nonionic intravenous contrast agents at pet/ct imaging:
Phantomandcaninestudies1,” Radiology,vol.227,no.3,pp.817–824,2003.
131
[74] S. A. Nehmeh, Y. E. Erdi, H. Kalaigian, K. S. Kolbert, T. Pan, H. Yeung,
O. Squire, A. Sinha, S. M. Larson, and J. L. Humm, “Correction for oral con-
trastartifactsinctattenuation-correctedpetimagesobtainedbycombinedpet/ct,”
Journal of Nuclear Medicine,vol.44,no.12,pp.1940–1944,2003.
[75] Y.-Y. Yau, W.-S. Chan, Y.-M. Tam, P. Vernon, S. Wong, M. Coel, and S. K.-
F. Chu, “Application of intravenous contrast in pet/ct: does it really introduce
significant attenuation correction error?,” Journal of Nuclear Medicine, vol. 46,
no.2,pp.283–291,2005.
[76] Y.Tai,K.Lin,M.Dahlbom,S.-C.Huang,andE.Hoffman,“Ahybridattenuation
correction technique to compensate for lung density in 3d total body pet,” in
NuclearScienceSymposiumandMedicalImagingConference,1994.,1994IEEE
Conference Record,vol.4,pp.1643–1647,IEEE.
[77] M. Xu, P. Cutler, and W. Luk, “Adaptive, segmented attenuation correction for
whole-bodypetimaging,”NuclearScience,IEEETransactionson,vol.43,no.1,
pp.331–336,1996.
[78] Y. Censor, D. E. Gustafson, A. Lent, and H. Tuy, “A new approach to the emis-
sioncomputerizedtomographyproblem: simultaneouscalculationofattenuation
and activity coefficients,” Nuclear Science, IEEE Transactions on, vol. 26, no. 2,
pp.2775–2779,1979.
[79] A.Welch,R.Clack,F.Natterer,andG.T.Gullberg,“Towardaccurateattenuation
correctioninspectwithouttransmissionmeasurements,”MedicalImaging,IEEE
Transactions on,vol.16,no.5,pp.532–541,1997.
[80] F.Natterer,“Determinationoftissueattenuationinemissiontomographyofopti-
callydensemedia,” Inverse Problems,vol.9,no.6,p.731,1993.
[81] A. Rezaei, J. Nuyts, M. Defrise, G. Bal, C. Michel, M. Conti, and C. Watson,
“Simultaneousreconstructionofactivityandattenuationintime-of-flightpet,”in
Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), 2011
IEEE,pp.2375–2382,IEEE,2011.
[82] A. Salomon, A. Goedicke, B. Schweizer, T. Aach, and V. Schulz, “Simultaneous
reconstruction of activity and attenuation for pet/mr,” Medical Imaging, IEEE
Transactions on,vol.30,no.3,pp.804–813,2011.
[83] J. Nuyts, P. Dupont, S. Stroobants, R. Benninck, L. Mortelmans, and P. Suetens,
“Simultaneous maximum a posteriori reconstruction of attenuation and activity
distributionsfromemissionsinograms,”MedicalImaging,IEEETransactionson,
vol.18,no.5,pp.393–403,1999.
132
[84] M.Defrise,A.Rezaei,andJ.Nuyts,“Time-of-flightpetdatadeterminetheatten-
uationsinogramuptoaconstant,”Physicsinmedicineandbiology,vol.57,no.4,
p.885,2012.
[85] R. Boellaard, M. Hofman, O. Hoekstra, and A. Lammertsma, “Accurate pet/mr
quantification using time of flight mlaa image reconstruction,” Molecular Imag-
ing and Biology,pp.1–9,2014.
[86] J.HamillandV.Panin,“Tof-mlaaforattenuationcorrectioninthoracicpet/ct,”in
Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), 2012
IEEE,pp.4040–4047,IEEE,2012.
[87] P. Iozzo, F. Geisler, V. Oikonen, M. M¨ aki, T. Takala, O. Solin, E. Ferrannini,
J. Knuuti, and P. Nuutila, “Insulin stimulates liver glucose uptake in humans: an
18F-FDGPETStudy,” J Nucl Med,vol.44,no.5,pp.682–689,2003.
[88] M.Hu,A.Han,L.Xing,W.Yang,Z.Fu,C.Huang,P.Zhang,L.Kong,andJ.Yu,
“Value of dual-time-point fdg pet/ct for mediastinal nodal staging in non-small-
cell lung cancer patients with lung comorbidity,” Clinical nuclear medicine,
vol.36,no.6,pp.429–433,2011.
[89] E. Asma and R. M. Leahy, “Temporally invariant uniform spatial resolution in
dynamic pet,” in Nuclear Science Symposium Conference Record, 2003 IEEE,
vol.5,pp.3092–3096,IEEE,2003.
[90] J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized-
likelihoodimagereconstruction: space-invarianttomographs,”ImageProcessing,
IEEE Transactions on,vol.5,no.9,pp.1346–1358,1996.
[91] E. Asma and R. M. Leahy, “Mean and covariance properties of dynamic pet
reconstructions from list-mode data,” Medical Imaging, IEEE Transactions on,
vol.25,no.1,pp.42–54,2006.
[92] A. Matthies, M. Hickeson, A. Cuchiara, and A. Alavi, “Dual time point 18f-
fdg pet for the evaluation of pulmonary nodules,” Journal of Nuclear Medicine,
vol.43,no.7,pp.871–875,2002.
[93] H. Young, R. Baum, U. Cremerius, K. Herholz, O. Hoekstra, A. Lammertsma,
J.Pruim,andP.Price,“Measurementofclinicalandsubclinicaltumourresponse
using [
18
F]-fluorodeoxyglucose and positron emission tomography: review and
1999 EORTC recommendations,” Eur J Cancer, vol. 35, no. 13, pp. 1773–1782,
1999.
[94] D. Delbeke, R. E. Coleman, M. J. Guiberteau, M. L. Brown, H. D. Royal, B. A.
Siegel,D.W.Townsend,L.L.Berland,J.A.Parker,K.Hubner,etal.,“Procedure
133
guideline for tumor imaging with 18F-FDG PET/CT 1.0,” J Nucl Med, vol. 47,
no.5,pp.885–895,2006.
[95] S.-C. Huang, “Anatomy of SUV,” Nucl Med Biol, vol. 27, no. 7, pp. 643–646,
2000.
[96] R. Boellaard, “Need for standardization of 18F-FDG PET/CT for treatment
response assessments,” J Nucl Med, vol. 52, no. Supplement 2, pp. 93S–100S,
2011.
[97] R. Boellaard, “Standards for PET image acquisition and quantitative data analy-
sis,” J Nucl Med,vol.50,no.Suppl1,pp.11S–20S,2009.
[98] U. Metser and E. Even-Sapir, “Increased
18
F-Fluorodeoxyglucose Uptake in
Benign, Nonphysiologic Lesions Found on Whole-Body Positron Emission
Tomography/Computed Tomography (PET/CT): Accumulated Data From Four
Years of Experience With PET/CT,” in Semin Nucl Med, vol. 37, pp. 206–222,
Elsevier,2007.
[99] P. D. Shreve, Y. Anzai, and R. L. Wahl, “Pitfalls in oncologic diagnosis with
FDG PET imaging: physiologic and benign variants,” Radiographics, vol. 19,
no.1,pp.61–77,1999.
[100] L.G.Strauss,“Fluorine-18deoxyglucoseandfalse-positiveresults: amajorprob-
lem in the diagnostics of oncological patients,” Eur J Nucl Med, vol. 23, no. 10,
pp.1409–1415,1996.
[101] S. M. Larson, Y. Erdi, T. Akhurst, M. Mazumdar, H. A. Macapinlac, R. D. Finn,
C. Casilla, M. Fazzari, N. Srivastava, H. W. Yeung, et al., “Tumor treatment
response based on visual and quantitative changes in global tumor glycolysis
usingPET-FDGimaging: thevisualresponsescoreandthechangeintotallesion
glycolysis,” Clin Positron Imaging,vol.2,no.3,pp.159–171,1999.
[102] M. Imbriaco, M. G. Caprio, G. Limite, L. Pace, T. De Falco, E. Capuano, and
M. Salvatore, “Dual-time-point 18f-fdg pet/ct versus dynamic breast mri of sus-
picious breast lesions,” American Journal of Roentgenology, vol. 191, no. 5,
pp.1323–1330,2008.
[103] H.Wu,A.Dimitrakopoulou-Strauss,T.O.Heichel,B.Lehner,L.Bernd,V.Ewer-
beck, C. Burger, and L. G. Strauss, “Quantitative evaluation of skeletal tumours
withdynamicFDGPET:SUVincomparisontoPatlakanalysis,”EurJNuclMed,
vol.28,no.6,pp.704–710,2001.
134
[104] N. M. Freedman, S. K. Sundaram, K. Kurdziel, J. A. Carrasquillo, M. Whatley,
J.M.Carson,D.Sellers,S.K.Libutti,J.C.Yang,andS.L.Bacharach,“Compar-
ison of SUV and Patlak slope for monitoring of cancer therapy using serial PET
scans,” Eur J Nucl Med and molecular imaging,vol.30,no.1,pp.46–53,2003.
[105] E. P. Visser, M. E. Philippens, L. Kienhorst, J. H. Kaanders, F. H. Corstens, L.-
F. de Geus-Oei, and W. J. Oyen, “Comparison of tumor volumes derived from
glucose metabolic rate maps and SUV maps in dynamic 18F-FDG PET,” J Nucl
Med,vol.49,no.6,pp.892–898,2008.
[106] A. Magri, A. Krol, W. Lee, E. Lipson, W. McGraw, and D. Feiglin, “A new
method to determine probability of malignancy using dynamic breast F-18-FDG
PETstudies,” J Nucl Med meeting abstract,vol.50,no.Suppl2,p.1445,2009.
[107] G. Wang, L. Fu, and J. Qi, “Maximum a posteriori reconstruction of the patlak
parametricimagefromsinogramsindynamicpet,”Physicsinmedicineandbiol-
ogy,vol.53,no.3,p.593,2008.
[108] G. Z. Ferl, X. Zhang, H.-M. Wu, and S.-C. Huang, “Estimation of the 18f-fdg
input function in mice by use of dynamic small-animal pet and minimal blood
sampledata,”JournalofNuclearMedicine,vol.48,no.12,pp.2037–2045,2007.
[109] K.-P. Wong, X. Zhang, and S.-C. Huang, “Improved derivation of input function
in dynamic mouse [18f] fdg pet using bladder radioactivity kinetics,” Molecular
Imaging and Biology,pp.1–11,2013.
[110] D.Kroon,“registrationtoolboxnon-rigidb-splinegridimageregistration,”2008.
[111] A. A. Joshi, D. W. Shattuck, P. M. Thompson, and R. M. Leahy, “Surface-
constrained volumetric brain registration using harmonic mappings,” Medical
Imaging, IEEE Transactions on,vol.26,no.12,pp.1657–1669,2007.
[112] A. A. Joshi, D. W. Shattuck, and R. M. Leahy, “A method for automated cortical
surface registration and labeling,” in Biomedical Image Registration, pp. 180–
189,Springer,2012.
[113] O. L. Munk, L. Bass, K. Roelsgaard, D. Bender, S. B. Hansen, and S. Keiding,
“LiverkineticsofglucoseanalogsmeasuredinpigsbyPET:importanceofdual-
inputbloodsampling,” J Nucl Med,vol.42,no.5,pp.795–801,2001.
[114] S.Y.Amy,H.-D.Lin,S.-C.Huang,M.E.Phelps,andH.-M.Wu,“Quantification
of cerebral glucose metabolic rate in mice using 18f-fdg and small-animal pet,”
Journal of Nuclear Medicine,vol.50,no.6,pp.966–973,2009.
135
[115] G. Z. Ferl, X. Zhang, H.-M. Wu, and S.-C. Huang, “Estimation of the 18f-fdg
input function in mice by use of dynamic small-animal pet and minimal blood
sampledata,”JournalofNuclearMedicine,vol.48,no.12,pp.2037–2045,2007.
[116] A. Rahmim and H. Zaidi, “Pet versus spect: strengths, limitations and chal-
lenges,” Nuclear medicine communications,vol.29,no.3,pp.193–207,2008.
[117] Y.Zhou,W.Ye,J.R.Braˇ si´ c,A.H.Crabb,J.Hilton,andD.F.Wong,“Aconsistent
andefficientgraphicalanalysismethodtoimprovethequantificationofreversible
tracerbindinginradioligandreceptordynamicpetstudies,”Neuroimage,vol.44,
no.3,pp.661–670,2009.
[118] K.-P. Wong, V. Kepe, M. Dahlbom, N. Satyamurthy, G. W. Small, J. R. Bar-
rio, and S.-C. Huang, “Comparative evaluation of logan and relative-equilibrium
graphical methods for parametric imaging of dynamic [¡ sup¿ 18¡/sup¿ f] fddnp
petdeterminations,” NeuroImage,vol.60,no.1,pp.241–251,2012.
[119] A. Rahmim, Y. Zhou, J. Tang, L. Lu, V. Sossi, and D. F. Wong, “Direct 4d para-
metricimagingforlinearizedmodelsofreversiblybindingpettracersusinggen-
eralized ab-em reconstruction,” Physics in medicine and biology, vol. 57, no. 3,
p.733,2012.
[120] C. Byrne, “Iterative algorithms for deblurring and deconvolution with con-
straints,” Inverse Problems,vol.14,no.6,p.1455,1998.
[121] G. Wang and J. Qi, “Direct estimation of kinetic parametric images for dynamic
pet,” Theranostics,vol.3,no.10,p.802,2013.
[122] W. Sureshbabu and O. Mawlawi, “Pet/ct imaging artifacts,” Journal of nuclear
medicine technology,vol.33,no.3,pp.156–161,2005.
[123] G. K. von Schulthess, Molecular anatomic imaging: PET-CT and SPECT-CT
integrated modality imaging. LippincottWilliams&Wilkins,2007.
[124] T.Beyer,D.W.Townsend,T.Brun,P.E.Kinahan,M.Charron,R.Roddy,J.Jerin,
J.Young,L.Byars,R.Nutt,etal.,“Acombinedpet/ctscannerforclinicaloncol-
ogy,” Journal of nuclear medicine,vol.41,no.8,pp.1369–1379,2000.
136

Abstract (if available)

Abstract We investigate using list‐mode PET data to perform Patlak and Logan modeling. We first propose an estimation method for irreversible tracers (where Patlak modeling is applied) for dynamic PET studies in which we compute voxel‐wise estimates of Patlak parameters using two frames of data. Our approach directly uses list‐mode arrival time for each event to estimate the Patlak parametric image. We use a penalized likelihood method in which the penalty function uses spatially variant weighting to ensure a count independent local impulse response. We evaluate performance of the method in comparison to fractional changes in standardized uptake value (%DSUV) between the two frames using Cramer‐Rao analysis and Monte Carlo simulation. Receiver operating characteristic (ROC) curves are used to compare performance in differentiating tumors from background based on the dynamic data sets. Using area under the ROC curve as a performance metric, we show superior performance of Patlak compared to %DSUV over a wide range of dynamic data sets and parameters. These results suggest that Patlak analysis can be used for dynamic dual‐frame PET data and may lead to better quantitative results relative to %DSUV method. ❧ We then extend our method to clinical wholebody patient data with ¹⁸F‐FDG, which is a commonly used irreversible tracer for oncologic applications. Several practical issues such as estimation of dynamic randoms and scatters rate functions, motion correction, and the estimation of blood input function are discussed and effectively handled. ❧ A new whole‐body PET scan protocol is proposed, where the patient is stepped through the PET scanner twice, after which the Patlak parametric image for each bed position is estimated from the listmode data of the two frames and stitched together to form a wholebody Patlak parametric image. Ten patients were scanned using this protocol. We then draw liver, muscle, and spine regions of interest (ROIs) for these patients and the mean and standard deviation (s.d.) for the Patlak slope (net influx rate) and intercept are computed. We also compute the SUV values and the s.d. for the same ROIs at the 2 corresponding frames time, after which the fractional SUV (%DSUV) is calculated. Results show that Patlak method has less inter-subject variation and higher tumor‐to‐background contrast. Overall, the Patlak estimates provide quantitative parameters better describing tracer uptake than SUV. This may be useful for cancer staging, treatment planning and assessing response to therapy. ❧ We also investigate the direct estimation of Logan parameters (where a reversible tracer is used) from listmode PET data. The ill‐conditioning nature of listmode reconstruction for PET data from a reversible radioactive tracer is explored and discussed. We derive a modified scheme improving the conditioning of the problem by utilizing data before the steady state. This method differs from a previously proposed sinogram based reconstruction methods by other researchers in that it follows the true likelihood of the data. Cramer‐Rao analysis and Monte Carlo simulation demonstrate the effectiveness of our method. ❧ Nowadays many PET/CT scans are performed with contrast agents aiming at enhancing soft tissue contrast in CT images. The contrast agent will increase the CT value significantly especially in the blood rich soft tissues such as liver and heart. However the liner attenuation coefficient of contrast agent is only slightly higher than water at 511 keV. As a result, using the same mapping from low to high energies for non‐contrast CT will result in over‐estimation of attenuation in regions contains contrast region for PET. ❧ This leads to inaccuracies in PET images. We describes a joint estimation algorithm for emission and attenuation maps using contrast enhanced CT as a prior. A comparison with simple rescaling scheme shows that our method can improve quantitative image accuracy when using CT with contrast to compute attenuation coefficients.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

PDF

Model based PET image reconstruction and kinetic parameter estimation

PDF

Measuring functional connectivity of the brain

PDF

Transmission tomography for high contrast media based on sparse data

PDF

New theory and methods for accelerated MRI reconstruction

PDF

Causality and consistency in electrophysiological signals

PDF

Automatic image and video enhancement with application to visually impaired people

PDF

Signal processing methods for interaction analysis of functional brain imaging data

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography

PDF

Fast flexible dynamic three-dimensional magnetic resonance imaging

PDF

Advanced coronary CT angiography image processing techniques

PDF

Statistical signal processing of magnetoencephalography data

PDF

Signal processing for channel sounding: parameter estimation and calibration

PDF

Functional connectivity analysis and network identification in the human brain

PDF

Efficient template representation for face recognition: image sampling from face collections

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Structure learning for manifolds and multivariate time series

PDF

Non‐steady state Kalman filter for subspace identification and predictive control

PDF

Facial age grouping and estimation via ensemble learning

Asset Metadata

Creator Zhu, Wentao (author)

Core Title Direct wholebody Patlak and Logan image estimation from listmode PET data

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 08/11/2014

Defense Date 05/20/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag direct,listmode,Logan,OAI-PMH Harvest,Patlak,wholebody

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Leahy, Richard M. (committee chair)

Creator Email zhuwentao.ee@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-459551

Unique identifier UC11287765

Identifier etd-ZhuWentao-2802.pdf (filename),usctheses-c3-459551 (legacy record id)

Legacy Identifier etd-ZhuWentao-2802.pdf

Dmrecord 459551

Document Type Dissertation

Format application/pdf (imt)

Rights Zhu, Wentao

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA