<%BANNER%>

Design and Evaluation of Reconstruction Methods Incorporating Estimated Motion in Gated Cardiac Emission Tomography

Permanent Link: http://ufdc.ufl.edu/UFE0022821/00001

Material Information

Title: Design and Evaluation of Reconstruction Methods Incorporating Estimated Motion in Gated Cardiac Emission Tomography
Physical Description: 1 online resource (102 p.)
Language: english
Creator: Parker, Jason
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2008

Subjects

Subjects / Keywords: cardiac, channelized, image, motion, spect
Nuclear and Radiological Engineering -- Dissertations, Academic -- UF
Genre: Nuclear Engineering Sciences thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: In medical imaging, artifacts due to patient motion can make accurate diagnoses more difficult. The cardiac and respiratory cycles are common sources of patient motion in cardiac emission tomography (ET). For a cyclical motion (such as cardiac or respiratory motion), a gating procedure may be used to freeze the acquisition at the different phases of the cycle, effectively removing the motion blur. However, this reduction in motion blur is at the cost of higher noise levels in the gated images. Our purpose was to develop and test new methods of image reconstruction in cardiac ET that use estimates of patient motion between image frames to improve the noise characteristics of the gated images. The primary method of reconstruction investigated here was one which simultaneously estimates the tomographic images and the frame-to-frame cardiac contraction motion vector fields in a mutually influential manner. The effect this simultaneous approach has on the accuracy of the estimated motion vector fields was tested against a conventional motion estimation method using a physical, dynamic phantom. The simultaneous method was found to have improved motion estimation accuracy compared to the conventional method. Furthermore, a new method for estimating and correcting the respiratory motion of the heart was developed and compared to several other methods proposed in the literature. The new method shows improved speed, accuracy, and robustness compared to other proposed methods. Finally, the simultaneous reconstruction method was evaluated against a standard reconstruction method in clinical use using a mathematical observer formulated to operate on both single- and multi-frame images, as well as a signal-to-noise ratio measurement using regions-of-interest. In both mathematical observer studies the difference between the simultaneous method and the standard method were not statistically significant. However, the signal-to-noise ratio calculation using regions-of-interest found the simultaneous method to have improved detection characteristics compared to the standard method. Results were found to be significant at the 1% level.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Jason Parker.
Thesis: Thesis (Ph.D.)--University of Florida, 2008.
Local: Adviser: Gilland, David R.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2008
System ID: UFE0022821:00001

Permanent Link: http://ufdc.ufl.edu/UFE0022821/00001

Material Information

Title: Design and Evaluation of Reconstruction Methods Incorporating Estimated Motion in Gated Cardiac Emission Tomography
Physical Description: 1 online resource (102 p.)
Language: english
Creator: Parker, Jason
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2008

Subjects

Subjects / Keywords: cardiac, channelized, image, motion, spect
Nuclear and Radiological Engineering -- Dissertations, Academic -- UF
Genre: Nuclear Engineering Sciences thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: In medical imaging, artifacts due to patient motion can make accurate diagnoses more difficult. The cardiac and respiratory cycles are common sources of patient motion in cardiac emission tomography (ET). For a cyclical motion (such as cardiac or respiratory motion), a gating procedure may be used to freeze the acquisition at the different phases of the cycle, effectively removing the motion blur. However, this reduction in motion blur is at the cost of higher noise levels in the gated images. Our purpose was to develop and test new methods of image reconstruction in cardiac ET that use estimates of patient motion between image frames to improve the noise characteristics of the gated images. The primary method of reconstruction investigated here was one which simultaneously estimates the tomographic images and the frame-to-frame cardiac contraction motion vector fields in a mutually influential manner. The effect this simultaneous approach has on the accuracy of the estimated motion vector fields was tested against a conventional motion estimation method using a physical, dynamic phantom. The simultaneous method was found to have improved motion estimation accuracy compared to the conventional method. Furthermore, a new method for estimating and correcting the respiratory motion of the heart was developed and compared to several other methods proposed in the literature. The new method shows improved speed, accuracy, and robustness compared to other proposed methods. Finally, the simultaneous reconstruction method was evaluated against a standard reconstruction method in clinical use using a mathematical observer formulated to operate on both single- and multi-frame images, as well as a signal-to-noise ratio measurement using regions-of-interest. In both mathematical observer studies the difference between the simultaneous method and the standard method were not statistically significant. However, the signal-to-noise ratio calculation using regions-of-interest found the simultaneous method to have improved detection characteristics compared to the standard method. Results were found to be significant at the 1% level.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Jason Parker.
Thesis: Thesis (Ph.D.)--University of Florida, 2008.
Local: Adviser: Gilland, David R.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2008
System ID: UFE0022821:00001


This item has the following downloads:


Full Text

PAGE 1

1 DESIGN AND EVALUATION OF RECONSTRUCTION METHODS INCORPORATING ESTIMATED MOTION IN GATED CARDIAC EMISSION TOMOGRAPHY By JASON G. PARKER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008

PAGE 2

2 2008 Jason G. Parker

PAGE 3

3 To the sick.

PAGE 4

4 ACKNOWLEDGEMENTS I thank m y Mom and Dad, my sisters Jessica and Julia, and my brother Jimbo, for their love and support. I thank the University of Fl orida for the rich educational opportunities it has presented me, and the State of Florida for developing a strong and influential institution. I thank Dr. Shukla for his insight and patience thr oughout many time-consuming experiments in his department. I attribute the majority of my clinical knowledge to the time I spent under his supervision at the Malcolm Randall Veterans Affa irs Medical Center in Gainesville, FL. I thank Dr. Hintenlang for his comments and direction on this work, as well as for three years of teaching in the fundamentals of Medical Physic s. His broad body of knowledge and expertise plays a key role in the University of Florida Medical Physics program. I thank Dr. Mair for his patient teaching style and encouragement. He has b een an invaluable contributor to my education and understanding of the scientific method. It has been an honor work beside him. Finally, I thank my advisor, mentor, and friend, Dr. Gilland. He has treated me fairly and with respect, while insisting on only my best performance. He challenged me in every aspect of my work and in doing so, spurred a wonderful period of intellectual growth for myself.

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGEMENTS .............................................................................................................4LIST OF TABLES ...........................................................................................................................7LIST OF FIGURES .........................................................................................................................8ABSTRACT ...................................................................................................................... .............10 CHAP TER 1 SPECIFIC AIMS ................................................................................................................. ...122 BACKGROUND AND SIGNIFICANCE .............................................................................. 14Simultaneous Image Reconstruction and Motion Estimation ................................................ 16Respiration-Induced Motion Artifacts ....................................................................................17Channelized Hotelling Observer .............................................................................................193 WALL MOTION ESTIMATION FO R GATED CARDIAC EMISSION TOMOGRAPHY: PHYSICAL PHANTOM EVALUATION ...............................................21Introduction .................................................................................................................. ...........21Methods ..................................................................................................................................22Data Acquisition .............................................................................................................. 22Marker Motion Estimation ..............................................................................................24Myocardium Motion Estimation .....................................................................................242-Frame Motion Estimation ............................................................................................27Motion Estimation Error ..................................................................................................27Summing Estimated Motion Vectors Across Frames .....................................................27Results .....................................................................................................................................288-Frame Data .................................................................................................................. .282-Frame Data .................................................................................................................. .29Summing Estimated Motion Vectors Across Frames .....................................................29Summary and Discussion .......................................................................................................294 IMPROVED RESPIRATORY MOTION CORRECTION IN GATED CARDIAC SPECT ......................................................................................................................... ...........36Introduction .................................................................................................................. ...........36Theory ........................................................................................................................ .............39Experimental Methods .......................................................................................................... ..41Algorithms .................................................................................................................... ...41Mathematical Phantom Evaluation .................................................................................. 43Respiratory-Gated Patient SPECT Images ...................................................................... 47

PAGE 6

6 Results .....................................................................................................................................48Motion Estimation on the Noise-Free Phantom Frames .................................................48Effect of Segmentation ....................................................................................................49Comparison of Methods at Optimal Segmentation ......................................................... 50Effect of Varying Level of Extra-Myocardial Activity ................................................... 51Human SPECT Image .....................................................................................................51Summary and Discussion .......................................................................................................525 SIMULTANEOUS IMAGE RECONSTRUCTION AND M OTION ESTIMATION IN GATED CARDIAC EMISSION TOMOGRAPHY: A MATHEMATICAL OBSERVER STUDY .............................................................................................................60Introduction .................................................................................................................. ...........60Background .................................................................................................................... .........61Motion Compensated Image Reconstruction ..................................................................61Simultaneous Image Reconstruction and Motion Estimation ......................................... 62Channelized Hotelling Observer .....................................................................................63Multiple-Array Channelized Hotelling Observer ............................................................ 65Methods ..................................................................................................................................67Simulation of Test Data ................................................................................................... 67Validation of the CHO .....................................................................................................68Image reconstruction ................................................................................................ 70ROC evaluation ........................................................................................................ 71Direct assessment of SNR using regions-of-interest across the ensemble data .......71Reconstruction without the effects of scatter attenuation, and detector response ................................................................................................................ 72Results .....................................................................................................................................72Validation of the CHO .....................................................................................................72Image Reconstruction ...................................................................................................... 73ROC Evaluation ............................................................................................................... 74Direct Assessment of SNR Using Regions-of-Interest acr oss the Ensemble Data ......... 75Reconstruction without the Effects of Atte nuation, Scatter, and Detector Response ..... 75Summary and Discussion .......................................................................................................766 CONCLUSIONS ................................................................................................................... .90APPENDIX A ROTATION MATRIX AND PARTIAL DERIVATIVES .................................................... 92B DETERMINATION OF ROTATION MATRIX ................................................................... 93C PARTIAL DERIVATIVES USED IN CGQ ALGORITHM ................................................. 94REFERENCES .................................................................................................................... ..........95BIOGRAPHICAL SKETCH .......................................................................................................102

PAGE 7

7 LIST OF TABLES Table page 3-1 Average estimated motion errors for 8-fram e data (pixel units) ........................................353-2 Average estimated motion errors fo r two-frame data (pixel units) ....................................353-3 Summed motion errors for OSEM-HS and RM ................................................................354-1 Relative organ activity concentrations for the simulated NCAT phantom ........................ 594-2 Respiratory motion in the mathematical phantom images relative to frame 1. Translations are given in pi xels, rotations in degrees. ....................................................... 594-3 Motion estimates from the noise-free phantom images ..................................................... 594-4 Average motion estimation errors for the 10 noisy simulated images with background activity at optimal segmentation .................................................................... 595-1 Relative organ activity concentrations per pixel for the simulated phantom ..................... 895-2 Temporal smoothing kernels ..............................................................................................89

PAGE 8

8 LIST OF FIGURES Figure page 3-1 Dynamic phantom on imaging bed with close-up of radioactive m arkers shown on right. Outer shell of phantom is removed. ........................................................................323-2 Projection images obtained from a si ngle angle of the marker acquisition. ......................323-3 Low noise projection images obtained from a single angle of the myocardium acquisition. .................................................................................................................. .......323-4 Single noise realization obtained after scaling and a dding Poisson noise to the projection images in Fig. 3-3. ............................................................................................323-5 Reconstructed images obtained from the marker acquisition ............................................333-6 Sum-of-Squared Errors vs Cut-off frequency for 5 (left) 10 (middle), and 50 (right) iterations of OSEM ............................................................................................................333-7 Transaxial slices of frame 1 of the low noise OSEM, noisy OSEM, and noisy RM reconstructions ...................................................................................................................333-8 Euclidean distance vs. Beta for OSEM-HS (left) and RM (right). ....................................343-9 Euclidean distance (ensemble average) vs frame number for markers 2, 4, 6, and 10. ....344-1 Select segmentation levels. ............................................................................................... .544-2 Transaxial slice of th e segmented human image used for motion estimation. .................. 544-3 Select slices of the phant om images at end-diastole. ......................................................... 554-4 Horizontal (left) and short (rig ht) axis profiles for Figure 4-3. .........................................554-5 Results for average PME found by varying the segmentation parameters. ....................... 564-6 Reconstructions of the noisy simulated images at end-diastole......................................... 564-7 Horizontal (left) and short (r ight) axis profiles for Fig. 4-6. .............................................574-8 Results for average PME found by varyi ng the background activity concentration at optimal segm entation. ........................................................................................................ 574-9 Reconstructions of the simulated data at extra-myocardial activ ity levels of 0.25 and 1.75 (fraction of average). ..................................................................................................584-10 Reconstructions of the human images with CGQ correcti on, principal axes correction, and no correction, respectively. ....................................................................... 58

PAGE 9

9 5-1 Four frequency bands used in the CHO. ............................................................................785-2 Spatial and frequency domain templates of the four bands used in the CHO. .................. 785-3 Phantom source distributions showing the shape and location of the defect. .................... 795-4 Noise-free reconstructions for the normal (left) and abnormal (right) datasets. ................ 795-5 Regions-of-interest used to calculate SNR directly. .......................................................... 795-6 D-prime vs defect contrast. ................................................................................................805-7 D-prime vs count level. ................................................................................................... ...805-8 D-prime vs number of channels used. ................................................................................ 815-9 D-prime vs number of training images used ...................................................................... 815-10 D-prime vs Butterworth cut-off frequency ........................................................................ 825-11 Select images us ed in Figure 5-10. ....................................................................................825-12 Results of varying the temporal smoothing kernel number. .............................................. 835-13 Example images of the five levels of temporal smoothing evaluated in Test 6. ............... 835-14 Full-width-half-max vs MLEM iteration number on noise-free data. ............................... 845-15 Full-width-half-max vs Butterworth cut-off frequency on noise-free data. ...................... 845-16 Full-width-half-max vs temporal convolution kernal on noise-free data. .........................855-17 Example slices of the PS-OSEM and RM reconstructed images. ..................................... 855-18 Motion vector fields for RM reconstructions in Fig. 5-15 ................................................. 865-19 Example RM and PS-OSEM images used in the ROC study. ........................................... 875-20 Reciever operating characteristiscs curve for RM and PS-OSEM using 2D CHO on a single lesion. ................................................................................................................ ......875-21 Reciever operating characteristiscs curve for RM and PS-OSEM using MACHO on a single lesion. ................................................................................................................ ......885-22 Results of SNRROI vs. 3D Butterworth cut-off frequency. ................................................885-23 Reconstructions of projection data obtained using a simple linear projection operator. ... 89

PAGE 10

10 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DESIGN AND EVALUATION OF RECONSTRUCTION METHODS INCORPORATING ESTIMATED MOTION IN GATED CARDIAC EMISSION TOMOGRAPHY By Jason G. Parker December 2008 Chair: David R. Gilland Major: Nuclear Engineering Sciences In medical imaging, artifacts due to patie nt motion can make accurate diagnoses more difficult. The cardiac and respiratory cycles ar e common sources of patient motion in cardiac emission tomography (ET). For a cyclical moti on (such as cardiac or respiratory motion), a gating procedure may be used to freeze the acqui sition at the different phases of the cycle, effectively removing the motion blur. However, this reduction in motion blur is at the cost of higher noise levels in the gated images. Our pur pose was to develop and test new methods of image reconstruction in cardiac ET that use estimates of patient motion between image frames to improve the noise characteristics of the gated images. The primary method of reconstruction investigated here was one wh ich simultaneously estimates the tomographic images and the frame-to-frame cardiac contraction motion vector fields in a mutually influential manner. The effect this simultaneous approach has on the accu racy of the estimated mo tion vector fields was tested against a conventional motion estimati on method using a physical, dynamic phantom. The simultaneous method was found to have improved motion estimation accuracy compared to the conventional method. Furthermore, a new method for estimating and correcting the respiratory motion of the heart was developed and compared to several other methods proposed in the

PAGE 11

11 literature. The new method shows improved spee d, accuracy, and robustness compared to other proposed methods. Finally, the simultaneous re construction method was evaluated against a standard reconstruction method in clinical use using a mathematical observer formulated to operate on both singleand multi-frame images, as well as a signal-to-noise ratio measurement using regions-of-interest. In both mathematical observer studies the difference between the simultaneous method and the standard method were not statistically significant. However, the signal-to-noise ratio calculati on using regions-of-interest found the simultaneous method to have improved detection characteristics compared to the standard method. Results were found to be significant at the 1% level.

PAGE 12

12 CHAPTER 1 SPECIFIC AIMS Single photon computed emission tomogra phy (SPECT) myocardial perfusion imaging (MPI) is an important tool in th e diagnosis and risk stratification of coronary artery disease. When using the common radiotracer technitium-9 9m-sestamibi, the 3D images are indicators of myocardial blood flow. The additi on of electrocardiographic (E CG) gating during the acquisition procedure allows further analysis of the beating motion of the heart, such as the determination of ejection fraction, wall motion abnormalities, and wall thickening. This motion information, together with the perfusion images, can lead to an assessment of viability in regions of decreased uptake. A negative impact of the gating procedure, however, is that the individual gated images have higher noise levels, impacting both the quality of the images and the accuracy of the motion measurements. The overall goal of this research is to develop and evaluate software methods which improve both the images and motion inform ation in gated SPECT MPI, while maintaining the clinically relevant characterist ics of speed, accuracy, and robustness. In SPECT MPI, there has been considerable research in both im age reconstruction and motion estimation; however, the two have rarely been viewed as a single problem. The primary reconstruction method investigated in this wo rk is one which simultaneously estimates the reconstructed image frames and the frame-to-fra me non-rigid contraction motion in a mutually influential manner. The impact of this appr oach on image quality and motion estimation, and ultimately diagnostic accur acy, were investigated. Specific Aim 1: To test the wall motion estimation accuracy of a simultaneous reconstruction and motion estimation method over the cardiac contraction cycle compared to a conventional motion estimation algorithm. The hypothesis that a joint reconstruction me thod is a better method than estimating the images and cardiac contraction motion vector fiel ds separately rests on the assumption that the

PAGE 13

13 motion estimates should improve as the images improve. This hypothesis was tested using point source markers on a physical, dynamic cardiac phantom. Specific Aim 2: To develop and test a new method for respiratory mo tion compensation which addresses the clinically releva nt issues of speed, accuracy, and robustness. Motion blur due to the respiratory motion of the heart can blur the myocardial walls and may decrease the diagnostic accuracy of the imaging procedure. In this work we developed an iterative method to correct this motion which treats the heart as a rigid body. The method was tested against other methods in the SPECT lite rature using both simulated and human data. Specific Aim 3: To test the detection characteristi cs of the simultaneous method compared to a standard reconstructi on method in clinical use. The simultaneous method was compared to a standard reconstruction method in clinical use using the channelized Hote lling observer (CHO). The CHO can be thought of as a measure of the signal-to-noise ratio and has been shown to closely model the detection capabilities of human observers.

PAGE 14

14 CHAPTER 2 BACKGROUND AND SIGNIFICANCE Heart disease is the leading cause of deat h in the United States [1], and is most commonly characterized by coronary artery disease (CAD). Common initial diagnostic procedures for patients with possible CAD in clude stress electrocardiography, also known as exercise tolerance testing (ETT), SPECT m yocardial perfusion imaging, and stress echocardiography. Analyses of these modalities gene rally finds the sensitivity and specificity of SPECT MPI and stress echocardiography to be compar able to each other, while both are superior to ETT [2]-[4]. Cost-effectiveness analyses ha ve found SPECT MPI to be a better initial diagnostic tool than ETT or stress echocardiography [5], and this has propelled SPECT MPI to be one of the most important tools in diagnosing patients with heart disease. Typical SPECT MPI procedure calls for the pa tient to be injected with a gamma emitting radionuclide, such as thallium-201 or technetium-99m (99mTc), and subsequently imaged using a rotating gamma camera. 99mTc is a more favorable radionuclide than thallium-201 in terms of detection, due to the higher energy (140 keV vs 68-80.3 keV) of the emitted gamma and shorter radiological half-life (6.0 hrs vs. 73.1 hrs). In the case of SPECT MPI with 99mTc, the sestamibi molecule is attached to the radionucl ide, forming the common radiotracer 99mTc-sestamibi. Shortly after injection of 99mTc-sestamibi, the radiotracer accumul ates within the myocardium in a distribution proportional to blood flow. Unlike thallium-201, 99mTc-sestamibi crosses cell membranes via passive transport, and thus ar eas of low uptake genera lly do not change over time, making 99mTc-sestamibi a good indicator of blood flow but less of an indicator of function [6]. An important improvement to SPECT MPI wa s the introduction of electrocardiographic (ECG) gating [7]. In this procedure, a live ECG si gnal from the patient is used to bin detected

PAGE 15

15 events into discrete phases of the cardiac contraction cycle. This eliminates much of the motion blur due to the beating motion, but more importantly enables quantitative determination of important functional information re garding the health of the hear t, specifically, ejection fraction (EF), regional wall motion abnormalities, and wall thickening. The addition of this functional information has been reported to increase the se nsitivity of SPECT MPI to greater than 95% and the specificity as high 94% [2]. A negative impact of ECG gating is that the individual gated images have a decreased SNR, due to the distribution of events over many time frames (typically 8). This decrease in SNR has motivated research into reconstruction met hods which take into account the relationship between the images in the time series. In th e simplest approach, one may view the time dimension as any other dimension, and the images as a single 4D image. It is then reasonable to perform any number of filter operations across the time dimension that one would use in the spatial dimension (e.g. Gaussian convolution, Fourie r filtering). This has th e effect of improving the noise characteristics of the images; however, it is at the ex pense of increased motion blur. Reconstruction followed by spatial and temporal smoothing is referred to as a post-smoothed reconstruction and is a typical approach in the clinical setting due to its simplicity. In Lalush and Tsui [8] several different dist ributions were evaluated for use as temporal Gibbs priors in a maximum a posteriori (MAP ) reconstruction. The authors found that the different functions and the diffe rent weighting schemes they applied controlled the extent to which contrast, noise, and motion blur were pres ent, but that ultimately the choice of these parameters was study specific. It was further shown in Lalush et al. [9] that by including a nonrigid motion estimation in the prior term, the image noise could be improved as in [8], but with a

PAGE 16

16 lower level of motion blur. A fundamental problem with this approach, however, is that the accuracy of the motion estimate suffers from the already high noise levels of the images. Simultaneous Image Reconstruc tion and Motion Estimation A m ore recent advancement in gated SPE CT MPI reconstruction research was the introduction of the simultaneous image recons truction and motion estimation approach. The simultaneous method rests on two assumptions: 1) a better motion estim ate leads to a better reconstructed image, and 2) a be tter reconstructed image leads to a better motion estimate. The method is implemented iteratively with alternating updates to the reconstructed images and estimation motion vector fields. In Gilland et al. [10] the first use of the simultaneous approach was demonstrated for a simple two-frame system. A conjugate gradient (CG) routine was used to minimize an objective function consisting of three terms: a negativ e log-likelihood term (reconstruction), an image matching term (motion), and a constraint term on the estimated motion consisting of a biomechan ical model of the strain en ergy of an elastic material undergoing deformation (here representing the m yocardium). This paper showed favorable results for the simultaneous method over standard maximum likelihood-expectation maximization (MLEM) [11] reconstructions in a sum-of-squared errors (SSE) sense. The accuracy of the estimated motion vectors was not investigated. A generalization of the method in [10] was presented in Mair et al. [12]. The method, now referred to as RM (for reconstruction and moti on), consisted of two steps: 1) the R-step, in which an image registration objective func tion was minimized while holding the motion constant, and 2) the M-step, in which the objective function was minimized holding the images constant. The method was generalized to any numb er of frames, and now had the effect of any one frame being influenced by the frame immedi ately before and the frame immediately after itself. The accuracy of the reconstructed images wa s evaluated again in an SSE sense, as well as

PAGE 17

17 qualitatively using visual insp ection and profile anal ysis. The study found that the simultaneous method produced images with better noise proper ties than those generated with a post-smoothed (PS) ordered subsets-expecta tion maximization (OSEM) [13] method. The RM-generated motion was evaluated quantitatively us ing the phantom matching error (PME), a measure of how well the frames matched after deformation by the es timated motion. The authors found that the PME did improve with iteration, conf irming the hypothesis that as the images improve the motion estimates improve. Despite these findings, several problems still remain with RM. The evaluation of the images in Mair et al. [12] does not confirm that the RM generated images are a better diagnostic tool than the current standard, as the SSE does not necessarily reflect the ability of a human to detect lesions. Furthermore, the PME is an indi rect measure of motion accuracy and its value is susceptible to noise, image artifacts, and inter polation errors. A direct assessment of the motion accuracy is needed if the RM motion is to be used in a diagnostic environment. Finally, the algorithm is inherently slow. In [12] the authors reported a co mputational time of 2h 50min (in comparison to a PS-OSEM computational time of ~30s ). This is not acceptable for clinical use, and methods for reducing the reconstruction time to < 1min need to be investigated. Respiration-Induced Motion Artifacts For imaging modalities requiring scan times longer than a single breath-hold (e.g. SPECT, PET, MR), the motion of the heart due to respiration can introduce blur into the reconstructed images [14]-[16]. Similar to ECG-gating of the cardiac contraction cycle, a gating procedure can be introduced to divide detect ed events up between pha ses of the respiratory cycle. If the respiratory motion of the heart between the gated im ages can be determined, much of the respiration-induced blur can be removed by registering each frame to a reference frame, and summing. Methods for estimating the frame-to-frame respiratory motion of the heart are an

PAGE 18

18 active area of research, with the proposed methods generally falling into one of two categories: external tracking systems, and data analysis methods. A brief review of several common methods is presented here. External tracking systems use hardware to track patient motion throughout an imaging procedure. One example of this approach was dem onstrated in Bruyant et al. [17] in which the authors proposed an external visual tracking system (VTS) to track patient motion during listmode SPECT acquisition. An elastic sleeve was out fitted with reflective spheres and fit to the patients chest. The motion of the spheres was tracked throughout the ac quisition by five video cameras. After the acquisition, the li st-mode data could be adjusted to remove much of the blur due to patient motion. The authors found that th e system was highly accurate and robust, and suggested its use for correcting the respirator y motion of the heart in gated SPECT imaging. However, it is known that the motion of mark ers placed outside the body may not reflect the motion of deeply-embedded structures within the body. Furthermore, external tracking systems may be expensive and require regular quality assu rance. Also, currently a VTS cannot be shared between multiple scanners, meaning each scanne r in a facility that will use VTS motion correction will require a co mplete VTS installation. In contrast, methods for motion correction ba sed on data analysis can be shared by many scanners, do not require maintenance, and are generally less expensive than external tracking systems. Furthermore, data analysis methods ha ve the advantage of meas uring the motion of the organ directly. An example of a data analysis me thod was detailed in Kovals ki et al. [18] where the authors proposed modeling the motion of th e heart as a rigid body motion during list-mode SPECT acquisition. The respirator y cycle was monitored concurre nt to the SPECT acquisition using a piezo-electric chest belt. The data we re binned post-reconstruction according to the

PAGE 19

19 respiratory signal and each bin was reconstructed. The motion between the resulting temporal frames was estimated analytically by fitting ellips oids to the myocardial activity in each frame. This motion was then used to shift each frame ba ck to a reference frame. The resulting images were found to have better spat ial resolution than images without correction. However, it was found that often the improvement was small and, fr om a clinical standpoint, negligible. It is known that the accuracy of data analysis met hods can be highly dependant on image noise, myocardial segmentation, and computation time, and the authors suggest that new methods designed to overcome these problems may lead to more significant respir atory motion correction in SPECT MPI. Channelized Hotelling Observer In evaluating new reconstruction methods it is perhaps most important to understand how a new method improves (or degrades) lesion dete ction compared to a standard method. The conventional approach to quantifying imaging sy stem detection performa nce is by the receiver operating characteristics (ROC) methodology [20]. In ROC analysis, human observers are asked to rate a large number of images based on how su re they are that each image contains a lesion. Their responses are used to generate plots of true positive fraction (TPF) vs. false positive fraction (FPF), often referred to as ROC curves. The area under the ROC curve may be used as an image quality figure of merit. A drawback to ROC analysis is that the studies can be time consuming and often rely on the pa rticipation of many observers. An alternative to performing human ROC studi es is to use a mathematical model of lesion detection to evaluate system performan ce. The Hotelling observer (HO) is one such mathematical model for binary detection tasks which has been shown to correlate well with human observer performance [19]. The HO uses se ts of training images to develop a covariance matrix which describes the statistical difference between signal absent (SA) and signal

PAGE 20

20 present (SP) images. After the ini tial training, the HO can be used to generate a test statistic for a test image which can then be compared to a deci sion threshold used to classify the current image as SA or SP. Thus, a large number of testing images may be used to mimic a human ROC study. In Yao and Barrett [21] the HO model was modified with a channel mechanism to imitate the frequency selective bands which human obser vers have been shown to operate under [29]. This extension of the HO, the channelized Hotelling observer (CHO) was found to be a better predictor of human observer perf ormance than the HO. The CHO has been used extensively in the SPECT MPI literature to both optimize sy stem performance and to compare new imaging methods.

PAGE 21

21 CHAPTER 3 WALL MOTION ESTIMATION FOR GATED CAR DIAC EMISSION TOMOGRAPHY: PHYSICAL PHANTOM EVALUATION1 Introduction In cardiac em ission tomography (ET), improved image quality has been demonstrated when estimated wall motion (non-rigid) is incorp orated into 4D recons truction methods [12], [23]-[24]. One approach is to estimate the mo tion vector field from an initial filtered backprojection image and use the estimate in the penalty term of a penaliz ed ML approach [24]. Several methods have been proposed for estimati ng the motion vector fields from reconstructed images, typically using optical flow or elastic de formation models [12], [ 25]-[29]. In Mair et al. [12] it was shown that simultaneously estimating the reconstructed images and motion vector fields in each iteration resulted in improved image noise. Since non-rigid motion estimation is play ing an important role in cardiac ET reconstruction, methods for evaluating motion estim ation accuracy have become important. The objective of this study is to evaluate the accur acy of wall motion estimation methods for gated cardiac ET using a physical, dynamic phantom. Th is extends our earlier work, which used simulated data of a mathematical phantom with a simple system model that excluded the effects of scatter, detector response, and attenuation. Our approach here is to us e high count images of point source markers attached to the myocardial surface of a dyna mic thorax phantom to provide an independent measure of the true myocardial motion. As in the earlier work, we compare estimated motion and reconstructed im age quality obtained from a simultaneous reconstruction/motion estimation method (RM method) [12] with th at obtained from a 1 IEEE 2008. Reprinted with permission from (J.G. Parker, D.R. Gilland, Wall motion estimation for gated cardiac emission to mography: physical phantom evaluation, IEEE Trans. Nuc. Sci. vol. 55, no. 1, pp. 531-536, 2008.).

PAGE 22

22 conventional optical flow method [2 9] applied to optimized OSEM [3 0] images. In this work, we will refer to this conventional approach as the standard method. One advantage of the physical phantom, relative to a mathematical phantom, is that the physics of photon propagation and detection are assure d to be realistically incorporated into the data. Also, the point marker vali dation approach enables the accu rate tracking of material points on the myocardium. One of the disadvantages of the physical phantom relates to the differences in motion characteristics compared with the typi cal human heart. For example, unlike the human heart, this phantom myocardium does not exhibit twisting, or wringing, moti on, and the base is stationary during contraction [31]. The wall thic kening pattern may not match the typical human case. A second disadvantage in this study is that, in order to ensure co -registration with the marker data, no extra-myocardial activity has been included (details in section II). Finally, myocardial defects were not avai lable for this phantom, and their effect on motion estimation could not be studied. It is clear th at the significance of this paper re sts on the premise that in spite of the shortcomings of the physical phantom, the relative performance of motion estimation methods with this phantom data is e quivalent to that with patient data. Methods Data Acquisition A physical dynam ic cardiac phantom2 was modified for moti on tracking in this study. The phantom consists of a water-filled torso containing two lungs (each filled with two parts Styrofoam beads and one part water, by volume), a Teflon spine, and a beating heart (Fig. 3-1). The heart assembly is composed of a 50 mL rounde d cylindrical latex membrane inserted into an identically shaped 100 mL latex membrane. The volume between the two membranes, 2 Data Spectrum Corp., Hillsbourough, NC 27278-2300

PAGE 23

23 representing the myocardium, may be filled with radioactive solution to mimic, for example, a myocardial perfusion study with 99mTc sestamibi. A computer contro lled pump is used to expand and contract the volume within the inner membrane while simu ltaneously sending an ECG signal to the imaging device. Values for stroke volum e, ejection fraction, and beats per minute (BPM) can be input to the controlling computer. The point source marker data acquisition preceded the myocardial activity data acquisition. Ten small point source markers were c onstructed out of latex tubing and attached to the outer (epicardial) surface of the myocardium with approximately uniform spacing, as shown in Fig. 3-1. A dry syringe was inserted in one end of each marker while another syringe containing 99mTc was inserted in the other end. The dr y syringe was used to apply a negative pressure inside the marker, thus removing a ny air and ultimately pulling activity out of the opposite, hot syringe. A pproximately 3 mCi of 99mTc was injected into each marker. The volume of 99mTc solution in each marker was approximately 0.5 mL. The phantom was programmed for a 50% ejection fraction and 60 BPM. Three consecutive 60 min. gated SPECT scans we re acquired on a commercial SPECT system3 using low energy high resolution para llel-hole collimators over 120 angles thr ough 360 with 8 gated frames per heart cycle, 3.56 mm pixel size, a nd 15% energy windows. The three projection sets were then summed. Projection images from a singl e angle obtained from the marker acquisition are shown in Fig. 3-2. In order to ensure registration of the marker and myocardial activity images, the phantom was not moved while the markers were allowed to decay to negligible activity levels. Through an access port to the myocardial compartment, 99mTc solution was injected. Including activity in 3 Triad 88, Trionix Research Laboratory, Twinsburg, OH 44087

PAGE 24

24 other regions of the phantom w ould have required moving the pha ntom, and so this was avoided for this study. A high activity level (73 mCi) a nd a long scan time (8x60 min.) were used. This high count, low-noise data provided a means of generating an ensemb le of datasets at any count level using a Poisson noise generator. Ten noisy da tasets were generated su ch that the count total in a mid-ventricular slice was approximately 215,000. This count level is in a range between gated 99mTc sestamibi SPECT and 18F-FDG cardiac PET. The ten datasets provided a means of estimating the statistical uncertainty in the resu lts. Low noise, projecti on images obtained from the myocardial acquisition from a single angle are shown in Fig. 3-3. One of the ten noise realizations is shown in Fig. 3-4. Marker Motion Estimation The m arker data were reconstructed using 50 iterations of OSEM using 12 projection subsets. Noise in the reconstruction was minimal due to the high-count density of the acquisition; therefore, no spatial smoothing was applied. A 3D ROI was specified for each marker in each frame and the pixel intensity center-of-mass (C OM) within each was determined. The vectors describing the motion of the COMs between succe ssive frames were then considered the true motion of the myocardium at the location of eac h marker. Reconstructed images obtained from the marker acquisition are shown in Fig. 3-5. Myocardium Motion Estimation The m otion of the myocardium was estimated from the noisy projection data of the myocardial activity usin g two methods: (1) the optical flow method of Horn and Schunck [29] applied to optimized OSEM [30] reconstructi ons (OSEM-HS), and (2) a simultaneous image reconstruction and motion estim ation method [12] (RM). OSEM-HS Method. In order to ensure high quality images for the OSEM-HS method, the OSEM reconstruction was optimized, in a sum-of-squared errors (SSE) sense, in terms of the

PAGE 25

25 iteration stopping point, 3D spatia l filter cut-off frequency, and tem poral filter kernel. It should be noted that the parameters which optimize the SSE may not be the best parameters for motion estimation. Values were determined for these pa rameters that minimized the SSE between the scaled low noise, OSEM reconstruction and a re construction of one of the noisy datasets: n jr j jff SSE1 2))( )(() ,( rr ff (3-1) where ) ,...,,(21 nfff f and ) ,..., ( 21 nfff f are the scaled, low-noise reconstructed image and the noisy reconstructed image, respectively, over n frames, and r is the 3D voxel index. To reduce the influence of extra-myocardial noise on the optimization, the index r was confined to pixels in a region of the myocardium. This was achi eved using a region-of-interest defined as pixels with intensity above 10% the maximum intensity in the low-noise image. The region defined this way contained all of the myo cardium and an approximately one pixel border around the myocardium. The OSEM reconstruction of the noisy data wa s tested at stopping it erations of 5, 10, and 50. For each of the three stopping points, a 3D post reconstruction Butterworth filter was applied over a range of cut-off frequency values (0.42-1. 7 cycles/cm). Finally, a three-point temporal convolution filter was applied to each of the resulting filtered images with several different weights ({0.1,0.8,0.1} applied once and twi ce and {0.25,0.5,0.25} applied once). The SSE was computed for each image and the results are shown in Fig. 3-6. The smallest SSE (277.51) occurred using 5 iterations of OSEM with a Bu tterworth cut-off freque ncy of 0.91 cycles/cm and a single application of the temporal c onvolution with weights {0.1,0.8,0.1}. These optimal parameters were then also used to rec onstruct the other nine noise realizations.

PAGE 26

26 The motion was estimated by applying the Horn and Schunck optical fl ow algorithm [29] to the optimized OSEM reconstructions. The Ho rn and Schunck optical fl ow algorithm requires a user-defined parameter, that weights the contribution of the smoothness constraint. The optimal was chosen based on the accuracy of the estimated motion (using a Euclidean distance metric described in section II.D below) for a single noise realization, and this optimal was then used with the other ni ne noise realizations. RM Method. The RM method is described in detail in Mair et al. [12]; an overview is presented here. RM is a method for simultane ous reconstruction and motion estimation that alternates between pixel intensity upda tes (R-step) and motion updates (M-step). The R-step is a penalized likelihood method th at incorporates the estimated motion into the penalty term. The M-step employs sequent ial quadratic programming and the conjugate gradient algorithm to compute motion given the current image estimate. The M-step used the simpler smoothness constraint of Horn and Schunck [29] rather than the strain energy constraint in Mair et al. [12]. For th is study, 50 RM iterations were used where each RM iteration comprised 2 R-step iterations and 1 M-step iteration. Two user-defined parameters, and are used in RM. Parameter is the scalar for the likelihood function and was chosen empi rically for this study. Parameter weights the contribution of the smoothness c onstraint as in the HS method, and was chosen by the same method used in the OSEM-HS method. The RM reconstructions were filtered with a Butterworth filter with cutoff frequency equal to 1.5 cycles/cm. Example transaxial sli ce images of low noise OSEM, noisy OSEM, and noisy RM are shown in Fig. 3-7.

PAGE 27

27 2-Frame Motion Estimation To further evaluate the motion estim ation methods at higher magnitudes of motion, 2frame datasets containing only the end-dias tolic (ED) and end-systolic (ES) frames (corresponding to frames 8 and 3, respectively) from each noise realization were created. Optimal reconstruction and motion estimation pa rameters were determined using the same methods employed with the 8-frame data. Motion Estimation Error The accuracy of the estim ated motion vector fields was determined by calculating the Euclidean distance (i.e. the magnitude of the v ector difference) between the estimated and the marker-based, true motion vector in the vicinity of each marker. In computing this error metric, the estimated motion was the average motion over the eight voxels nearest the marker center-ofmass. The motion error was calculated for each marker in each frame and in each noise realization, for both the 8-frame and 2-frame data sets. The mean error and standard deviation ( ) over the 10 realizations was com puted at each marker location. Summing Estimated Motion Vectors Across Frames The estim ated motion vectors at each marker location were added across several systolic frame-to-frame intervals (8 2, 83) and compared to the net movement of the markers through corresponding frame intervals. The purpose for th is analysis was two-fold: (1) to determine whether the motion errors accumulate over frames or if one or both of the angle and magnitude components of the errors cancel across frames and, (2) to compar e the accuracy of the total systolic motion estimation using e ither a series of small frame intervals or a single large frame interval. This frame-to-frame vector summing i nvolved 3D linear interpol ation of the estimated motion vectors in the vicinity of the markers at each frame.

PAGE 28

28 Results 8-Frame Data The results of the optimal search f or OSEM-HS and RM for the 8-frame data are shown in Fig. 3-8. The figure shows the motion error as a function of averaged across all markers and all frames for a single noise realization. The optimal was 0.15 for both OSEM-HS and RM. This value of when applied to the ensemble data, resulted in an average motion error (including all markers and all frames) of 0.15 pixels ( = 0.001) for OSEM-HS and 0.13 pixels ( = 0.001) for RM. A more detailed presentation of these results is given in Fig. 3-9, which shows plots of the frame-by-frame, ensemble av erage motion error for four selected markers. The plots emphasize the relationship between moti on magnitude and estimated motion error (i.e. an increase in motion magnitude yields an incr ease in estimated motion error). Frames 1 and 4 are the frames of peak systolic and diastolic motion, respectively, and th us, the motion error is often greatest in these frames. The errors we re similar for both methods except frame 1, the frame of greatest motion magnitude, where the RM error was significantly smaller for 7 out of 10 markers. Table 3-1 summarizes the motion errors (ensemble average) at each frame by averaging across all markers. Also included is the true ma gnitude for each frame, again to emphasize the relationship between motion magnitude and motion estimation error. The significance of the differences in th e estimated motion error between the two methods was evaluated in each frame based on pair ed t-tests using the av erage marker error for each noise realization. The differences were significant at the 5% level for all but frame 7. In spite of the statistical signifi cance of the differences, the practi cal significance is questionable given the generally small magnitude of the differences.

PAGE 29

29 2-Frame Data The optim al for the 2-frame data was found to be 0.06 for OSEM-HS and 0.03 for RM. These values of applied to the ensemble data resulted in an average motion error including all markers and all frames of 0.55 pixels ( = 0.005) for OSEM-HS and 0.37 pixels ( = 0.004) for RM. Table 3-2 gives the frameby-frame motion error and true motion magnitude averaged across all markers. The differences between methods for both frames were found to be statistically significant at the 5% level. For this 2-frame data, which has greater frame-to-frame motion magnitude than the 8-frame data, the better performance of RM wa s more pronounced. The practical significance of this differe nce requires furthe r investigation. Summing Estimated Motion Vectors Across Frames The summe d estimated motion errors (average d over all ten marker s over all ten noise realizations) for both OSEM-HS and RM are given in Table 3-3. Since the angle error decreases as the motion vectors are summed over frames, whereas the magnitude error increases, these limited observations suggest that the angle error cancels over fr ames, and the magnitude error accrues over frames. The summed motion for the comp lete systolic range (Table 3-3, row 2) is only slightly worse than the motion obtained by estimating the motion directly between the two frames (Table 3-2, row 1). Further investigation is needed to determine if there are situations in which the method of summing over smaller intervals could produce improved accuracy over direct estimation of motion between two frames. Summary and Discussion This paper describes a study that used a phys ical phantom with point source markers to evaluate the accuracy of wall motion estimation for gated cardiac ET. Of the two estimation methods that were evaluatedthe established optical flow method of Horn and Schunck and a

PAGE 30

30 simultaneous reconstruction and motion estima tion methodthe simultaneous method resulted in a statistically significant lower overall error, and the difference was greatest for the largest true motion magnitude. The absolute difference between the methods overall was relatively small, and it is reasonable to qu estion the practical signifi cance of the difference if the methods were to be used, for example, within a motion-compensated reconstruction algorithm. Aside from assessing the relative pe rformance of these particular motion estimation methods, a more general contribution of this work has been the development of a physical phantom-based method for evaluating motion estimation methods. Several limitations of this study should be note d, and the first concerns the regions of the myocardium that were sampled for this study. The 10 point markers do not sample the entirety of the myocardium, and so the regional variation may have been missed. More specifically, the markers were located only on the epicardial, or ou ter, surface of the phant om myocardium. It is known that the motion magnitude on this surface is less than that on the endocardial, or inner, surface, and it is reasonable to que stion if the results would be the same for the more challenging endocardial region. Results from this study suggest that the relative advantage of the simultaneous method increases as the true mo tion magnitude increases. Finally, the markers could not be located within the myocardium (i.e. away from either surface), which would be expected to be more challenging than the high image contrast regions located along the surfaces. Another limitation of this study concerns the differences in motion characteristics of the phantom compared with the typi cal human heart. One prominent difference is the degree to which the true motion is tangential to the myocardi al surface. In human hear ts, this component is substantial and coincides with the wringing moti on that has been observed [31]; in the phantom,

PAGE 31

31 this component is relatively small. It is known that a weakness of the motion estimation methods tested here is their ability to estimate this tangential component of motion. The results of this physical phantom study are generally consistent w ith an earlier study that used a mathematical phantom [12]. Specifi cally, both studies demonstrated more accurate motion estimation with the simultaneous RM me thod compared with an independent motion estimation method applied to fixed, reconstructed images. In the earlier work, the independent estimation method was the M-step of the RM me thod applied to optimized OSEM images. The earlier work also demonstrated improved im age quality with RM compared with OSEM; however, the study here showed only minor di fferences in image quality with the two reconstruction methods. In terms of and different values were found appropriate for this data compared with the mathematical phantom data. For this difference is due to the different smoothness constraints used in the two studies. For this difference is due to the extramyocardial (i.e. liver) activity included in the mathematical phantom, which greatly affects the likelihood function value, and, therefore, the choice of

PAGE 32

32 Figure 3-1. DCP on imaging bed with close-up of radioactive markers shown on right. Outer shell of phantom is removed. Figure 3-2. Projection images obtained from a single angle of the marker acquisition. Figure 3-3. Low noise projecti on images obtained from a singl e angle of the myocardium acquisition. Figure 3-4. Single noise realization obtained after scaling an d adding Poisson noise to the projection images in Fig. 3-3.

PAGE 33

33 Figure 3-5. Reconstructed images obtained from the marker acquisition, shown with slices 1434 summed. Figure 3-6. SSE vs. Cut-off frequency for 5 (l eft), 10 (middle), and 50 (right) iterations of OSEM, where 4-D OSEM 1 = single applic ation of {0.1,0.8,0.1} temporal filter, 4-D OSEM 2 = two applications of {0.1,0.8,0.1} te mporal filter, and 4-D OSEM 3 = single application of {0.25,0.5,0.25} temporal filter. Figure 3-7. Transaxial slices of frame 1 of the low noise OSEM, noisy OSEM, and noisy RM reconstructions. Profiles for the oblique line shown in the low noise image are normalized and shown in the plot.

PAGE 34

34 Figure 3-8. Euclidean distance vs. Beta for OSEM-HS (left) and RM (right). Figure 3-9. Euclidean distance (ensemble average) vs. frame number for markers 2, 4, 6, and 10. The plots emphasize the relationship between motion magnitude (greatest in frames 1 and 4) and estimated motion error.

PAGE 35

35 Table 3-1. Average estimated motion errors for 8-frame data (pixel units) Frame number True magnitude OSEM-HS Error RM Error 1 0.939 0.429 0.328 2 0.286 0.100 0.087 3 0.217 0.077 0.069 4 0.463 0.150 0.141 5 0.307 0.119 0.116 6 0.194 0.102 0.101 7 0.211 0.110 0.105 8 0.166 0.098 0.100 Table 3-2. Average estimated motion errors for two-frame data (pixel units) Frame number True mag. OSEM-HS Error RM Error 1 (ED ES) 1.347 0.462 0.351 2 (ES ED) 1.347 0.642 0.386 Table 3-3. Summed motion errors for OSEM-HS and RM OSEM-HS RM Frame interval Motion error (pixels) Angle error (degrees) Mag. error (pixels) Motion error (pixels) Angle error (degrees) Mag. error (pixels) 8 2 0.441 17.6 -0.3680.34715.6-0.252 8 3 0.487 14.5 -0.4010.38613.0-0.288

PAGE 36

36 CHAPTER 4 IMPROVED RESPIRATORY MOTION CORREC TION IN GAT ED CARDIAC SPECT Introduction In cardiac im aging modalities requiring scan times longer than a si ngle breath-hold (e.g. PET, SPECT, MR), the motion of the heart due to respiration introduces blur into the tomographic images [14]-[16], [33] -[34]. By sensing the phase of the respiratory cycle, a gating procedure may be used to produce a time series of images with reduced motion blur. When combined with electrocardiographic (ECG) ga ting of the cardiac contraction cycle, the acquisition is referred to as a dual-gated study [35] -[36] and produces a matrix of projection sets. A drawback to respiratory-gating, similar to EC G-gating, is higher noise levels in the gated images [37]-[38]. To combat the increased noise levels of respiratory-gated images in emission tomography, motion correction using an estimate of the heart motion between the image frames combined with temporal summing, has been proposed [17]-[18 ],[36],[39]-[46]. Gating of the respiratory cycle may be accomplished by affixing a piezo-electric elasticized belt ar ound the patients chest. As the pati ent breathes, a signal is generated corresponding to the tension in the belt. The si gnal may be input to the imaging device as a representation of the respiratory cycle, similar to ECG-gating of the cardiac contraction cycle. Image data may be binned in terms of respiratory amplitude or phase, where amplitude binning has less in-frame motion blur th an phase binning, but phase binning has bette r noise level uniformity across bins [47]. Given the gated image frames, many technique s currently exist for estimating respirationinduced heart motion. Physical devices have been used that track the motion of markers placed on the patients chest throughout the acquisition [17],[39]-[40] These devices offer high precision; however, the motion of deeply-embedde d structures within th e body may be different

PAGE 37

37 than that of the surface of the body. Furthermore, implementation of phys ical tracking methods requires hardware integration and calibration fo r each scanner in a facility, which may be expensive and require regular quality assurance. Motion estimation methods based on image data analysis alone are easier to implement and have the advantage of estimating the motion of the organ directly. Rigid-body [18],[41]-[43], affine [44],[45], and deformable models [16],[36] of the respiratory heart motion have all been shown to successfully decrease respiration-indu ced artifacts in cardiac imaging. In SPECT, the rigid-body model of the hear t [48] has been used widely due to spatial resolution constraints, which place an upper-bound on the extent of the improvement possible using motion estimation and temporal summing. A drawback of image data analysis methods is th at their accuracy is highly influenced by image noise, reconstruction artifacts, and available computational resources. Within the data analysis approach, both analytical and iterative methods have been investigated in estimating the respiratory he art motion from tomographic data. Analytical methods [18],[41]-[43] are genera lly faster than iterative methods [44]-[46], but their accuracy may be highly dependent on the level of noise in the images. Iterative methods, aside from their high computational costs, often require a good initial estimate to reach convergence [49]. In Klein et al. [44], the affine moti on for a given frame-pair used an initial estimate defined as the motion estimate from the immediately preceding frame-pair. A cost function was implemented which penalized deviations from this starti ng estimate. The method showed convergence, but required a user-defined weighting parameter on the cost function, which may complicate the use of the method in a clinical environment. Furthermore, the author s reported an average computation time of 720s. The method in Dey et al. [45] used spatio-temporal penalty terms

PAGE 38

38 which penalized deviations from the simple-harmonic oscillation equation. The method demonstrated convergence but again required user-d efined weighting parameters on the penalty term. We previously investigated methods for es timating and correcting the respiratory motion of the heart iteratively without user-defined parameters or pena lty terms [46]. The method used a rigid-body model with rotation pa rameterized by Euler angles and iteratively minimized an image-registration function using an optimized conjugate gradient method. The Euler angles correspond to a rotation about th e x-axis followed by a rotation a bout the y-axis, followed by a second rotation about the x-axis. In Parker et al. [46], the mo tion was estimated after summing the ECG-gated frames. This summing adds motion blur to the images, but increases the SNR. The method was slow due to repeated calculati on of the objective functi on and its derivatives which required extensive tri gonometric evaluations. On modern computer architectures, trigonometric functions are approximated using si mple iterative equations [50], and thus may be ill-suited for use in deeply-embedded loops. To improve the computational efficiency of the method in [46] while maintaining high accuracy and convergence, here we investigate the use of the quaterni on parameterization of rotation [51] in iteratively estimating the rigid body respiratory motion of the heart. Quaternions have been used extensively in computer anima tion and gaming industries due to their simplicity and improved physical characteristics [52] comp ared to Euler representations. For example, gimbal lock, a condition resulting when the second Euler angle is equal to 90 and a degree of freedom is lost, is avoided using the quatern ion representation. Furt hermore, a quaternionparameterized rotation matrix is free of trigon ometric functions, and we hypothesize that this simplification will reduce the computational burden of repeated calculations of the objective

PAGE 39

39 function and its partial derivativ es. We compare this method quantitatively in terms of speed, accuracy, and robustness to several other iterative and analytical techniques, using simulated data of a mathematical phantom. The quaternion met hod is further evaluated on a respiratory phase binned human SPECT image in terms of computat ion time and image quality after correction by the estimated motion. It is shown that the use of the quaternion representation enables a fast, accurate, and robust model for respiratory motion co rrection, which has practical applications not only in SPECT but also in other imaging m odalities in which breath-hold acquisition is impractical, such as PET and MR. Theory We represen t the image space by a discrete space variable r = ( x,y,z ), and assume a series of tomographic image frames 12,,, J f ff We assume that the frame-to-frame motion of an object can be modeled as a rigid body. We cons ider the problem then of forming a single composite image by registering and summing all frames to the reference frame, 1 f using an estimate of the rigid-body motion between each frame and the reference frame. The rigid body motion between 1 f and j f is estimated by minimizing the image registration objective function: 2 1Q, Q+jjEff rbrrb (4-1) where b is the 3D translation motion, j f is the image at frame j r is the discrete 3D spatial coordinate for each voxel center, and Q is the 3D rotation matrix. The sum is taken over all voxels. A total of 6 parameters (3 translati onal, 3 rotational) are needed to compute j E In contrast to other iterative methods no user-specified parameters or penalty terms are used. Since

PAGE 40

40 the rigid body motion is a smooth vector fiel d which is completely specified by few (6) parameters, and the images contai n orders of magnitude more voxe ls that are being matched, we do not see the need to include any additiona l penalty terms in the objective function. The rotation matrix Q may be parameterized by twelve different se ts of three Euler angles, where each angle represents a rotation ab out a Cartesian axis. In each case, specific rotation matrices have multiple distinct parame terizations. This non-uniqueness is due to the dependence of these representations on trigonomet ric functions. This motivated us to use the quaternion representation which has a unique representation, and onl y involves simpler algebraic functions [51]. For the readers convenience we give a brief outline of quaternion-parameterized rotation. See Horn [51] for a more complete treatment. A quaternion q may be viewed as the sum of a scalar and a 3D vector as follows 0123 qijk (4-2) where ,, ijk are unit vectors in 3-space and 0123,,, are scalars. q is a unit quaternion if 2222 01231 (4-3) The unit quaternion cos(/2)sin(/2)qn represents the rota tion through the angle about the axis parallel to the unit vector n (with initial point the origin), where the direction of rotation and n form a right-handed system (i.e. is the direction of ro tation of a right-handed screw being tightened with the tip moving in direction of n). Thus, we impose the condition 00 which guarantees the unique quaternion para meterization for rotation matrices. So, from equation (4-3), we have 222 01231 Thus, by using this representation, each

PAGE 41

41 rotation matrix is determined by the unit quaternion vector 123 ,,q which is explicitly stated in Appendix B, equation (4-B.1). Thus we minimize the objective function in equation (41) with respect to the variables 123123,,,,,bbb where 123,,bbb b. This minimization was performed by using th e modified conjugate gradient algorithm CG_DESCENT developed by Hager and Zhang [53]-[54]. This is a globally convergent nonlinear conjugate gradient met hod with guaranteed descent and a fast and highly accurate line search technique. As a result, we refer to the proposed method as the CGQ (conjugate gradient quaternion) method. Experimental Methods The proposed motion estimation method was co mpared with three existing methods on two different datasets using both quantitative metrics and visual inspection. In this section we discuss details of implementation of the mo tion estimation algorithms, how these motion estimates were used to determine a single summed image, the generation of the simulated data and a description of the patient data, and the evaluation methods for the motion obtained from all four methods applied to the two datasets. Algorithms We compared our algorithm (CGQ) with the principal axes method [55],[56], the generalized centers of mass method [42], and one which uses the same CG_DESCENT algorithm to minimize the objective function in equation (4-1) but re presents the rotation matrix by Euler angles [46] instead of quaterni ons. Each algorithm was used to minimize j E in equation (4-1) for j = 2 to j = J, where J is the total number of image frames. The partial derivatives of (1) used in th e CGQ algorithm are given in Appendix C. The CGQ algorithm was initialized with each component of q set to zero, and each component of b

PAGE 42

42 found by calculating the 3D center of mass in each frame. In the rare case that 0 became imaginary, a very large ne gative number (approaching ) was returned for each rotational component of the gradient. The conjugate gr adient algorithm (CG_DESCENT) was stopped when the change in the objective function betw een successive iterations was less than 0.l%. For comparison to the quaternion method, we implemented three methods that have been used previously for rigid-body motion estimation in medical imaging. These include an iterative conjugate gradient approach a nd two analytical approaches: pr incipal axes and generalized centers-of-mass. Instead of using quaternions, methods ha ve been developed based on Euler angle representations of the rotation matrix Q in equation (4-1) [44]-[46]. Here we compare our method with one which uses th e pitch-yaw-ro ll system ( ,, ) corresponding to rotations about the x-axis, the y-axis, a nd the z-axis, in that order. The rotation matrix and partial derivatives for this system are given in Appe ndix A. The algorithm was initialized with each rotational component ( ,, ) set to zero, and each component of b found by calculating the 3D center of mass in each frame. The stopping rule was id entical to that used in CGQ. We refer to this as the CGPYR (conjugate grad ient pitch-yaw-roll) algorithm. In the principal axes method the translationa l motion is first estimated by the difference in the centers of mass (COMs) between the refe rence and deformed frames. Then the deformed frame is translated so that both frames have the same center of mass. The rotational motion is then estimated by the rotation that aligns the pr incipal axes of both frames. This rotation is estimated from analytic relationships between the singular value decompositions (SVDs) of the inertia matrices for the reference and translated deformed frame. This method is not based on iterative algorithms for minimizing an objective function as the CGQ and CGPYR methods, so is

PAGE 43

43 much simpler to implement, but there is no comp ensation for noise in reconstructed images. The accuracy and speed of this algorithm are very de pendent on how the SVDs are computed. In this paper, we used the SVD algorithm of the LAPA CK [57] software library (driver routine DSYEVD). The CGQ algorithm is also compared with the generalized center-of-mass (GCOM) algorithm developed by Feng et al. [42]. By considering the intensity map as a distribution of masses, the center of mass of an image can be computed in te rms of the moments of these masses about the coordinate axes. The GCOM me thod is based on the concept of generalized center of mass points which are defined by replac ing the (order 1) moment in the usual COM formula, by moments of higher order, 1,2,3, For any triplet of orders, the resulting generalized center of mass points are used to estimate the motion by a least squares fitting algorithm. The motion estimates obtained from a fi nite number of different triplets are then compared using the root-mean-squared distance (RMSD) between the reference and deformed (by the motion estimate) frames. The motion that minimizes this distance (from the finitely many possible motion estimates) is designated to be the motion estimate. Here, we used the code provided by Dr. Bing Feng which is based on the SVD routine take n from the Numerical Recipes software library [58]. Mathematical Phantom Evaluation Image generation: Tests were performed on an ensemble of 10 noisy realizations of a dual gated (respiratory and ECG binning) ca rdiac SPECT scan generated by the 4D NURBSbased Cardiac-Torso (NCAT) phantom [59]. Rela tive organ activity concentrations were based on rest 99mTc sestamibi studies taken from Gilland et al. [60] and are given in Table I.

PAGE 44

44 To generate realistic moti on blur within the ECG gated frames, ECG gating was simulated with 32 frames for one cardiac contracti on cycle, and then averaged to a final sampling of 8 frames. Respiratory motion was simulated by rotating and shifting each of the ECG-frames for a total of 32 respiratory-frames per ECG gate These frames were then averaged to a final sampling of 8 respiratory fram es per ECG gate. These transf ormations resulted in an 8x8 (respiratory cycle bins x ECG bins) matrix of phantom images. The motion used to generate the respiratory frames is given in Table II. These values were taken from human MR studies in McLeish et al. [16] and represen t a slightly larger than averag e displacement. In the table, bx represents the lateral translation, by represents the anteriorposterior translation, bz represents the cranio-caudal translation, represents the rotation about the lateral axis, represents the rotation about the anteriorposterior axis, and represents the rotation about the cranio-caudal axis. This notation is identical to that used in CGPYR. Non-rigid de formation of the heart due to respiration was not modeled. Projection data fo r each of the dual-gated phantom distributions were simulated using the SIMI ND Monte Carlo program [61]. S catter, attenuation, detector stopping power, and detector re sponse were modeled assuming a gamma energy of 140 keV, crystal thickness of 1.27 cm, pixel size of 0. 3125 cm, 64 projections over 180 from the leftposterior-oblique to right-an terior-oblique angles, 128 de tector elements, 15% energy resolution, 20% energy window, and a low energy ge neral purpose parallel-hole collimator. The simulations were allowed to run until noise in the projection data was negligible ( 2%). The data were then scaled to approximately 14,000 total counts in a 3 mm mi d-ventricular slice per frame. This count level is typical of a 99mTc-sestimibi study. Finally, Poisson noise was simulated in each projection set a total of 10 times generating an ensemble of noise realizations. Each of the noise realizations was reconstructed using 5 it erations of ordered-subsets

PAGE 45

45 expectation-maximization (OSEM) [13] with 8 s ubsets without attenuation correction. The dualgated reconstructed images were summed over th e ECG frames to generate the 8 respiratorygated frames upon which the motion estimation methods were tested. Evaluation metrics: For any motion estimate determined by the rotation matrix Q and the translation vector b we define the phantom matching error (PME): 2 1 2QQ +J PPjjj jPME ff,, r,b rrb (4-4) where Pj f is the original phantom image without extra-myocar dial activity, at frame j The rotational error is defined as the absolute di fference between the estimated rotation and true rotation, in the Euler representation, in degree units. The translational error is defined as the absolute difference between the estimated transl ation and true translation, in pixel units. Each of the motion estimation methods was tested on the original, noise-free phantom images without extra-myocardial activity, and on the noisy reconstructed images. Their accuracy was quantified in terms of translational, rotati onal, and phantom matching errors, each averaged over all eight frames. In the case of the motion es timated from the original phantom images, the PME is a measure of the optimal registration accuracy of each method. The total computation times over all frames were also compared. The noise-free phantom frames were correct ed by translating and rotating frames 2 through 8 based on the estimated motion, and summi ng with frame 1. The resulting images were evaluated by visual inspection a nd profile analysis. For comparison purposes, three other images were generated: 1) an uncorrect ed image, i.e. a summation over the respiratory frames without any motion estimation, 2) an image corrected with the known, true motion, and 3) an ideal image

PAGE 46

46 without any respiration-induced motion (as if the patient held their breath for the duration of the acquisition). Effect of segmentation: Estimating the respiratory motion of the heart on images with extra-myocardial activity requires segmentation of the myocardial activity due to the nonuniform displacement of organs in the thorax du ring respiration. To inve stigate the effects of segmentation on the motion estimation methods, we segmented the left ventricular (LV) wall of the noisy reconstructions over a range of segmentation levels. Th e different segmentations were created by first varying a 3D Butterworth filter cutoff frequency from 0.32 to 0.64 cycles/cm with a step size of 0.04 cycles/cm. Then, an inte nsity threshold was varied from 0 to 25% of the maximum pixel intensity with a step size of 2.5%. Pixel values below this threshold were set to zero. This resulted in 99 different segmentations for each noise realization (990 total segmented images). Fig. 4-1 shows 16 of the 99 segmentations for the first noise realization. We can see that as the threshold is increased, less of the extra-my ocardial activity is included. The optimal level of segmentation for each method was defined as th at which generated the lowest average PME. Comparison of methods at optimal segmentation: The methods were then compared at the optimal segmentation level for each. The errors were computed as averages over all 10 noise realizations and all 8 respirat ory frames. A single noise reali zation was corrected using the estimated motion from each method, followed by a 3D Butterworth filte r of cutoff frequency 0.48 cycles/cm. This level of smoothing was chosen empirically to be comparable to that which we may expect in a clinical sett ing. The resulting images were evaluated by visual inspection and profile analysis. Again for comparison, an uncorrected, true-motion corr ected, and ideal image were also generated.

PAGE 47

47 Effects of varying levels of extra-myocardial activity: In a clinical setting, we would expect to see a high variability in the level of extra-myocardial activity present among patients. The maximum standard deviation of extra-myocardi al activity in Gilland et al. [60] was 75%. To test the robustness of the methods under these variable conditions, we calculated the average PME over ten noise realizations for each me thod at optimal segmentation as the extramyocardial activity was varied fr om 0 to 300% (2 standard devi ations from average). A single noise realization for 25% and 175% background activ ity scaling levels was corrected using the motion estimation methods, followed by a 3D Bu tterworth filter of cutoff frequency 0.48 cycles/cm. Respiratory-Gated Patient SPECT Images The methods were further evaluated with a respiratory-gated ac quisition of a human subject. The patient was imaged at rest after injection of 99mTc-sestimibi with a matrix size of 128, pixel size of 0.467 cm, with 68 angles ove r 204 about the left-ant erior-oblique angle. Gating of the respiratory cycle wa s performed using a chest belt as described in section 4.1, with 8 gates per respiratory cycle binned post-acquisi tion according to phase. The patient dataset was stripped of patient identifiers in complianc e with the Health Insu rance Portability and Accountability Act (HIPPA)1. The data were first reconstructed using 5 itera tions of OSEM with 4 subsets. To estimate the respiratory motion, a segmentation of the LV wall was created by first smoothing with a 3D Butterworth filter of cutoff frequency 0.4 cycles/cm, followed by an intensity threshold of 27% of the maximum pixel intensity. A user specifie d 3D elliptical region-of-interest (ROI) in each frame was drawn around the LV, and all voxels outs ide this ROI were set to zero. The images 1 IRB Protocol #157-2004.

PAGE 48

48 were then cropped to a final matrix size of 24 with 24 slices. These parameters were chosen empirically, and a transaxial slic e of the segmentation at end-e xpiration is shown in Fig.4-2. The motion estimation methods were applied to the data, and the resulting motion estimates were used to correct the original, unf iltered reconstruction. This was followed by a 3D Butterworth filter of cut-off frequency 0.37 cycles/cm. This level of smoothing was chosen empirically to be comparable to that which we ma y expect in a clinical setting. For comparison, an uncorrected image was also created by summi ng the original, unfiltered reconstructions across the respiratory frames, and smoothing with an identical 3D Butterworth filter of cut-off frequency 0.37 cycles/cm. Since the true moti on is not known in these images, the resulting images were evaluated only by visual inspection of myocardial wall blur and uniformity, average objective function value over all frames, and computation time. Results All motion estimation methods were implemen ted on a standard Linux workstation with dual AMD Opteron 250 (2.4 GHz) processors and 4 Gb of memory per processor. In testing, the methods were restricted to ex ecute on a single processor. Motion Estimation on the Noise-Free Phantom Frames Results for the motion estimates on the noise-fre e phantom frames are given in Table III. These results indicate that in terms of accuracy all correction methods have similar registration and translational errors, but the ro tational error for CGQ a nd principal axes are better than all the other methods. The principal axes method is the computationally fastest followed by GCOM, CGQ, and CGPYR. Note that using the same CG algorithm, but only changing the parameterization of the rotation matrix from the Euler angle representati ons used in CGPYR and our previous work [46], resulted in a 50% reduction in the computation time.

PAGE 49

49 Fig. 4-3 shows selected short and horizont al views of the noise -free phantom without respiratory motion (ideal), with respiratory motion corrected e ither by the true motion, CGQestimated, or principle axes estimated motion, and with respiratory mo tion with no correction. CGQ and principal axes corrected images are displayed here because they represent the fastest of the iterative and analy tical methods, respectively. The imag es are shown at the end-diastolic phase of the cardiac contraction motion. The imag es corrected with CGQ and principal axes are more similar to the image corr ected with the true motion in terms of uniformity and sharpness throughout the LV and RV walls, than the uncorrected image. We can see that none of the corrected images achieves the le vel of myocardial uniformity a nd sharpness of the ideal image. This may be attributed to the presence of w ithin-frame motion blur as well as interpolation effects inherent in th e registration process. Profiles for the lines shown in Fig. 4-3 are gi ven in Fig. 4-4 for the ideal, CGQ-corrected, and uncorrected images. Profiles for the images co rrected with the true motion and principal axes motion were nearly identical to the CGQ profile and thus are not displaye d here. The plots show that the corrected images are mo re similar to the ideal images than the uncorrected images. The corrected profiles are ta ller, indicating better c ontrast, and thinner, re presenting better spatial resolution, than the uncorrected image profiles. Based on the similar accuracy but significantly worse computation time of CGPYR and GCOM to CGQ and principal axes, respectively, these methods were excluded from any further evaluation. Effect of Segmentation The results of the segmentation investigation on the reconstructions of the simulated data are shown in Fig. 4-5. The figure shows for the principal axes and CGQ methods a plot of the average PME as a function of intensity threshold and 3D Butterworth cu toff frequency. We can see that CGQ generates a relatively flat and lo w surface, indicating accurate motion estimation

PAGE 50

50 across a wide range of segmentation levels. The optimal segmentation for CGQ was at threshold = 17.5% and cutoff freq. = 0.44 cycles/cm, genera ting an average PME of 343.26. Principal axes was highly unstable across the range of segmenta tions, with optimal segmentation at threshold = 12.5% and cutoff freq. = 0.48 cycles/cm, genera ting an average PME of 757.23. The differences between the two plots are most profound for th resholds greater than 15%. These levels of segmentation correspond to the top two rows of Fig. 4-1. Comparison of Methods at Optimal Segmentation The average motion estimation results fo r principal axes and CGQ at optimal segmentation over all 10 noise re alizations are given in Table IV. CGQ performed significantly better than principal axes in te rms of PME accuracy and redundancy (standard deviation) with an average PME of 343.26 100.52. However, the average computation time of CGQ was 55.10 11.13s. The average PME for principal axes was 757.23 158.38, but it again executed on average in under 1s. Fig. 4-6 shows selected short and horizont al views of the noise -free phantom without respiratory motion (ideal), with respiratory motion corrected e ither by the true motion, CGQestimated, or principle axes estimated motion, and with respiratory mo tion with no correction. The images are shown at the end-diastolic phase of the cardiac contraction motion. The images corrected with CGQ and principal axes are more similar to the image corrected with the true motion in terms of sharpness through out the LV walls as well as LV blood pool contrast, than the uncorrected image. Profiles for th e lines shown in Fig. 4-6 are gi ven in Fig. 4-7 for the ideal, CGQ-corrected, and uncorrected images. Profiles fo r the images corrected with the true motion and principal axes motion were again nearly id entical to the CGQ profile and thus are not displayed here. The profiles for the corrected im ages are again more similar to the profiles for

PAGE 51

51 the ideal images than the uncorrected images. Sp ecifically, the corrected image profiles are taller than the uncorrected image prof iles, indicating better contrast. Effect of Varying Level of Extra-Myocardial Activity In a clinical setting we would expect a large variation in inter-patient extra-myocardial activity levels. Therefore, he re we examine how the motion estimation methods behave (at optimal segmentation) as the extra-myocardial activity level is varied. The results for principal axes and CGQ as the background activity levels ar e varied from 0% to 300% are shown in Fig. 4-8. We can see that the accuracy of CGQ improves as the level of background activity decreases, and worsens as the background levels are in creased. However, even as the background activity is increased to 300%, the CGQ-generated motion estimate is still more accurate than a motion estimate of zero for each component (dashed line at 4447.94). This shows good robustness on the part of CGQ over a large range of background activities, as we may expect to see in a clinical setting. In contrast, the princi pal axes method is highly unstable across the range of background activity levels. Images obtained using principal axes and CGQ at background activity scaling factors of 0.25 and 1.75 ( 1 stan dard deviation) are shown in Fig. 4-9. The CGQ-corrected images are more similar to the imag es corrected with the true motion in terms of LV wall uniformity and definition, than the uncorrected images. The images corrected with the principal axes estimated motion are less similar to the images corrected with the true motion than the uncorrected images. An examination of the principal axes motion es timates found that for several frames at each bac kground activity scaling level, the method generated rotational estimates with errors greater than 90. Human SPECT Image The CGQ-corrected, principal axes corrected and uncorrected human images are shown in Fig. 4-10. We can see that the CGQ-corrected images are more uniform and have decreased

PAGE 52

52 blurring in areas running perpendicular to the ax ial direction (marked with arrows) compared to the uncorrected image. The principal axes co rrected image does not resemble the expected activity distribution. Examination of the motion estimation obtained from principal axes found that two frames had rotational estimates greater than 90. The average objective function at convergence for CGQ was 431.16. An average ob jective function value calculated using the principal axes estimated motion was 1368.48. Thes e may be compared to an average objective function of 631.32 calculated using a motion estimate of zero for all components (i.e. no correction). Contrary to the previous studies using the phantom and noisy reconstructions, the computational time for CGQ was 6s. This may be at tributed to the larger pixel size of the human study that enabled a small segmented ROI (24x24x24) to be used for the motion estimation. The computation time for principle axes was less than 1s. Summary and Discussion In this work, we have developed a new method for respiration-induced cardiac motion correction using a rigid body model with a rotati on matrix parameterized by a unit quaternion. The method minimizes an image registration func tion using an optimized conjugate gradient routine. The implementation uses no user-defined input parameters or prior terms, simplifying the use of the method in a clin ical setting. Images corrected with the CGQ-generated motion were found to be more similar to images correct ed with the known, true motion compared to uncorrected images on both phantom and simula ted data. Our method was also tested on a respiratory-gated study of a human subject. The corrected images have increased uniformity and decreased motion blur in areas of the myocardium running perpendicular to the axial direction. Furthermore, the method was shown to be relatively insensitive to changes in segmentation and extra-myocardial activity level, a key requirement for clinical use. The computation time ranged from 61s on a phantom image to 55.1s on reconstr uctions of simulated da ta to 6s on a human

PAGE 53

53 SPECT image. More human studies over a range of detector configurations are needed to determine a typical computation time in routine clin ical use, but all of the times reported here are within an acceptable range, especially when ha rdware acceleration techniques are considered. It should be noted that although consistent improveme nt between corrected and uncorrected images was found, the difference was not always significant, and it is unclear at this time how respiratory motion correction may impact diagno stic accuracy in cardiac SPECT. Future work will focus on a mathematical observer study as we ll as a human receiver operating characteristic study on the effects of correcting cardiac respiratory motion blur in cardiac SPECT.

PAGE 54

54 Figure 4-1. Select segmentati on levels. The images are summations over the cardiac cycle and shown at end-expiration. Thre shold is given in units of percent of maximum pixel intensity, cutoff frequency is gi ven in units of cycles/cm. Figure 4-2. Transaxial sli ce of the segmented human image used for motion estimation (shown at end-expiration).

PAGE 55

55 Figure 4-3. Select slices of th e phantom images at end-diastole. Figure 4-4. Horizontal (l eft) and short (right) axis profiles for Figure 4-3.

PAGE 56

56 Figure 4-5. Results for average PME found by va rying the segmentation parameters. Threshold is given in percentage of maximum pixel intensity, cutoff frequency is given in cycles/cm. Figure 4-6. Reconstructions of the noi sy simulated images at end-diastole

PAGE 57

57 Figure 4-7. Horizontal (l eft) and short (right) axis profiles for Fig. 4-6. Figure 4-8. Results for average PME found by va rying the background activity concentration at optimal segmentation. Error bars indicate th e standard deviation of the PME over the ten noise realizations, for each scale factor. The dashed line at 4447.94 represents the PME of an uncorrected image.

PAGE 58

58 Figure 4-9. Reconstructions of the simulated data at extra-myocardial ac tivity levels of 0.25 and 1.75 (fraction of average). C GQ-corrected images are near ly indistinguishable from images corrected with the true motion. Figure 4-10. Reconstructions of the human images with CGQ correction, principal axes correction, and no correction, respectively. The CGQ-corrected image has increased uniformity and sharpness in areas perpendicular to the axial motion (marked with arrows) compared to the uncorrected image.

PAGE 59

59 Table 4-1. Relative organ activity concen trations for the simulated NCAT phantom Organ Relative activity Myocardium 75 Heart blood pool 6 Liver 13 Gall bladder 324 Lung 6 Kidney 45 Spleen 45 Bowel 37 Background (body) 6 Table 4-2. Respiratory motion in the mathema tical phantom images relative to frame 1. Translations are given in pi xels, rotations in degrees. Frame interval bx b y bz 1 1 0.0 0.0 0.0 0.0 0.0 0.0 1 2 -0.28 -0.864 -0.96 -1.675 -1.65 -0.375 1 3 -0.56 -1.728 -1.92 -3.35 -3.3 -0.75 1 4 -0.84 -2.592 -2.88 -5.025 -4.95 -1.125 1 5 -1.12 -3.456 -3.84 -6.7 -6.6 -1.5 1 6 -0.84 -2.592 -2.88 -5.025 -4.95 -1.125 1 7 -0.56 -1.728 -1.92 -3.35 -3.3 -0.75 1 8 -0.28 -0.864 -0.96 -1.675 -1.65 -0.375 Table 4-3. Motion estimates from the noise-free phantom images Method Translational error (pixels) Rotational error (degrees) PME Computation time (s) CGQ 0.15 0.01 100.0361.00 CGPYR 0.15 0.15 100.03119.72 Principal axes 0.15 0.01 100.810.13 GCOM 0.14 0.36 107.9654.59 No correction 1.4 2.5 4447.94Table 4-4. Average motion estimation errors fo r the 10 noisy simulated images with background activity at optimal segmentation Method Translational error (pixels) Rotational error (degrees) PME Computation time (s) CGQ 0.24 .05 1.05 .16 343.26 100.52 55.10 11.13 Principal axes 0.31 .07 2.25 .42 757.23 158.38 0.08 .04 No correction 1.4 2.5 4447.94

PAGE 60

60 CHAPTER 5 SIMULTANEOUS IMAGE RECONSTRUCTION AND M OTION ESTIMATION IN GATED CARDIAC EMISSION TOMOGRAPHY: A MATHEMATICAL OBSERVER STUDY Introduction In this chapter, we evaluate the image quality of the simultaneous image reconstruction and motion estimation method compared to a st andard reconstruction method in terms of a binary lesion detection task in SPECT myo cardial perfusion imag ing. ECG-gated source distributions of a normal human torso and a myocardial perfusion defect were simulated using a mathematical cardiac torso phantom. Projection data were obtained from the source distributions using a Monte Carlo code which modeled the ef fects of scatter, attenuation, and detector response. A total of 500 noise re alizations for each projecti on set (500 normal, 500 abnormal) were reconstructed using the simultaneous im age reconstruction and motion estimation method and a post-smoothed ordered-subs ets expectation-maximization ( PS-OSEM) [13] reconstruction. A channelized Hotelling observer [2 1],[62] was formulated to operate on single short axis slices (CHO), as well as on the complete multi -frame images (MACHO) [63]-[64]. During the study we noticed an inconsistenc y between the CHO and human observers in that the optimal level of smoothing for the CHO corresponded to Butterwor th cutoff frequencies less than 0.05 cycles/pixel. Because we felt that an evaluation of the methods at such high levels of smoothing would not necessarily predict the relative performance of the methods in a clinical setting, we chose to adjust th e reconstruction parame ters for RM and PS-OSEM such that the methods produced images with spat ial resolution characteristics t ypical of clinical images. The CHO was then used to generate an ensemble of test statistics for each reconstruction method, and a receiver operating characteristics methodology was used to anal yze the statistical ensembles. The area under the ROC curve (AUC) was used as an image quality figure of merit, and the statistical significance of the difference be tween the ROC estimates was determined.

PAGE 61

61 The signal-to-noise ratio (SNR) for each met hod was also calculated directly using a regions-of-interest technique across the ensemb le data. The statistical significance of the difference in SNR between RM and PS-OSEM using the regions-of-interest method was evaluated using a paired t-test. Background Motion Compensated Image Reconstruction It has been shown that incl uding non-rigid motion estimation in the penalty term of a maximum a posteriori (MAP) reconstruction can im prove the noise properties of reconstructed gated images while retaining a low level of motion blur [9],[10],[12],[24],[ 65]. The idea was first investigated in Lalush et al. [9] where the auth ors modeled the motion of a mathematical beating heart as an ellipsoid undergoing affine transfor mations. The process involved two steps: 1) motion estimation based on general known properties of the human heart followed by, 2) a reconstruction of the gated images using the estim ated motion vectors as Gibbs priors in a MAP reconstruction. The authors found that clique weights based on the estimated motion vectors generated by the affine motion model enforced better spatial resolution and comparable noise levels compared to several other no-motion families of Gibbs priors. In Gravier et al. [24] a similar approach wa s used except that the authors generated their estimated motion vector field by using the 2D opt ical flow method of Horn and Schunck [29] on initial filtered backprojection (FBP) reconstructions of the gated imag es. In Gravier et al. [65] the authors extended their method to 3D and evalua ted the resulting images in terms of lesion detection using a channelized Hotelling observer. It was found that the 4D motion compensated images had superior lesion detection properties co mpared to images without temporal smoothing. However, in both [24] and [65] the authors noted that the accuracy of the non-rigid motion

PAGE 62

62 vector fields used in their MAP reconstructions suffered from the already high noise levels of the FBP images. Simultaneous Image Reconstruc tion and Motion Estimation Motivated by the inherent problem of highnoise motion estimation in [24] and [65], a method which simultaneously estimates both th e tomographic images and non-rigid motion vector fields has been proposed [10],[12],[ 66]. The method is described by the following objective function: f,m(f)f,mmISELEE (5-1) where f is the 4D tomographic image, m is the 4D motion vector field, 11log fp J j jj jiLH figiH fi (5-2) is the negative log likelihood of obtaining the detector data from the reconstructed image f, Hf ( i ) is the projection operator, g is the measured detector data vector, J is the total number of frames, p is the total number of detector bins, 2 1 1 J Ijjj jEff rf,m rrmr (5-3) is an image registration term taken over all voxels r, 222 1 J sjjj jEuvwrmrrr (5-4) is a priori information about the obj ect motion (in this case a smoothness constraint on the estimated motion vector field), u,v,w are the components of the estimated motion vector field m, and is the gradient operator. In eq. 5-1 and are user-defined scalars which control the influence of the reconstruction and motion smoothi ng terms, respectively. It should be noted that Es can be used to enforce more stringent models of the object motion such as incompressibility,

PAGE 63

63 consistency with a prediction fiel d, or (as in Mair et al. [12]) the strain energy of a material undergoing elastic deformation. This simultane ous method has been shown to have improved motion estimation properties compared to a stan dard motion estimation method [67], as well as improved image quality compared to a standard r econstruction method in an SSE (eq. 3-1) sense [12]. However, an investigation of the detection properties of the simultaneous method has been notably absent from previous studies. Channelized Hotelling Observer Receiver operating characteristics studies usi ng human observers are considered the gold standard in evaluating the detection propertie s of imaging modalities [20]. However, these studies are time consuming and require the partic ipation of multiple observers. An alternative to human studies are mathematical observers designed to mimic the detectio n properties of humans [68]. The channelized Hotelling observer is one such model obs erver with a channel operator developed to model the human visual system. The method has been shown to correlate well with human observers in relatively difficult imagi ng tasks with both corre lated noise and random backgrounds [62]. A brief description of the channelized Hotelling observer is given here; a more detailed description can be found in Yao and Barrett [21] and Myers and Barrett [62]. The CHO uses a set of training images to bui ld a covariance matrix which describes the statistical differences between si gnal absent (SA) and signal presen t (SP) classes. It has been suggested that the number of images used for trai ning should be at least e qual to the order of the covariance matrix. After the covariance matrix has been generated, the CHO can calculate a test statistic for a given image which can be compared to a decision threshold to classify the image as normal or abnormal. What separates the CHO fr om other mathematical observers (e.g. the nonprewhitening observer) is a pre-filter process which models the human visual system by operating on the power of select frequency bands, rather than the full frequency spectrum. This

PAGE 64

64 pre-filter stage is referred to as a channel mechanism and its inclusion into the Hotelling observer has demonstrated good correlat ion with human observers [69]-[71]. We represent each of the reconstructed images as a 1D vector,g of length N. The CHO test statistic, CHO is defined as: 1 2 T CHOSPSAS (5-5) where SP is the mean feature vector of a set of SP images, SA is the mean feature vector of a set of SA images, T is the transpose operator, 1 2S is the inverse of the average inter-class scatter matrix (2S ), and is the feature vector of the current image. The average inter-class scatter matrix is defined by: 2,,1 2SPSASSS (5-6) where ,iS is the average intra-class scatter matrix, defined by: T iii iS (5-7) where the operator represents an average over class i. The CHO operates using a number of frequency bands, L. The feature vector is defined by: Ug (5-8) where U is the LxN channel operator matrix. We use 4 rotationally symmetric frequency bands or channels, with each band one octave wi de ranging from 0.03125 to 0.0625 cycles/pixel, 0.0625 to 0.125 cycles/pixel, 0.125 to 0.25 cycles/p ixel, and 0.25 to 0.5 cycles/pixel. Fig. 5-1 shows a plot of the four channe ls. These channels have been shown to correlate well with human observers [71]. The feature vector v can be calculated in either the frequency or spatial domain; here we use the spatial domain approach to av oid taking a large number of Fourier transforms.

PAGE 65

65 Spatial domain templates for U may be calculated by taking the i nverse Fourier tr ansform of the frequency bands. The spatial domain templates may then be used each time v is calculated. The frequency bands and their correspon ding spatial domain templates are shown in Fig. 5-2. In this paper, the center of the defect must be aligned with the center of the spatial domain template. All images were scaled to integer intensities betw een 0 and 255 to mimic the conditions of a human ROC study. The CHO SNR d may be used as an image quality figure of merit: 22'varvarSPSA CHO SPSAd (5-9) where SP is the mean of CHO given class SP, SA is the mean of CHO given class SA, varSP is the variance of CHO within class SP, and varSA is the variance of CHO within class SA. The area under the RO C curve (AUC) is related to d by: 11 AUCerf 222 d (5-10) where erf is the Gaussian error function for d : 202 erf'(')d tdedt (5-11) Multiple-Array Channelized Hotelling Observer In Chen et al. [63] it was hypothesized that in a clinical settin g, human observers would consider both multiple slices and multiple view s in the detection process, and that a CHO modeling this behavior would reflect human obser vers more accurately. The authors provided a mathematical framework to extend the CHO to op erate on 4D images. However, the authors did

PAGE 66

66 not find a significant difference in correlation with human observers when using the 4D approach compared to the 2D approach. More recently our group has developed an exte nsion of the CHO for use in gated imaging [64] that is based on the framework provided by Chen et al. [63]. The difference between our method and that of [63] is that we operate on mu ltiple frames, as opposed to multiple slices and multiple views. We propose the use of this method for ECG-gated tomographic image sets because it takes into account the fu ll series of gated images, more closely modeling a clinical evaluation. We refer to the method as the multiple-array channelized Hotelling observer (MACHO). The first step of the MACHO is to calculate the usual CHO test statistic for each 2-D slice of each image frame. This results in a vector of test statistic, CHO of length N for each noisy reconstructed image. In th is study a single short axis slice through th e center of the lesion location in each image frame was used, and thus N was equal to the number of ECG-gated frames. The second step uses a Hotelling observer (HO) to generate an ov erall test statistic, MACHO, for each noisy reconstruction, from th e collection of CHO test statistics CHO The HO is similar to the CHO, but w ithout the pre-fi ltering stage: 1 ,, 2CHOT M ACHOCHOSPCHOSA CHOS (5-12) where CHOSP is the mean vector of CHO over all SP training images, CHOSA is the mean vector of CHO over all SA training images, 1 2CHOS is the inverse of the aver age inter-class scatter matrix (2CHOS), and CHO is the vector containing a CHO for each slice of the current image. The average inter-class scatter matrix is defined by:

PAGE 67

67 2,,1 2CHO CHO CHOSP SASSS (5-13) where ,CHOiS is the average intra-class scatter matrix, defined by: ,, ,CHOT CHOCHOiCHOCHOi i iS (5-14) The MACHO SNR, d, is defined as: 22'varvarSPSA MACHO SPSAd (5-15) Methods Simulation of Test Data The 4D NURBS-based Cardiac Torso (NCAT) phantom [59] was used in this study. An activity distribution representing an ECG-gated 99mTc-Sestamibi rest SPECT scan of a normal (defect-free) human torso was ge nerated with 32 frames over one complete cardiac contraction cycle. Relative organ activity concentrations were taken from Gilland et al. [60] and are given in Table 5-1. The 32 frames were then down-sampled to 8, providing realistic within-frame motion blur in the ECG gated images. Additionally, an activity distri bution representing a myocar dial perfusion defect was generated. The defect, shown in Fig. 5-3 (after subtracting from the norma l activity distribution), had a circumferential width of 120, length of 31 .25 mm, and transmural thickness equal to the thickness of the myocardium. This defect is cons idered to be a large sized defect [72]. The cardiac contraction motion for the defect was sa mpled at 32 frames and down-sampled to 8. The SIMIND Monte Carlo program [61] was used to generate realistic projection data from the normal activity distribution and the defe ct activity distribution, separately. After the simulation, the defect projection data could be subtracted from the normal projection data to

PAGE 68

68 create a projection data set with a myocardial perfusion defect. The simulations of the normal and defect projection data were carried out sepa rately to avoid simulating the extra-myocardial activity (a time consuming procedure) more th an once. Also, separate simulations made it possible to adjust the defect contrast in the ev aluation of the CHO code (section 5.2.C). Scatter, attenuation, detector stopping power, and detector response were modeled in both simulations assuming a gamma energy of 140 keV, crystal th ickness of 1.27 cm, pixel size of 0.3125 cm, 64 projections over 180 (LPO to RAO), 128x128 detector elements, and a low energy general purpose parallel hole collimator (3mm across hexagonal holes with 4.2mm pitch). The simulations were allowed to run until noise in the projection data was negligible ( 2%). OSEM reconstructions of the noise -free projection data with and w ithout the defect are shown in Fig. 5-4. In the figure, the pixel intensity in th e location of the defect is scaled to 95% of the normal myocardial wall pixel intens ity (referred to as a 5% defect). Validation of the CHO Prior to the com parison of RM and PS-OSEM, we tested the CHO across a wide variety of reconstructed images. The purpose of this step was to ensure the code could be used as an indicator of image quality. Six test s were performed. Each test us ed 500 noise realizations of SP and 500 noise realizations of SA data with reconstruction usin g 5 iterations of OSEM with 8 subsets. Test 1: The CHO SNR (eq. 5-9) acro ss a variety of defect intensities was tested. The defect projection set was scaled from 50% to 97. 5% of the intensity of the normal myocardial wall and subtracted from the normal projection se t, creating projection se ts with myocardial perfusion defects ranging from 2.5% to 50%. The intensities of the projection sets were then scaled to 14,000 counts/frame within a 3mm mi d-ventricular slice, and Poisson noise was simulated. This count level is representati ve of a clinical re st SPECT scan using 99mTc-Sestamibi

PAGE 69

69 and for this paper will be referred to as the clinical count level. The images were smoothed with a 3D Butterworth filter with cut-off freque ncy equal to 0.15 cycles/pixel. This level of smoothing was chosen to be representative of that which we would find in a clinical setting [73]. The response of the CHO across the range of defect contrasts was analyzed using eq. (5-9). Test 2: The CHO SNR across a variety of noise levels was tested. A projection set with a single 25% defect was generated, and the total intensity was s caled (from the clinical count level) by factors of 0.1 to 1.5. This scali ng was followed by simulated Poisson noise. At each level of scaling the images were smoothed with a 3D Butterworth filter with cut-off frequency equal to 0.15 cycles/pixel, and the CHO SNR was calculated using eq. (5-9). Test 3: The response of the CHO SNR to the numb er of channels used was tested. A set of reconstructed images with a sing le 25% defect at the clinical count level we re used to test the CHO with the number of channels used ranging from 1 to 4. Test 4: The response of the CHO SNR to the numb er of training images used was tested. This test was used to ensure that the number training images used in the ROC evaluation would yield statistically significant resu lts. A set of reconstructed imag es with a single 25% defect at the clinical count level were us ed to test the CHO with traini ng image quantities ranging from 10 to 250. Test 5: The response of the CHO SNR to varying levels of spatial smoothing was tested. A set of reconstructed images w ith a single 25% defect and the cl inical count level was tested with a 3D Butterworth filter with cut-off fr equency ranging from 0.0 to 0.3 cycles/pixel. Test 6: The response of the CHO SNR to vary ing levels of temporal smoothing was investigated. A set of reconstructed images with a single 25% defect and th e clinical count level was tested with a 3D Butterworth filter with cut-off frequency 0.15 cycl es/pixel and a 3-point

PAGE 70

70 temporal convolution over a range of kernel weights. The temporal smoothing kernels are defined in Table 5-2. Image reconstruction For the com parison of the RM and PS-OS EM methods, a set of 500 SA and 500 SP projection sets at the clinical count level were created from the noise-free simulations. The SP sets had a 5% defect. Based on the results of Test 5 and Test 6 in section 5.2.C, we chose to evaluate the RM and PS-OSEM methods at equa l levels of spatial resolution. The spatial resolution was quantified as the full-width at half-maximum (F WHM) of a profile through the left-ventricular wall of a noise-free, defect-free reconstruction in the short axis view. To ensure this was a valid measure of spatial resolution, it was tested on MLEM images as a function of MLEM [11] iteration number, 3D Butterworth cu t-off frequency, and temporal convolution kernel. PS-OSEM reconstruction: For the PS-OSEM reconstructions, a spatial resolution level similar to clinical images was defined as 5 iter ations of OSEM with 8 subsets followed by a 3D Butterworth filter with a cut-off frequency of 0.15 cycles/pixel [73] and a 3-point temporal convolution with weights {.25,.5,.25} (temporal c onvolution kernel #4). The FWHM of a noisefree PS-OSEM image at this spatial resolution was calculated. These settings were then used to reconstruct the 500 noise realizations of SA and SP projection sets. RM reconstruction: In this paper the reconstructionstep (R-step) of RM was a MAP reconstruction which used a modification of Greens one step late algorithm [74]. The implementation of the R-step is detailed in Mair et al. [12]. The motion-step (M-step) in this paper was a generalization of the Horn and Sc hunck optical flow algorithm into a sequential quadratic framework [75]. The parameters and were determined empirically using a noise-

PAGE 71

71 free SA projection set to generate a spatial resolution equal to the PS-OSEM images. These settings were then used to reconstruct the 500 noise realizations of SA and SP projection sets. ROC evaluation The CHO was used to generate a set of test statistics for both reconstruction m ethods. The CHO test statistics were calculated using a central slice (through the defect center) at the end-diastolic frame of the gate d reconstructions. The ROCKIT so ftware package developed at the University of Chicago [20] was used to generate a maximum likelihood estimation of the parameters of the bivariate bino rmal ROC model for the paired data as well as to calculate the statistical significance of the difference be tween the binormal ROC curve estimates. An identical study was also performed us ing the MACHO in place of the CHO. The MACHO test statistics were calculated using a central slice of all frames (8) of the gated reconstructions. Direct assessment of SNR using regions-of-interest across the ensemble data The SNR for both RM and PS-OSEM was also calculated directly using a regions-ofinterest technique [76] across the SP ensemble data. The average pixel intensity in the area of the defect, L, was calculated using a region-of -interest defined by the margins of the defect (Fig. 5-5) in the central slice of the end-diastolic frame. Then, in three consecutive slices of the enddiastolic frame, two regions-of-interest with sh ape identical to the defect ROI were defined within the myocardial activity but outside the area of the defect. The average pixel intensity of each lesion-free ROI, bi, was calculated (where i=1,2,,6). The SNR was then defined as: iROI b B L SNR (5-16) where B is the average of bi, and ib is the standard deviation of b across all i. The final value of the SNRROI for both RM and PS-OSEM was computed as an average over the ensemble of SP

PAGE 72

72 data (500 noise realizatio ns). The significance of the difference in SNRROI between the RM and PS-OSEM methods was evaluated using a paired t-test. For comparison to the results of Test 5 using the CHO, the response of the SNRROI was also tested across a ra nge of 3D Butterworth cutoff frequencies. Reconstruction without the effects of scatter attenuation, and detector response We also compared the RM and PS-OSEM methods qualitatively using projection data obtained from the phantom by a simple linear pr ojector that ignored the effects of scatter, attenuation, and detector response. Poisson noise at the clinical count level was simulated in the projection set and RM and PS-OSEM reconstruction parameters were optimized using the SSE (eq. 3.1) between the original, noise-free phantom and the reconstructed images. Results Validation of the CHO Test 1: Fig. 5-6 shows the results found by testi ng the CHO at different levels of defect contras t. The figure shows the CHO SNR increases as the defect contrast is increased. This is in agreement with previous studies [72] and validat es the codes ability to rank image quality in terms of defect detection. Test 2: Fig. 5-7 shows the results found by testi ng the CHO at varying count levels. The plot shows that the CHO SNR increases as the to tal counts increase. This is in agreement with our expectation that detection should improve as image noise improves. The study further supports the codes ability to rank image quality in te rms of defect detection. Test 3: Fig. 5-8 shows the results of varying the number of channels used by the CHO from 1 to 4. As expected as the number of cha nnels increases, the CHO SNR also increases. This is due to the inclusion of more info rmation as more channels are used.

PAGE 73

73 Test 4: Fig. 5-9 shows the results found by varying the number of images used for training. There is a sharp incline in SNR betw een 10 and 30 images, after which point the plot levels out. We would expect a sharp incline betw een 0 and 16 images, which is the order of our covariance matrix. However, several authors have pointed out that althou gh this rule of thumb guarantees the existence of the inverse scatter matrix, it often underestimates the number of images needed for a statistically significant study [77]-[80], which may explain the continued increase up to ~30 training images. Test 5: Fig. 5-10 shows the results of the CHO SNR versus Butterworth cut-off frequency. There is a sharp increase in SNR between cut-off frequencies of 0. and 0.025 cycles/pixel, followed by a gradual and consistent decrease. Images representing four sample points in Fig. 5-10 are shown in Fig. 5-11. It is clear from the images that the level of smoothing which generates the highest CHO SNR (cut-off = 0.025 cycles/pixel) is considerably smoother than typical clinical images. These results do not agree with those reported in Frey et al. [72]; however, a similar discrepancy was noted by anothe r author [81]. Because the level of spatial smoothing is not normally optimized in papers using the CHO, at this time it is unclear if the results of this test are unique. Test 6: Fig. 5-12 shows the results of the CHO SNR versus temporal smoothing kernel. The results display a steady increase in CHO SN R as temporal smoothing weights are increased. This is in agreement with the findings of Test 5. Images representing the 5 levels of temporal smoothing are shown in Fig. 5-13. Image Reconstruction The results of the FW HM versus MLEM ite ration number investigation for noise-free data are shown in Fig. 5-14. As expected the spatial resolution of the image improves rapidly with the first few MLEM iterations (iterations 110) but then begins to level off (iterations 20-

PAGE 74

74 50). Fig. 5-15 shows the results of the FWHM versus Butterworth cut-off frequency. The spatial resolution again improves rapidly at first but then begins to level off. Figure 5-16 shows the results of FWHM versus temporal smoothing kernel The spatial resolution degrades as the level of temporal smoothing is raised. These experime nts validated the use of the FWHM of the LV wall to quantify the level of spatial resolution in the images. The FWHM of the PS-OSEM image after 5 iterations of OSEM with 8 subsets, Butterworth filter with cut-off frequency 0.15 cycles/pixel, and temporal convolution kernel number 4 was 60.5 pixels. The RM parameters found to yield an equivalent FWHM were =0.13 and =0.05. Example images for the PSOSEM and RM reconstructions of the noise free data and their correspond ing profiles are shown in Fig. 5-17. The estimated motion vector fiel ds generated by RM on the noise-free data are shown in Fig. 5-18. The vector fields dem onstrate a good balance between motion magnitude and motion smoothness. Example reconstructed images for the RM and PS-OSEM methods used in the CHO and MACHO study are shown in Fig. 5-19. The PS-OSEM images appear to have better contrast and uniformity of the activity distribution compared to the RM images. This is in contrast to results found in previous papers which used simpler mo dels of the system matrix and phantom object and higher spatial resolution in the reconstructed images [10],[12],[66]. ROC Evaluation The ROC curve generated using the CHO is s hown in Fig. 5-20. The curve for PS-OSEM stays above RM at all areas of the plot, representing superior de tection regardless of threshold effects. The AUC for RM was 0.69, com pared to 0.72 for PS-OSEM. The p-value calculated from the correlated bivariate chi-square test st atistic was 0.54 (46% conf idence level). This result demonstrates that when using a si ngle short-axis slice, it cannot be said that the detection of the RM and PS-OSEM methods differs.

PAGE 75

75 Fig. 5-21 shows the ROC curve generated usi ng the MACHO. In contrast to Fig. 5-20, the RM curve stays above PS-OSEM at all ar eas of the plot. The AUC for RM was 0.85, compared to 0.83 for PS-OSEM. The p-value calcula ted from the correlated bivariate chi-square test statistic was 0.80 (20% confidence level). This result demonstrates that when considering all frames of a gated image, it cannot be said that the detection characteristics of the RM and PSOSEM methods differ. Direct Assessment of SNR Using Regionsof-In terest across the Ensemble Data The signal-to-noise ratios (SNRROI) calculated using the regions-of-interest technique were 0.50 for RM and 0.44 for PS-OSEM. The difference between the two methods was found to be significant at the 1% le vel. The response of SNRROI to cut-off frequency is shown in Fig. 522. The results indicate that SNRROI is more dependent on defect contrast than image noise. This is in contrast to the behavior of the CHO, which seems to be heavily dependent on image noise rather than defect contrast. Reconstruction without the Effects of Atte nuation, Scatter, and Detector Response The reconstructions of the simplified linear projection simulation are shown in Fig. 5-23. In contrast to the results of sec tion 5.3.C, the RM reconstructions appear to have better contrast and uniformity of the activity distribution th roughout the LV wall compared to the PS-OSEM images. The findings of this simplified study are in agreement with previous papers [10],[12],[66]. The results suggest that RM may improve upon PS-OSEM in idealized situations, but that these improvements are less dramatic wh en scatter, attenuation, and detector response are included. Similar results were found in Ch 3 (Fig. 3-7) which used a physical dynamic phantom that modeled the mo tion of a beating heart.

PAGE 76

76 Summary and Discussion In this work we evaluated the detecti on characteristics of a simultaneous image reconstruction and motion estimation method in term s of lesion detection compared to a standard reconstruction method. A channelized Hotelling observer was used to evaluate reconstructed images from both methods on a central slice in th e short-axis view at th e end-diastolic frame. The difference in the AUC calculated from the CHO test statistics for the simultaneous and standard methods was not considered to be stat istically significant. The methods were further evaluated using a CHO modified for use on gate d images. When considering all 8 frames, the difference between the AUC for the simultane ous method and the standard method was not found to be statisti cally significant. The SNR of the methods was also evaluated directly using a regions-of-interest method. Using this method it was found that th e simultaneous method had a higher SNRROI than the standard method. The results were found to be significant at the 1% level. There is substantial data relating detectability to SNR [82]-[85] and it may be said that the findings of this ROI SNR investigation show that RM is a better method for lesion detection in cardiac SPECT than the conventional PS-OSEM. However, the results of th is study conflict with the results of the CHO studies. We believe that based on the results of Test 5 and Test 6 (section 5.3.A and Figs. 5-9 and 5-11) the CHO does not correlate well with the de tection characteristics of human observers in the case of varying levels of image smoothness. Th is is an important finding because it indicates that the CHO cannot be used as a measur e of image quality when the methods under investigation have different levels of smoothing. An area of future work will be to further develop and compare the ROI method to human observers. At least one other author [72] has demonstrated results whic h conflict with the findings of Test 5 and Test 6. There are several possible reasons for this discrepancy. The phantom

PAGE 77

77 population used in Frey et al. [72] was based on mathematically defined anatomies (ours was based on the Visible Human Project) and the simula tion of projection data was carried out using an analytical operator, compared to the Monte Carlo method used here. The defects in [72] were also smaller than in our study, and it may be reas onable to assume that larg er defects can lead to higher levels of optimal smoothing using the CHO due to the larger difference in average image intensity. This discrepancy will be the subject of future work.

PAGE 78

78 .03 .06 .125 .25 .5 1 Frequency (cycles/pixel)Intensity (arb.) Figure 5-1. Plot of the four fr equency bands used in the CHO. Figure 5-2. The four frequency bands used in the CHO with the zero frequency located at the center of the image (top) a nd their corresponding spatial domain templates (bottom).

PAGE 79

79 Figure 5-3. Phantom source dist ributions showing the shape a nd location of the defect. The images are shown without ex tra-myocardial activity. Figure 5-4. Noise-free reconstructions for the normal (left) and abnormal (right) datasets. Figure 5-5. ROIs used to calculate SNR directly.

PAGE 80

80 0 5 10 15 20 25 30 35 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 d primeDefect contrast (%) Figure 5-6. d vs. defect contrast. 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 d primeCount level (fraction of clinical level) Figure 5-7. 'd vs. count level.

PAGE 81

81 1 2 3 4 0 0.5 1 1.5 2 2.5 3 3.5 d primeNumber of channels Figure 5-8. 'd vs. number of channels used. 0 50 100 150 200 250 0 0.5 1 1.5 2 2.5 3 3.5 d primeNumber of training images Figure 5-9. d vs. number of training images used

PAGE 82

82 0 0.05 0.1 0.15 0.2 0.25 0.3 0 1 2 3 4 5 6 d primeButterworth cutoff frequency (cycles/pixel) Figure 5-10. 'd vs. Butterworth cut-off frequency Figure 5-11. Select imag es used in Figure 5-10.

PAGE 83

83 1 2 3 4 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 d primeTemporal smoothing kernal Figure 5-12. 'd vs. temporal smoothing kernel number. The weights for each kernel are defined in section 5.2.C. Figure 5-13. Example images of the five levels of temporal smoothing evaluated in Test 6.

PAGE 84

84 5 10 15 20 25 30 35 40 45 50 50 60 70 80 90 100 110 120 130 140 IterationsFWHM (pixels) Figure 5-14. FWHM vs MLEM ite ration number on noise-free data. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 50 60 70 80 90 100 110 120 130 140 Butterworth cutoff freq. (cycles/pixel)FWHM (pixels) Figure 5-15. FWHM vs Butterworth cut-off frequency on noise-free data.

PAGE 85

85 1 1.5 2 2.5 3 3.5 4 4.5 5 60.6 60.8 61 61.2 61.4 61.6 61.8 62 Temporal convolution kernalFWHM (pixels) Figure 5-16. The FWHM vs temporal convolution kernal on noise-free data 0 20 40 60 80 100 120 140 0 0.2 0.4 0.6 0.8 1 PixelIntensity (normalized) OSEM RM Figure 5-17. Example slices of the A) PS-OSEM and B) RM reconstructed images, and C) their corresponding profiles. A B C

PAGE 86

86 Figure 5-18. Motion vector fields for RM reconstructions in Fig. 5-15 at max-systole (top) and max-diastole (bottom). Areas enclosed with yellow boxes are zoomed in on the right.

PAGE 87

87 Figure 5-19. Example RM and PS-OSEM images used in the ROC study. 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 False Positive FractionTrue Positive Fraction RM (AUC=0.69) PS-OSEM (AUC=0.72) Figure 5-20. ROC curve for RM and PS-OS EM using 2D CHO on a single lesion.

PAGE 88

88 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 False Positive FractionTrue Positive Fraction RM (AUC=0.85) PS-OSEM (AUC=0.83) Figure 5-21. ROC curve for RM and PS-OSEM using MACHO on a single lesion. 0 0.05 0.1 0.15 0.2 0.25 0.3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 SNRROIButterworth cutoff frequency (cycles/pixel) Figure 5-22. SNRROI vs. 3D Butterworth cut-off frequency.

PAGE 89

89 Figure 5-23. Reconstructions of projection data obtained using a simple linear projection operator. RM images appear to have impr oved contrast and uniformity of the activity distribution throughout the LV wall. Table 5-1. Relative organ activity concentrations per pixel for the simulated phantom Organ Relative activity Myocardium 75 Heart blood pool 6 Liver 13 Gall bladder 324 Lung 6 Kidney 45 Spleen 45 Bowel 37 Background (body) 6 Table 5-2. Tempor al smoothing kernels Kernel number Weights Numb er of applications 1 {0.0,1.0,0.0} 2 {0.1,0.8,0.1} 1 3 {0.1,0.8,0.1} 2 4 {0.25,0.5,0.25} 1 5 {0.33,0.34,0.33} 1

PAGE 90

90 CHAPTER 6 CONCLUSIONS This dissertation investigat ed image reconstruction methods incorporating estimated motion in gated imaging. The primary met hod investigated was a simultaneous image reconstruction and motion estimation appro ach. In Chapter 3 a physical, dynamic cardiac phantom was modified by attaching point source markers to the phantom epicardial wall. The markers were used to provide an independent measure of the heart wall motion. Acquisitions of the phantom myocardium were r econstructed and the frame-to-frame motion was estimated with the simultaneous method, as well as with a conventional method. It was found that the simultaneous method produced motion estimates with greater moti on accuracy than the conventional method. As mentioned in Chapter 3, this study ignored the wringing, or twisting, motion of the heart. A future area of investiga tion in this area will be an adaptation of the physical phantom to mimic a level of wringi ng motion found in the human heart. Also, the development of a heart assembly with greater sampling (upwards of 100 marker locations) at both the inner and outer surface of the myocardial chamber will be considered. This will require the development of smaller point source markers to account for the limited space at the inner wall of the phantom myocardium. Also, imaging m odalities with greater spatial resolution, such at CT and MR, will be investigated in determining the true marker motion between frames. In Chapter 4 a new rigid-body motion estim ation technique for use in estimating respiration-induced heart motion was developed. The new method used rotation parameterized by a unit quaternion, and the method was compared to several other methods which have been proposed in the literatu re. The quaternion method was found to be faster, more robust, and more accurate than other techniques across a wide ra nge of reconstructed images. Motion-corrected images using the quaternion method demonstrated improved spatial resolution over uncorrected

PAGE 91

91 images. A future direction of this study will be a receiver operating characteristics analysis comparing the CGQ method with a standard reco nstruction method which i gnores the effects of respiratory blurring of the heart. In Chapter 5 the simultaneous image reconstruction and motion estimation method was compared to a standard reconstruction method in clinical practice in te rms of lesion detection using a channelized Hotelling observer formulated to operate on singleand multiple-frame reconstructed images. The signalto-noise ratio of the simultane ous method was also calculated directly using a regions-of-inte rest method and compared to th e standard method. In both the singleand multiple-frame observer studies the differences in detection between the simultaneous method and the standard method were not statistically significant. However, using the regions-of-interest method, the simultaneous method was found to have superior detection properties compared to the sta ndard method. The difference was found to be significant at the 1% level. A future direction of this final study will be an investigation into the discrepancy between the CHO and human observers under varyi ng levels of image smoothness. We will also further investigate the SNR calculation using regions-of-interest as a model of the detection characteristics of human observers.

PAGE 92

92 APPENDIX A ROTATION MATRIX AND PART IAL DERIVATIVES In the Euler representation, Q is the product of the rotation matrices about each axis: 100cos0sincossin0 Q(,,)QQQ0cossin010sincos0 0sincossin0cos001xyz (A-1) The orthogonal (rotation) matrix Q corresponding to a rotation by ( ) is: coscos cossin sin Q(,,)sinsincoscossinsinsinsincoscoscossin cossincossinsincossinsinsincoscoscos (A-2) In computing the partial derivatives of our objective function, the products QQQ x yz, QQQ x yz and QQQ x yz are needed: sincos coscos 0 QQQcoscossinsinsincossinsinsincos0 cossinsincossincoscossinsinsin0xyz (A-3) sincossinsincos QQQcoscossincossinsinsinsin coscoscoscossincossincosxyz (A-4) 000 QQQcoscossinsinsincossinsinsincoscoscos cossinsincossincoscossinsinsinsinsincosxyz (A-5)

PAGE 93

93 APPENDIX B DETERMINATION OF ROTATION MATRIX The orthogonal (rotation) m atrix R corresponding to a rotation by q is: 2222 012312031302 2222 0123120301232301 2222 1302230101232()2() Q(,,,)2()2() 2()2() (B-1) With the substitution 222 01231 in place, the part ial derivatives of Q are: 13 12 23 00 2 13 1 213 100 2 121 301 00022 Q 242 224 (B-2) 2 23 2 210 00 23 21 13 200 2 21 2 032 00422 Q 202 224 (B-3) 2 332 301 00 2 331 032 300 3231 12 00422 Q 242 220 (B-4)

PAGE 94

94 APPENDIX C PARTIAL DERIVATIVES USED IN CGQ ALGORITHM The partial d erivatives of E with respect to each of the parameters are: 1 11(Q,) Q 2()(Q+)Q+jjE fff rb rrbrbr (C-1) 1 22(Q,) Q 2()(Q+)Q+jjE fff rb rrbrbr (C-2) 1 33(Q,) Q 2()(Q+)Q+jjE fff rb rrbrbr (C-3) 1(Q+ (Q,) 2()(Q+)j j xf E ff bx rrb) b rrb (C-4) 1(Q+ (Q,) 2()(Q+)j j yf E ff by rrb) b rrb (C-5) 1(Q+ (Q,) 2()(Q+)j j zf E ff bz rrb) b rrb (C-6)

PAGE 95

95 REFERENCES [1] A.M. Minino, M.P. Heron, S.L. Murphy, a nd K.D. Kochankek, Deaths: final data for 2004. National vital statistics reports, National Center for Health Statistics, vol. 55, no. 19, 2007. [2] C.Y. Loong and C. Anagnostopoulos, Diagnos is of coronary artery disease by radionuclide myocardial perfusion imaging, Heart, vol. 90, pp. 2-9, 2004. [3] A.F. Schinkel, J.J. Bax, A. Elhendy, et al., Long-term prognostic value of dobutamine stress echocardiography compared with myocardial perfusion scanni ng in patients unable to perform exercise tests, Am. J. Med., vol. 117, pp. 1-9, 2004. [4] T. Marwick, B. Willemart, A.M. DHondt, et al., Selection of the optimal nonexercise stress for the evaluation of ischemic regi onal myocardial dysfunction and malperfusion. Comparison of dobutamine and aden osine using echocardiography and 99mTc-MIBI sinlge photon emission computed tomography, Circulation, vol. 87, pp. 345-354, 1993. [5] R.D. Des Prez, L.J. Shaw, R.L. Gillespie, W. A. Jaber, G.L. Noble, P. Soman, D.G. Wolinsky, and K.A. Williams, Cost-effectiv eness of myocardial perfusion imaging: A summary of the currently available literature, J. Nucl. Cardiol., vol. 12, pp. 750-759, 2005. [6] B.H. Duncan, A.W. Ahlberg, M.G. Levine, et al., Comparison of electrocardiographic gated technetium-99m sestimibi single-photon emission computed tomographic imaging and rest-redistribution thal lium-201 in the prediction of myocardial viability, Am. J. Cardiology, vol. 85, pp. 680-684, 2000. [7] T.M. Bateman, D.S. Berman, G.V. Heller, et al., American Society of Nuclear Cardiology position statement on electrocardiographic gating of myocardial perfusion SPECT scintigrams, J. Nucl. Cardiol., vol. 6, pp. 470-471, 1999. [8] D.S. Lalush and B.M.W. Tsui, Simulation ev aluation of Gibbs prior distributions for use in maximum a posteriori SPECT reconstructions, Trans. Med. Imag., vol. 11, pp. 267275, 1990. [9] D.S. Lalush, L. Cui, and B.M.W. Tsui, A priori models for four-dimensional reconstruction in gated cardiac SPECT, Nuc. Sci. Symp. Conf. Rec., vol. 3, pp. 1923-1927, 1996. [10] D.R. Gilland, B.A. Mair, J.E. Bowsher, and R.J. Jaszczak, Simultaneous reconstruction and motion estimation for gated cardiac ECT, IEEE Trans. Nucl. Sci., vol. 49, no. 5, pp. 2344-2349, 2002. [11] L.A. Shepp and Y. Vardi, Maximum likeli hood reconstruction for emission tomography, IEEE Trans. Med. Img., vol. MI-1, no. 2, pp. 113-122, 1982.

PAGE 96

96 [12] B.A. Mair, D.R. Gilland, and J. Sun, Estimation of images and nonrigid deformations in gated emission CT, IEEE Trans. Med. Imag., vol. 25, no. 9, pp. 1130-1144, 2006. [13] H. Hudson and R. Larkin, Accelerated imag e reconstruction using ordered-subsets of projection data, IEEE Trans. Med. Imag., v. 13, no. 4, pp. 601-609, 1994. [14] M.M. Ter-Pogossian, S.R. Bergmann, and B.E. Sobel, Influence of cardiac and respiratory motion on tomographic reconstr uctions of the heart: implications for quantitative nuclear cardiology, J. Comp. Asst. Tomo., vol. 6, pp. 1148-1155, 1982. [15] Y. Wang, S.J. Riederer, and R.L. Ehman, R espiratory motion of the heart: Kinematics and the implications for the spatial resolution in coronary imaging, Magn. Reson. Med., vol. 33, pp. 713-719, 1995. [16] K. McLeish, D.L.G. Hill, D. Atkinson, J.M. Blackall, and R. Razavi, A study of the motion and deformation of the heart due to respiration, IEEE Trans. Med. Img., vol. 21, no. 9, pp. 1142-1150, 2002. [17] P.P. Bruyant, M.A. Gennert, G.C. Speckert, R. D. Beach, J.D. Morgenstern, N. Kumar, S. Nadella, and M.A. King, A robust visual tracking system for patient motion detection in SPECT: Hardware solutions, IEEE Trans. Nucl. Sci., vol. 52, pp. 1288-1294, 2005. [18] G. Kovalski, O. Israel, Z. Keidar, A. Frenkel, J. Sachs, and H. Azhari Correction of heart motion due to respiration in clinical myocar dial perfusion SPECT scans using respiratory gating, J. Nuc. Med., vol. 48, no. 4, 630-636, 2007. [19] R.D. Fiete, H.H. Barrett, W.E. Smith, and K.J. Myers, Hotelling trace criterion and its correlation with human-observer performance, J. Opt. Soc. Am. A., vol. 4, no. 5, pp. 945953, 1987. [20] C.E. Metz, Basic principles of ROC analysis, Seminars in Nucl Med., vol. 8, pp. 283298, 1978. [21] J. Yao, and H.H. Barrett, Predicting hu man performance by a channelized Hotelling observer model, Mathematical Methods in Medical Imaging, SPIE, vol. 1768, pp. 161168, 1992. [22] R.L. DeValois and K.K. DeValois, Spatial vision, Oxford University press, New York, 1988. [23] D.S. Lalush and B.M.W. Tsui, Block iterative techniques for fast 4D reconstruction using a priori motion models in gated cardiac SPECT, Phys. Med. Biol., vol. 43, no. 4, pp. 875 886, 1998. [24] E. Gravier and Y. Yang, Motion-compensa ted reconstruction of tomographic image sequences, IEEE Nuclear Science Symposium Conference Record., vol 4, pp. 2927, 2003.

PAGE 97

97 [25] S.M. Song and R.M. Leahy, Computation of 3D velocity fields from 3D cine CT images of a human heart, IEEE Trans. Med. Imag., vol. 10, pp. 295, 1991. [26] G.J. Klein, B.W. Reutter, and R.H. Hues man, Nonrigid summing of gated PET via optical flow, IEEE Trans. Nucl. Sci., vol. 44, no. 4, pp. 1509, 1997. [27] G.J. Klein and R.H. Huesman, Elastic mate rial model mismatch effects in deformable motion estimation, IEEE Trans. Nucl. Sci., vol. 47, no. 3, pp. 1000, 2000. [28] G.J. Klein and R.H. Huesman, Four-dimensi onal processing of deformable cardiac PET data, Medical Image Analysis, vol. 6, pp. 29, 2002. [29] B.K. Horn and B.G. Schunck, Determining optical flow, Artif. Intell., vol. 17, pp. 185203, 1981. [30] H.M. Hudson and P. Larkin, Accelerated EM reconstruction using ordered subsets, Trans. Med. Img., vol. 13, no. 4, pp. 601-609, 1994. [31] A.A. Young and L. Axel, Three-dimensional motion and deformation of the heart wall: Estimation with spatial m odulation of magnetization a model-based approach, Radiology, vol. 185, pp. 241-247, 1992. [32] A. Mitiche and P. Bourthemy, Computation a nd analysis of image motion: A synopsis of current problems and methods, Int. J. Comput. Vision, vol. 19, no. 1, pp. 29-55, 1996. [33] L. Livieratos, K. Rajappan, L. Stegger, K. Schafers, D.L. Bailey, P.G. Camici, Respiratory gating of cardiac P ET data in list-mode acquisition, Eur. J. Nucl. Med., vol. 33, no. 5, pp. 584-588, 2006. [34] G. Shechter, C. Ozturk, J.R. Resar, E.R. McVeigh, Respiratory motion of the heart from free breathing coronary angiograms, IEEE Trans Med Imag., vol. 23, no. 8, pp. 1046 1056, 2004. [35] G.J. Klein, B.W. Reutter, M.H. Hol, J.H. Reed, R.H. Huesman, Real-time system for respiratory-cardiac gating in positron tomography, IEEE Trans Nucl Sci., vol. 45, no. 4, pp. 2139, 1998. [36] A. Martinez-Moller, D. Zikic, R.M. Botnar, R.A. Bundschuh, W. Howe, S.I. Ziegler, N. Navab, M. Schwaiger, S.G. Nekolla, Dual car diacrespiratory gated PET: implementation and results from a feasibility study, Eur. J. Nuc. Med. Mol. Img., vol. 34, no. 9, pp. 14471154, 2007. [37] E.J. Hoffman, M.E. Phelps, G. Wisenberg, H.R. Schelbert, D.E. Kuhl, Electrocardiographic ga ting in positron emission computed tomography, J. Comput. Assist. Tomogr., vol. 3, pp. 733-739, 1979. [38] N.C. Detorie, A.L. Kesner, T.D. Solberg, M. Dahlbom, Evaluation of image noise in respiratory gated PET, IEEE Trans. Nuc. Sci., vol. 54, no. 1, pp. 66-70, 2007.

PAGE 98

98 [39] R.D. Beach, P.H. Pretorius, G. Boening, P.P. Bruyant, B. Feng, R.R. Fulton, M.A. Gennert, S. Nadella, and M.A. King, Feasibility of stereo-infrared tracking to monitor patient motion during cardiac SPECT imaging, IEEE Trans. Nucl. Sci., vol. 51, no. 5, pp. 2693, 2004. [40] R.D. Beach, H. Depold, G. Boening, P.P. Bruyant, B. Feng, H.C. Gifford, M.A. Gennert, S. Nadella, M.A. King, An adaptive appro ach to decomposing patient-motion tracking data acquired during cardiac SPECT imaging, IEEE Trans. Nuc. Sci. vol. 54, no. 1, pp. 130-139, 2007. [41] P.P. Bruyant, M.A. King, P.H. Pretorius, Co rrection of the respiratory motion of the heart by tracking of the center of ma ss of thresholded projections : A simulation study using the dynamic MCAT phantom, IEEE Trans. Nucl. Sci., vol. 49, no. 5, pp. 2159-2166, 2002. [42] B. Feng, P.P. Bruyant, P.H. Pretorius, R.D. Beach, H.C. Gifford, J. Dey, M.A. Gennert, M.A. King, Estimation of the rigid-body moti on from three-dimensional images using a generalized center-of-ma ss points approach, IEEE Trans. Nucl. Sci., vol. 53, no. 5, pp. 2712-2718, 2006. [43] L. Livieratos, L. Stegger, P.M. Bloomfield, K. Schafers, D.L. Bailey, P.G. Camici, Rigidbody transformation of list-mode projection da ta for respiratory motion correction in cardiac PET, Phys. Med. Biol., vol. 50, pp. 3313-3322, 2005. [44] G.J. Klein, B.W. Reutter, R.H. Huesman, F our-dimensional affine registration models for respiratory-gated PET, IEEE Trans. Nucl. Sci., vol. 48, no. 3, pp. 756-760, 2001. [45] J. Dey, B. Feng, K.L. Johnson, J.E. McNamara, P.H. Pretorius, M.A. King, Respiratory motion correction in cardiac SPECT using a ffine and free-form defo rmation registration with temporal and spatial constraints, Proc. 9th Int. Meeting Fully3D Image Rec. Rad. Nucl. Med., pp. 201-204, 2007. [46] J.G. Parker, B.A. Mair, D.R. Gilland, M. Mahoney, Cardiac emission tomography with 3D respiratory motion correction, Proc. 9th Int. Meeting Full y3D Image Rec. Rad. Nucl. Med., pp. 413-416, 2007. [47] W. Lu, P.J. Parikh, J.P. Hubenschmidt, J.D. Bradley, D.A. Low, A comparison between amplitude sorting and phase-angle sorting usi ng external respiratory measurement for 4D CT, Med. Phys., vol. 33, no. 8, pp. 2964-2974, 2006. [48] H.W. Korin, R.L. Ehman, S.J. Riederer, J.P. Felmlee, R.C. Grimm, Respiratory kinematics of the upper abdominal organs: a quantit ative study, Magn Reson. Med., vol. 23, pp. 172-178, 1992. [49] P. Thevenaz, U.E. Ruttimann, M. Unser, A pyramid approach to subpixel registration based on intensity, IEEE Trans. Image Proc., vol. 7, no. 1, pp. 27-41, 1998. [50] V. Kantabutra, On hardware for computi ng exponential and trigonometric functions, IEEE Trans. Comp., vol. 45, no. 3, pp. 328-39, 1996.

PAGE 99

99 [51] B.K.P. Horn, Closed-form solution of abso lute orientation usi ng unit quaternions, J. Opt. Soc. of America A., vol. 4, pp. 629-642, 1987. [52] R. Mukundan, Quaternions: From classica l mechanics to computer graphics, and beyond, Proc. 7th Asian Tech. Conf. in Mathematics. Invited paper. Melaka, Malaysia, Dec. 18-21, 2002. [53] W.W. Hagar and H. Zhang, A new conjugate gradient method with guaranteed descent and an efficient line search, SIAM J. Optimization. vol. 16, pp. 170-192, 2005. [54] W.W. Hagar and H. Zhang, Algorithm 851: CG_DESCENT, A c onjugate gradient method with guaranteed descent, ACM Trans. Math. Soft. vol. 32, pp. 113-137, 2006. [55] T.L. Faber and E.M. Stokely, Orientation of 3-D structures in medical images, IEEE Trans. Pattern Anal. and Mach. Intel., vol. 10, no. 5, pp. 626-633, 1988. [56] N.M. Alpert, J.F. Bradshaw, D. Kennedy, J.A. Correia, The principal axes transformation A method for image registration, J. Nucl. Med. vol. 31, pp. 1717-1722, 1990. [57] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users Guide Third Edition, Knoxville, University of Tennessee, 1999. [58] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, Ca mbridge University Press, 1993. [59] W.P. Segars, D.S. Lalush, B.M.W. Tsui, "A Realistic spline-based dynamic heart phantom," IEEE Trans. Nucl. Sci. vol. 46, no. 3, pp. 503-506, 1999. [60] D.R. Gilland, R.J. Jaszczak, M.W. Hanson, K.L. Greer, R.E. Coleman, An experimental phantom based on quantitative SPECT analys is of patient MIBI biodistribution, J. Nucl. Med., vol. 37, no. 154P, 1996. [61] M. Ljungberg and S.E. Strand, A Monte-Ca rlo program simulating scintillation camera imaging, Comp. Meth. Progr. Biomed. vol. 29, pp. 257-272, 1989. [62] K.J. Myers and H.H. Barrett, Addition of a channel mechanism to the ideal-observer model, J. Opt. Soc. Amer., vol. 4, pp. 2447-2457, 1987. [63] M. Chen, J.E. Bowsher, A.H. Baydush, K.L. Gilland, D.M. DeLong, R.J. Jaszczak, Using the Hotelling observer on multislice and multiv iew simulated SPECT myocardial images, IEEE Trans. Nuc. Sci., vol. 49, no. 3, pp. 661-667, 2002. [64] Z. Cao, D.R. Gilland, B.A Mair, B.M.W. Tsui, An observer model evaluation of simultaneous reconstruction and motion estimation for emission tomography, IEEE Int. Symp. Biomed. Img.: Nano to Macro., vol. 1, pp. 940-943, 2004.

PAGE 100

100 [65] E. Gravier, Y. Yang, M.A. King, and M. Jin, Fully 4D motion-compensated reconstruction of cardiac SPECT images, Phys. Med. Bio., vol. 51, pp. 4603-4619, 2006. [66] Z. Cao, D.R. Gilland, B.A. Mair, R.J. Jaszczak, Three-dimensional motion estimation with image reconstruction for gated cardiac ECT, IEEE Trans. Nuc. Sci., vol. 50, no. 3, pp. 384-388, 2003. [67] J.G. Parker, D.R. Gilland, Wall motion estimation for gated cardiac emission tomography: physical phantom evaluation, IEEE Trans. Nuc. Sci., vol. 55, no. 1, pp. 531536, 2008. [68] H.H. Barrett, J. Yao, J.P. Rolland, K.J. Myers, Model observers for assessment of image quality, Proc. Nat. Acad. Sci. USA, vol., 90, pp. 9758-9765, 1993. [69] H.C. Gifford, M.A. King, D.J. de Vries, E. J. Soares, Channelized Hotelling and human observer correlation for lesion det ection in hepatic SPECT imaging, J. Nuc. Med., vol. 41, pp. 514-521, 2000. [70] K.J. LaCroix, B.M.W. Tsui, E.C. Frey, Rot ationally symmetric vs. oriented frequency channels for the Hotelling observer: A comparison with human observers, IEEE Nuc. Sci. Symp. Conf. Rec., vol. 3, pp. 1402-1406, 1999. [71] S.D. Wollenweber, B.M.W. Tsui, D.S. Lalush, E.C. Frey, K.J. LaCroix, G.T. Gullberg, Comparison of Hotelling observer models and human observers in defect detection from myocardial SPECT imaging, IEEE Trans. Nuc. Sci., vol. 46, no. 6, pp. 2098-2103, 1999. [72] E.C. Frey, K.L. Gilland, B.M.W. Tsui, Application of Task-Based Measures of Image Quality to Optimization and Evaluation of Three-Dimensional Reconstruction-Based Compensation Methods in Myo cardial Perfusion SPECT, IEEE Trans. Med. Img., vol. 21, no. 9, pp. 1040-1050, 2002. [73] D.R. Gilland, B.M.W. Tsui, W.H. McCartney, J.R. Perry, J. Berg, Determination of the Optimum Filter Function for SPECT Imaging, J. Nuc. Med., vol. 29, pp. 643-650, 1988. [74] P. Green, Bayesian reconstructions from emission tomography using a modified EM algorithm, IEEE Trans. Med. Img., vol. 9, no. 1, pp. 84-93, 1990. [75] D.R. Gilland, B.A. Mair, J.G. Parker, Motion estimation for cardiac emission tomography by optical flow methods, Phys. Med. Bio., vol. 53, pp. 2991-3006, 2008. [76] M. Chen, R.J. Jaszczak, J.E. Bowsher, D.R. Gilland, An evaluation of cardiac uniformity, contrast, and SNR with dual-head 180 and triple-head 360 SPECT scans, IEEE Trans. Nuc. Sci., vol. 48, no. 4, pp. 1428-1434, 2001. [77] S.D. Wollenweber, B.M.W. Tsui, D.S. Lalush, E.C. Frey, G.T. Gullberg, Evaluation of myocardial defect detecti on between parallel-hole and fan-beam SPECT using the Hotelling trace, IEEE Nuc. Sci. Symp. Conf. Rec. vol. 2, pp. 1428-1432, 1998.

PAGE 101

101 [78] C.K. Abbey, H.H. Barrett, M.P. Eckstein, P ractical issues and methodology in assesment of image quality using model observers, SPIE Med. Img. Symp., vol. 3032, pp. 182-194, 1997. [79] K. Fukunaga, R.R. Hayes, Effects of sample size in classifier design, IEEE Trans. Pattern Analysis Mach. Int., vol. 11, pp. 873-885, 1989. [80] R.F. Wagner, H.P. Chan, B. Sahiner, N. Petrick, J.T. Mossoba, Finite-sample effects and resampling plans: Applications to linear cl assifiers in computer-aided diagnosis, SPIE Img. Proc. Conf., vol. 3034, pp. 467-477, 1997. [81] M.V. Narayanan, H.C. Gifford, M.A. King, P.H. Pretorius, T.H. Farncombe, P. Bruyant, M.N. Wernick, Optimization of Iterative Reconstructions of 99mTc Cardiac SPECT Studies Using Numerical Observers, IEEE Trans. Nuc. Sci. vol. 49, no. 5, pp. 2355-2360, 2002. [82] H. Devries, The quantum char acter of light and its bearing upon threshold of vision, the differential sensitivity and visual acuity of the eye, Physica (Utrecht), vol. 10, pp. 553564, 1943. [83] A. Rose, The sensitivity performance of the human eye on an absolute scale, J. Opt. Soc. Am., vol. 38, pp. 196-208, 1948. [84] A. Rose, Television pickup tube s and the problem of vision, Adv. Electron., vol. 1, p. 131, 1948. [85] L.M. Biberman, ed., Perception of displa yed information, Plenum, New York, 1973.

PAGE 102

102 BIOGRAPHICAL SKETCH Jason Glenn Parker was born on October 26, 197 9, in Orlando, FL. He was the oldest of four children and grew up in Venice, FL. He ear ned an A.A. at Manatee Community College in 2002, a B.S. at Rensselaer Polytechnic Institut e in 2005, and a Ph.D. from the University of Florida in 2008. After graduating from UF he will begin working as a Senior Scientist in Dayton, OH, supervising research into advanced magnetic resonance imaging techniques.