FUNCTIONAL DYNAMIC SPECT IMAGING USING A SINGLE SLOW CAMERA ROTATION By Troy Farncombe B.Sc, University of Calgary, 1995 M.Sc., University of British Columbia, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF PHYSICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October 2000 © Troy Farncombe, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics The University of British Columbia 6224 Agricultural Road Vancouver, B.C., Canada V6T 1Z1 Date: A B S T R A C T Dynamic single photon emission computed tomography (SPECT) is a relatively new imaging method that uses radioactive tracers and tomographic data acquisition tech-niques in order to quantify temporal changes in regional radiotracer concentrations within a patient. This is important as the rate of change in tracer concentration within an organ often can be related to the functional ability of that organ. In this work, a new method is presented that is able to determine these kinetic rates while using a conventional single or multiple detector SPECT camera, and more importantly, a single, slow camera rotation in the data collection process (herein this reconstruction method will be referred to as dSPECT). This reconstruction method is based on the fact that a temporal change in the activity concentration at a given location can be represented by a linear inequality constraint over time. Two iterative reconstruc-tion algorithms, constrained least squares (CLS) and dynamic expectation maximization (dEM), have been tested using this approach with a variety of computer simulations and phantom experiments. In simulations involving a slow dynamic change of activity, results indicate that the dSPECT reconstruction procedure typically produces kinetic parameter estimates with a 7% error when using projection data acquired with a single 180° rotation of a triple headed SPECT camera system. This error increases to about 15% for data acquired with a dual head system and further increases to about 50% for single detector acquisitions. When simulated with faster dynamic parameters, errors increased slightly to about 8%, 12% and 55% for acquisitions involving triple, dual and single head systems respectively. ii T A B L E O F C O N T E N T S Abstract ii Table of Contents iii List of Figures viii List of Tables xiv Acknowledgments xvii 1 Outline of this Work 1 1.1 Thesis Organization 2 1.2 Thesis Conventions and Abbreviations 4 2 Introduction to Nuclear Medicine and Image Reconstruction 5 2.1 Nuclear Medicine 5 2.1.1 Radionuclides , 6 2.1.2 Clinical Applications of Nuclear Medicine 6 2.2 Equipment Used in Nuclear Medicine 11 2.2.1 The Gamma Camera 12 2.3 Image Reconstruction Procedures 16 2.3.1 Analytical Reconstruction Methods 17 2.3.2 Iterative Reconstruction Methods 18 2.4 Chapter Summary 25 3 Dynamic S P E C T Applications and Current Methods 26 3.1 Applications of Dynamic Functional Imaging 26 iii 3.1.1 Myocardial Metabolism of Fatty Acids 27 3.1.2 Myocardial Perfusion 28 3.1.3 Renal Function 30 3.1.4 Pharmacokinetic (PK) Modelling 31 3.1.5 Reduction of Dynamic Artifacts in SPECT Imaging 33 3.2 Dynamic Functional Imaging Methods 34 3.2.1 Planar Dynamic Imaging 34 3.2.2 Positron Emission Tomography 36 3.3 Current Dynamic SPECT Methods 38 3.3.1 Fast Rotation Image Based SPECT 39 3.3.2 Slow Rotation Dynamic SPECT Methods 42 3.3.3 Linear Methods 44 3.4 Chapter Summary 46 4 The d S P E C T Method 48 4.1 dSPECT Requirements 48 4.2 Theory of dSPECT Image Reconstruction 51 4.2.1 Linear Inequality Constraints 52 4.3 Components of the dSPECT Method 56 4.3.1 The Activity Vector 56 4.3.2 The Probability System Matrix . . . .' 57 4.3.3 The Data Vector 59 4.3.4 The Linear Inequality Difference Tensor 59 4.4 The dSPECT Reconstruction 63 4.4.1 Determine Object Bounds 63 4.4.2 Determination of Peak Activity Position 65 iv 4.4.3 Determination of Dynamic and Static Voxels 67 4.4.4 Initial Object Estimate 70 4.4.5 Additional Shape Constraints 71 4.4.6 Specify Regularizer 72 4.4.7 Dynamic Image Reconstruction 73 4.4.8 a posteriori Fit 78 4.5 Chapter Summary 80 5 Computer Simulations 81 5.1 Simulation Methods 81 5.2 Data Analysis Methods 84 5.2.1 Relative Standard Deviation . . 84 5.2.2 Kinetic Parameter Data Analysis 86 5.2.3 Peak Activity Position Error Analysis . 86 5.3 Analytical Simulations 87 5.3.1 Mono Exponential Washout 88 5.3.2 Dual Exponential Washout 91 5.3.3 Dual Exponential with Slow Washin and Washout 94 5.3.4 Dual Exponential with Fast Washin and Washout 95 5.4 Discretized Simulations 97 5.4.1 Mono Exponential Washout 98 5.4.2 Dual Exponential Washout 100 5.4.3 Dual Exponential with Washin and Washout 102 5.5 Chapter Summary 104 6 Phantom and Volunteer Experiments 105 v 6.1 The Dynamic Heart-in-Thorax Phantom 105 6.1.1 Dynamic Heart Insert 107 6.2 Reconstruction Accuracy Assessment Methods 110 6.3 Dynamic Phantom Experiments I l l 6.3.1 Non-attenuated Experiments I l l 6.3.2 Attenuated Experiments 115 6.4 Dynamic Renal Experiments 118 6.4.1 MAG3 Experiments 118 6.5 Chapter Summary 120 7 Conclusion 122 7.1 General Comments 122 7.2 Future Work 123 7.2.1 Fully Four Dimensional Reconstruction 123 7.2.2 Photon Scatter Compensation . : 124 7.2.3 Myocardial Gating 124 A Likelihood Functions 125 A . l Gaussian Likelihood Objective Function 125 A.2 Poisson Likelihood Objective Function 127 A. 3 The Static Expectation Maximization (EM) Algorithm 128 B d S P E C T Reconstructions 130 B. l Analytical Simulation Results 130 B . l . l Mono Exponential Washout 130 B.l.2 Dual Exponential Washout 138 B.l.3 Dual Exponential with Slow Washin and Washout 145 vi B.1.4 Dual Exponential with Fast Washin and Washout 152 B.2 Discretized Simulation Results 159 B.2.1 Monoexponential Washout 159 B.2.2 Dual Exponential Washout 166 B.2.3 Dual Exponential with Washin and Washout 173 B.3 Experimental Results 180 B.3.1 Nonattenuated Experiments 180 B.3.2 Attenuated Experiments 184 Bibliography 187 vii L I S T O F F I G U R E S 2.1 A schematic diagram of a typical SPECT scintillating detector SPECT camera 12 2.2 A schematic diagram of a parallel hole collimator 13 2.3 A photomultiplier position logic circuit connected to the back side of the scintillating detector. The amount of light produced in each quadrant ( X + , X - , Y + and Y~) are combined in the summing matrix circuit (SMC) to obtain the X and Y positions on the detector surface. The total light output is given by Z 15 2.4 The reordering of projection data into a sinogram for use in tomographic reconstruction 18 2.5 The steps involved in a filtered backprojection reconstruction of tomo-graphic data. The original projection data is first ID Fourier transformed and then multiplied by appropriate filters prior to being inverse trans-formed. This modified projection data is then backprojected at each angle in order to yield the transaxial image. 19 2.6 The reordering of the object pixel indices into a vector 20 2.7 Examples of parallel hole system matrices acquired with different size cam-era bins. . 21 2.8 An example of the geometrical contribution of an 8x8 voxel object to a single camera bin. For each camera bin, the resultant matrix gets reordered into a single row of the system matrix C. (Values shown are approximate.) 22 viii 2.9 System matrices modelling the case of a) collimator blurring in the transax-ial direction only (2D), and b) true camera response with blurring in all directions (3D). 23 3.1 A simple one compartment model used in PK modelling. Following in-jection of the drug into the compartment, the amount of drug within the compartment (A(t)), is reduced in an exponentially decreasing manner with decay rate k 32 3.2 Shown on the right is a reconstructed image using the standard filtered backprojection routine of an object undergoing a changing activity distri-bution. On the left is the true object. Notice the lack of definition in the lung regions and the streaks emanating from the leftmost heart region in the reconstructed image 34 3.3 Examples of a dynamic planar image sequence of a renogram acquired at 10 second intervals. Images shown correspond to the times 1 min, 5 min, 10 min and 15 min after radiotracer injection. On the right is a plot depicting the amount of activity present in each kidney as a function of time (a time activity curve) 35 3.4 Examples of kinetic models used in dynamic positron emission tomography (PET) imaging 38 4.1 Diagram depicting similar effects yielded in projection data due to a chang-ing tracer distribution or because of attenuation effects with a static dis-tribution. In a) a dynamic object is present in the absence of attenuation, while b) represents a static object with attenuation. In both cases, 180° projection measurements were acquired in a counterclockwise rotation with the camera starting at the 6 o'clock position 49 ix 4.2 Schematic representation of the Siemens Profile transmission imaging ge-ometry. Transmission imaging is performed simultaneously with emission imaging through the use of radioisotopes with different emission energies. 51 4.3 A simplified representation of an expected myocardial time activity curve for fatty acids 53 4.4 The reordering of object voxels from a single transaxial slice acquired at 16 time frames into a spatio-temporal vector for use in the dSPECT image reconstruction process. Notice voxel indices proceed over time first, followed by spatially 57 4.5 A flowchart depicting the various steps involved in the dSPECT recon-struction method 64 4.6 a) The true time activity curve (TAC) of a single voxel for a dynamic increasing/decreasing simulation, b) The reconstructed TAC for the sim-ulation assuming only increasing, and only decreasing activity 67 4.7 An example of the dynamic mask used in dSPECT reconstructions. A value of -1 signifies a null voxel, values of 0 correspond to static voxels, while other values correspond to the positions of peak voxel activity. . . . 70 4.8 An example of the initial object activity values used in a dSPECT recon-struction. On the left is the static EM reconstruction shown and on the right are the initial time activity curves for the two heart sections for the dynamic mask shown in Figure 4.7 71 4.9 Examples of more stringent inequality constraints to produce increasing and decreasing activities. On the right is depicted an example of a com-bined increasing/decreasing activity with convexity and concavity con-straints 73 x 4.10 The three possibilites of the peak activity position when used with the modified difference matrix above 80 5.1 A diagram of the annulus phantom used in the first set of simulations. Projection data was collected into 64 transaxial camera bins over a 20 minute acquisition time with various camera geometries 82 5.2 Left) Activity distribution map shown at t=5 min. Right) Nonuniform attenuation map phantom used in the attenuated simulations. Attenuation coefficients shown are in units of m m - 1 83 5.3 True annulus object shown at t = 5 min for mono exponential decay sim-ulation 89 5.4 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation 92 5.5 True image at t = 5 min along with time activity curves for each dynamic region for the slow dual exponential simulation 95 5.6 True image at t = 5 min along with time activity curves for each dynamic region for the rapid dual exponential simulation 96 5.7 True image at t = 5 min along with time activity curves for each dynamic region for the mono exponential decay simulation 99 5.8 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation. . 101 5.9 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential increase and decrease simulation. . . . . . 103 6.1 The dynamic heart-in-thorax phantom shown with the surrounding het-erogeneous attenuating material. The dynamic heart insert can be seen located centrally in the thorax 106 xi 6.2 Cutaway diagram of the mixing apparatus used to ensure uniform mixing of the inflowing water with the radioactivity 108 6.3 A schematic of the experimental configuration of the phantom and the ex-ternally situated syringe pump and driving motor. Arrangement of equip-ment is shown from the top 109 6.4 An attenuation map of the dynamic phantom thorax acquired with the Siemens Profile transmission system. Attenuation coefficients are in units of c m - 1 115 6.5 Top) Images acquired using a dynamic planar imaging technique. Bottom) Reprojected dSPECT images into planar images 120 6.6 Images of transaxial slices reconstructed using the dEM reconstruction algorithm on data acquired with a single 90° acquisition protocol. The plateau seen in the TA curve for the right kidney may represent a second point of increasing activity due perhaps to a second pass of radiophar-maceutical within the kidney. As the reconstruction algorithm can only reconstruct TA curves with a single peak position at this time, the second peak appears as a plateau ; 121 B . l True annulus object shown at t = 5 min for mono exponential decay sim-ulation 131 B.2 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation 138 B.3 True image at t = 5 min along with time activity curves for each dynamic region for the slow dual exponential simulation 145 B.4 True image at t = 5 min along with time activity curves for each dynamic region for the rapid dual exponential simulation 152 xii B.5 True image at t = 5 min along with time activity curves for each dynamic region for the mono exponential decay simulation 159 B.6 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation 166 B.7 True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential increase and decrease simulation 173 xiii L I S T O F T A B L E S 2.1 A summary of the common radionuclides used in nuclear medicine imaging. 6 2.2 Table of common applications of nuclear medicine imaging 7 2.3 Common scinitillating materials used in nuclear medicine imaging. Nal is used, primarily for low energy photons, while BGO and B a F 2 are used for high energy (511 keV) photons as their higher densities results in higher stopping power of the higher energy photons. Note the two light photons produced in the scintillation of BaF 2 14 2.4 Examples of some of the regularizing functions used in least squares image reconstruction methods 24 3.1 Examples of possible dynamic processes in the body and current radio-pharmaceuticals that may be used in dynamic imaging 27 5.1 Summary of simulation acquisition parameters used for all simulations. . 88 5.2 True annulus phantom dynamic parameters used in simulations of mono exponential decay 89 5.3 Summary of overall accuracy values for the mono exponential decreasing activity simulations 90 5.4 True dynamic parameters used in simulations involving dual exponential decay in the annulus phantom 93 5.5 Summary of overall accuracy values for dual exponential decreasing simu-lations using the annulus phantom 93 5.6 True dynamic parameters used in simulations involving a combined in-creasing and decreasing dual exponential function 94 xiv 5.7 True dynamic parameters used in simulations involving a combined in-creasing and decreasing dual exponential function with fast dynamic com-ponents 95 5.8 Summary of overall accuracy values for slow washin / washout simulations. 97 5.9 Summary of overall accuracy values for fast washin / washout simulations. 97 5.10 Summary of the acquisition protocols used in discretized thorax simulations. 98 5.11 Summary of the true dynamic parameters used in mono exponential decay simulations of the thorax phantom 98 5.12 Summary of overall accuracy values for mono exponential decreasing tho-rax simulations. (Accuracy values appear large due to large errors in the low count lung regions.) 100 5.13 Summary of the true dynamic parameters used in dual exponential decay simulations of the thorax phantom. 101 5.14 Summary of overall accuracy values for dual exponential decreasing thorax simulations. 102 5.15 Summary of the true dynamic parameters used in dual exponential washin/washout simulations of the thorax phantom 102 5.16 Summary of overall accuracy values for dual exponential increasing and decreasing thorax simulations 103 6.1 Summary of experimental acquisition protocols used for experiments of the dynamic heart phantom with no attenuating material. 112 6.2 True projection time activity curve parameters for a mono exponential washout of activity acquired with a single head camera. 112 6.3 True projection time activity curve parameters for a mono exponential washout of activity acquired with a triple head camera. 113 xv 6.4 True projection time activity curve parameters for a dual exponential change of activity acquired with a dual head camera. 113 6.5 Reconstructed dynamic paramters of non-attenuated, single head dynamic phantom acquisition 114 6.6 Reconstructed dynamic paramters of non-attenuated, triple head dynamic phantom acquisition 114 6.7 Reconstructed dynamic paramters of non-attenuated, dual head dynamic phantom acquisition 115 6.8 Summary of the acquisition protocols for experiments involving attenua-tion and a combined increase/decrease 116 6.9 Preset washout times for experiments involving both increase/decrease be-haviour, and nonuniform attenuation. Half-life times have been calculated based upon the preset water flow rates 117 xvi A C K N O W L E D G M E N T S This dissertation would not have been possible without the assistance of numerous par-ties. I would first like to acknowledge the financial support that I have received from the Natural Sciences and Engineering Research Council of Canada (NSERC). Additional support has been provided by the Government of France in the means of travel assistance for travelling to France to work with collaborators on this project. I thank both of these organizations for their financial assistance in the completion of this project. On a personal basis, I first have to thank my supervisor, Dr. Anna Celler for all the ideas and support that she has provided over the past several years. It has been great to work with a person who truly believes in the research being undertaken and has put together such a strong group to achieve a common goal. Second, this thesis has been aided greatly from the wisdom of Dr. Ronald Harrop. His insight into not only this research project, but also into matters unrelated to this work has made the preparation of this thesis possible. Additionally the help given by collaborators Dr. Jean Maeght and Dr. Dominik Noll from the Universite Paul Sabatier (Toulouse, France) has been instrumental in the completion of this work, and I thank them sincerely for all their thoughts and ideas. Many of the simulations performed in this work have been provided by Caitlin Bever and I thank her for her help in this, as well as being a good friend around the office. Finally I would like to thank my parents George and Sheila Farncombe for all the support they have given me over the years leading to the completion of this degree. xvii C H A P T E R 1 O U T L I N E O F T H I S W O R K In this work, a new image reconstruction method for nuclear medicine imaging is pre-sented that is able to reconstruct a time series of three dimensional images of a changing radiotracer distribution within the body. This method, known as dynamic single photon emission computed tomography (dSPECT), is unique in the fact that the reconstruction procedure uses projection data acquired using a single slow rotation of a conventional SPECT camera. The advantages of this method are as follows. • It allows dynamic SPECT imaging to be performed using any commonly available SPECT camera (ie., single, dual, and triple head). • It improves upon the low signal to noise ratio (poor contrast) of other dynamic SPECT methods using a fast rotation acquisition. • It is inherently quantitative as it is an iterative method and includes patient specific attenuation correction and accurate geometrical modelling. • It allows for easy clinical implementation using current computer hardware. 1 Chapter 1. Outline of this Work 2 1.1 T H E S I S O R G A N I Z A T I O N In the presentation of this reconstruction method, the work is divided into the following chapters: C H A P T E R 2 - INTRODUCTION TO N U C L E A R M E D I C I N E AND I M A G E RE C O N S T R U C T I O N In this chapter, a short overview of nuclear medicine imaging is presented. A description of some common clinical procedures is given along with information on the equipment used for such studies. An introduction to conventional SPECT image reconstruction follows where the object is to determine the three dimensional spatial distribution of a radiopharmaceutical within the body by means of a series of two dimensional projection measurements acquired at various angles around the patient. C H A P T E R 3 - D Y N A M I C S P E C T APPLICATIONS AND C U R R E N T M E T H O D S Chapter 3 presents some potential clinical applications of dynamic SPECT imaging along with a description of some of the current methods used to obtain dynamic information related to organ function from these studies. C H A P T E R 4 - T H E D S P E C T M E T H O D In this next chapter, the dSPECT image reconstruction method is described in detail. Stated briefly, the dSPECT method represents the change in activity concentration over time within each pixel by linear inequalities. The subsequent dSPECT image recon-struction procedure then uses these inequality relations in a mathematical optimization routine. Here the general dSPECT image reconstruction method will be described along with implementation details which aid in either reconstruction accuracy or speed. Chapter 1. Outline of this Work 3 C H A P T E R 5 - C O M P U T E R S I M U L A T I O N S To test the dSPECT image reconstruction process described in the previous chapter, a series of computer simulations is presented using two different dynamic objects modelling different dynamic behaviours and acquired with various camera acquisition geometries. This chapter addresses the accuracy of the dSPECT method by using a series of quanti-tative accuracy measurements that will be described. C H A P T E R 6 - P H A N T O M A N D V O L U N T E E R E X P E R I M E N T S As a continuation to the simulations performed in the previous chapter, this chapter further validates the dSPECT method by using experimentally obtained dynamic data created using a dynamic heart-in-thorax phantom. Again, accuracy measurements will be given for these dynamic reconstructions. An initial trial was also performed on a healthy volunteer subject at this time to evaluate the clinical utility of dSPECT imaging. C H A P T E R 7 - G E N E R A L C O M M E N T S A N D C O N C L U S I O N S The thesis concludes with general comments on the dSPECT reconstruction method, and the presentation of some future work to be conducted in order to further validate and improve the method. In addition to the seven chapters described, two appendices are also included. They consist of: Chapter 1. Outline of this Work 4 A P P E N D I X A - A D E R I V A T I O N O F L I K E L I H O O D F U N C T I O N S A N D T H E E X P E C T A T I O N M A X I M I Z A T I O N A L G O R I T H M A derivation of the two different likelihood functions (Gaussian and Poisson) used in the dSPECT reconstruction method is given. As well, starting from the Poisson like-lihood function, a well known standard iterative reconstruction algorithm (expectation maximization, or EM) is derived. A P P E N D I X B - S I M U L A T I O N R E S U L T S Reconstructed images and accuracy measurements are presented for a sampling of camera acquisitions and dynamic behaviours. 1.2 T H E S I S C O N V E N T I O N S A N D A B B R E V I A T I O N S Throughout this thesis various symbols and shorthand notation are used. These include extensive use of the following. i - object pixel index j - camera bin index k - camera projection index N - total number of object voxels S - number of projection angles C - static projection operator matrix C - dynamic projection operator matrix x - static voxel values x - voxel differences Ai - difference matrix for voxel i A - difference tensor d - measured projection data £R - regional error value ST - total error value C H A P T E R 2 INTRODUCTION TO NUCLEAR MEDICINE AND IMAGE RECONSTRUCTION The following chapter presents a brief synopsis of nuclear medicine imaging and dis-cusses the important role it plays in diagnostic medicine. Examples of nuclear medicine imaging procedures will be presented along with details of how three dimensional tracer distribution information is created and displayed for interpretation. In the latter half of this introduction, the basic three dimensional reconstruction problem will be outlined in mathematical terms and techniques will be presented for reconstructing the underlying three dimensional activity distribution from a series of two dimensional projection data. 2.1 N U C L E A R M E D I C I N E Nuclear medicine imaging is often regarded as a functional imaging technique because resultant images are able to depict information related to the organ function rather than only corresponding to anatomy. In order to do this, pharmaceuticals are labelled with a traceable radioactive element and introduced into a patient. When this radioactive element decays, the emitted radiation (usually 7 photons), is detected by an external 5 Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 6 detector. Since the location of the radioactive decay can be found through image recon-struction techniques, it is then possible to know where the pharmaceutical travelled in the body, and hence infer information related to the physiological condition of the body. 2.1.1 R A D I O N U C L I D E S Most radionuclides used in nuclear medicine emit low energy (« 70-511 keV) gamma radiation and have relatively short physical half-lives. The most common radionuclide used in nuclear medicine is 9 9 m T c due to its short half-life and ease of production from a " M o generator. Summarized in Table 2.1 are some of the more common radionuclides used for nuclear medicine imaging. Single Photon Emitters Positron Emitters Nuclide Half-life Nuclide Half-life 99m rpc 6.02 hrs n C 20.3 min 123j 13.3 hrs i s N 10.0 min m I n 67.2 hrs i s 0 2.0 min 1 3 3 X e 125.8 hrs 18p 110 min 201rp| 72.9 hrs 8 2 Rb 1.27 min 6 7 G a 78.1 hrs Table 2.1: A summary of the common radionuclides used in nuclear medicine imaging. 2.1.2 C L I N I C A L A P P L I C A T I O N S O F N U C L E A R M E D I C I N E Almost every part of the body has been imaged to some degree in nuclear medicine [1]. Table 2.2 presents some of the more common radiopharmaceuticals used in nuclear medicine along with the typical clinical applications of these pharmaceuticals. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 7 Application Area Pharmaceutical Purpose Central Nervous System 1 2 3 I - Isopropyl-p-iodoamphetamine (IMP) Cerebral perfusion and metabolism (decreased uptake in infarcted tissue) 99mTc-Hexamethylpropylene Amine Oxime (HMPAO) Cerebral perfusion 9 9 m Ethy l Cyteinate Dimer (ECD) Cerebral perfusion 18F-Fluorodeoxyglucose (FDG) Glucose metabolism Thyroid 123I-Sodium Iodide Thyroid uptake 9 9 m T c pertechnetate Thyroid size Respiratory 1 3 3 X e and 1 2 7 X e gas Detection of airway obstructions 9 9 m T c DTPA aerosol Ventilation 9 9 m T c macroaggregated albumin (MAA) Lung perfusion Heptobiliary 9 9 m T c albumin colloid Tumour detection 6 7 G a gallium citrate Detection of abscesses and tumours 9 9 m Tc-IDA derivatives (HIDA, PIPIDA, etc) Gallbladder imaging Renal 99mTc-Mercaptoacetyl-glycylglycylglycine (MAG3) Renal plasma flow 9 9 m Tc-DTPA Glomerular filtration Skeletal 99mTc-phosphate Bone blood flow and bone formation rate Myocardial 9 9 mTc-Sestamibi Blood perfusion 99mTc-Teboroxime Perfusion 2 U iTl-Thallous Chloride Perfusion Table 2.2: Table of common applications of nuclear medicine imaging. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 8 M Y O C A R D I A L I M A G I N G One of the most common applications of nuclear medicine imaging is in the assessment of myocardial perfusion. The functional ability of the myocardium can be related to the amount of blood supplied to it in the coronary arteries [2]. If the flow of blood is impaired, then the heart works less efficiently and may have difficulty in supplying the blood to the extremities and other vital organs. Typical radiopharmaceuticals used for myocardial imaging (eg., 9 9 m T c sestamibi and 201thallous chloride) are injected into a patient and travel throughout the body in the bloodstream. A small percentage of the injected pharmaceutical will be taken up by regions of the heart that are perfused with blood. Once in the heart, the pharmaceutical will remain there for an extended period of time thus permitting imaging to take place. By interpreting the resultant images produced through this procedure, regions with low blood perfusion can be differentiated from regions with higher perfusion based upon the relative amounts of tracer present within these regions [3] [4] [5] [6]. B R A I N I M A G I N G Another common application of nuclear medicine imaging is that of brain imaging. One of the most common forms of brain imaging uses positron emission tomography (PET) with the radiopharmaceutical 18F-fluorodeoxyglucose (FDG). This pharmaceutical is a glucose analog that travels in the bloodstream and is taken up by cells in the body undergoing glucose metabolism. Once taken up by these cells, the FDG molecule gets phosphorylated into various metabolic products and one of these products, FDG-6-phosphate, is retained for a long period of time within the cells. Because of this, it is possible to image the location of the attached 1 8 F and thus infer the location of the glucose metabolism in the body. Regions of the body with increased glucose metabolism (eg., stimulated brain Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 9 centres, tumours) will have higher amounts of FDG-6-phosphate present and so will appear as bright regions in the subsequent images. Other radiopharmaceuticals also used for brain imaging are the single photon emitting pharmaceuticals 123I-isopropyl-p-iodoamphetamine (IMP), 99mTc-hexamethylpropylene amine oxime (HMPAO) and 9 9 mTc-ethyl cyteinate dimer (ECD). In regions of the brain with less perfusion (due perhaps to stroke, infarction, etc.), the uptake of these pharma-ceuticals will be diminished relative to normally perfused regions [7] [8]. L U N G I M A G I N G Another study routinely performed in many nuclear medicine departments is that of lung imaging. This type of imaging can further be broken down into two subtypes; perfusion and ventilation. In perfusion imaging, 9 9 m T c macroaggregated albumin (MAA) is injected into the patient and is trapped in some of the peripheral pulmonary arterioles and capillaries thus reflecting regions of the lungs perfused with blood at the time of injection. Thus, from this scan it may be possible to visualize deficiencies in the blood supply to the lungs due to pulmonary emboli (PE). In order to more accurately diagnose a pulmonary embolism, usually a lung ventilation scan using a 9 9 m T c labelled aerosol is performed in conjuction with the perfusion scan. This scan is useful in itself to depict any airway obstructions since the gas is only capable of reaching alveoli that are not obstructed. By performing both a perfusion and a ventilation scan (termed a V / Q scan) it is possible to diagnose a PE since a normal ventilation scan combined with a decreased perfusion scan may often indicate the presence of a PE [9] [10]. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 10 GI T R A C T I M A G I N G Organs such as the liver, spleen, stomach, and intestinal tract are also routinely imaged in nuclear medicine. Through the use of 9 9 m T c albumin colloid, it is possible to evaluate the effectiveness of the liver's reticuloendothelial tissue at removing particulate matter from the bloodstream. Compromises in this tissue are reflected by abnormal uptake patterns within the liver[ll]. Other GI imaging procedures include the measurement of gastric accomodation [12], and bleeding within the intestinal tract [13]. B O N E I M A G I N G Another of the more commonly performed diagnostic procedures used in nuclear medicine is bone imaging. Often small bone fissures are not possible to visualize clearly using more conventional imaging techniques (eg, x-rays). In these cases, nuclear medicine imaging using phosphate and diphosphonate complexes (eg., methylene diphosphonate, or MDP) can often be used to visualize the location and severity of these fissures [14]. This is because these substances are used extensively in the formation of new bone material and so regions in the body undergoing active bone formation will utilize these substances to a greater extent than other regions. Therefore by imaging with these agents, areas of abnormal bone growth can be visualized. Additionally bone tumour metastases can also be clearly visualized due to the increased concentration of these pharmaceuticals [15]. O N C O L O G I C A L I M A G I N G Cancer is the rapid, uncontrolled growth of cells. Owing to their fast growth, the intake and turnover of physiological materials is quite high. This fact can be exploited in nuclear medicine imaging in order to determine the location and classification of cancerous lesions. If a compound such as glucose is labelled with a radionuclide and injected into a patient Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 11 with a cancerous tumour, the material will be taken up preferentially by those cells that are undergoing mitosis (ie, the cancerous cells). As a result, the local concentration of this material would be increased in the vicinity of such cancerous lesions and hence by detecting the emitted radiation from the radionuclide, it is possible to determine the location of such tumours [16] [17]. R E N A L I M A G I N G Another common application of nuclear medicine is in the determination of renal function. An impairment of blood flow through the kidneys limits their effectiveness to remove waste from the body. By using a blood perfusion imaging agent that gets trapped in renal arterioles and capillaries, it is possible to determine regions of the kidneys with impaired blood flow. Another measure of renal function measures the rate at which blood is filtered as there is a direct connection between renal filtration and renal function. As in the case of most other nuclear medicine procedures, in renal imaging, a radiolabeled pharmaceutical is injected into the blood. Some of this pharmaceutical will travel to the kidneys where it gets filtered and eventually passed out of the body in the urine. By measuring the rate at which the tracer enters and exits the kidneys, a good indication can be obtained as to the functional state of the kidneys [18] [19]. 2.2 E Q U I P M E N T U S E D IN N U C L E A R M E D I C I N E With the exception of the positron emission applications (eg., FDG brain imaging), the equipment used for the nuclear medicine procedures described consists of the gamma camera that detects the photons emitted in the radioactive decay and the associated computers that either control the movement of the camera or that display and analyze Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 12 the collected data. 2.2.1 T H E G A M M A C A M E R A The gamma camera (shown in Figure 2.1) is the main piece of hardware used in nu-clear medicine imaging and is best discussed in the terms of its individual components. Proceeding from the exterior of the camera closest to the patient these parts consist of: Display Terminal Acquisition Computer Pulse Height Analyzer ' Photomultipliers ^ Detector Crystal Collimator Figure 2.1: A schematic diagram of a typical SPECT scintillating detector SPECT cam-era. 1) Lead collimator (Figure 2.2) - The purpose of the lead collimator on the gamma camera is to limit the incidence angle of the photons being detected by the detector. Collimators used in nuclear medicine imaging include fanbeam, pinhole and conebeam, but the most commonly used collimator is the parallel hole collimator. This type of collimator typically consists of a stack of 0.2 mm thick lead foils pressed together in order to create a series of hexagonal holes. These holes allow only photons travelling Position Logic Circuits Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 13 at small angles perpendicular to the camera face to pass through and be incident upon the camera surface. By limiting the acceptance angle of the incident photons, a priori knowledge regarding the originating position of the incoming photons is available. Figure 2.2: A schematic diagram of a parallel hole collimator. The size of the holes varies depending on the collimator used, but in a typical high resolution, low energy collimator, the holes would be approximately 1.4 mm in diameter with the 0.2 mm lead foil thickness separating them (called the septa). The overall thickness for such a collimator would typically be about 2.5 cm. This combination of hole size, collimator thickness and septal thickness allows for an efficiency of approximately 1 in 104 photons emitted from the patient being able to pass through the collimator and be detected. Collimators with larger holes would obviously have a higher photon detection efficiency, yet would also allow photons emitted from larger angles to pass through thus resulting in a blurring of the source distribution and thus lowering the detector spatial resolution. Conversely, collimators with smaller holes have lower detection efficiency, but increased spatial resolution. 2) Detector Crystal - Once the photons have passed through the collimator, it is necessary to detect them. Depending on the collimator used, a limited number of photons Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 14 emitted from the patient will pass through the collimator and become incident upon the detector surface. Because of the relatively low efficiency of the collimators, it is important to utilize an efficient photon detection system. The detector used most often in gamma cameras is the Nal scintillating crystal. The Nal scintillator is a highly efficient photon detector that converts an incoming gamma photon into visible light photons. Nal detectors have a high detection efficiency and as seen in Table 2.3, have a high light output compared to other scintillating detec-tors. However, Nal crystals suffer from a relatively poor energy resolution and because of light spread within the crystal itself, they also have a relatively poor intrinsic spatial resolution. Material Density Relative Light Decay Wavelength (g/cm3) Output Constant (ns) of Emission (nm) Nal 3.67 1.00 230 410 Bismuth germanate 7.13 0.12 300 480 (BGO) Barium fluoride 4.89 0.05 0.6-0.8 220 (BaF2) 0.15 630 310 Table 2.3: Common scinitillating materials used in nuclear medicine imaging. Nal is used primarily for low energy photons, while BGO and B a F 2 are used for high energy (511 keV) photons as their higher densities results in higher stopping power of the higher energy photons. Note the two light photons produced in the scintillation of BaF 2 . 3) Photomultiplier Array (Figure 2.3) - While it has been seen that Nal is a highly efficient detector with high light output, the amount of light produced within these detectors is still quite small. In order to detect this light, photomultipliers are used on the back of the crystal to convert the light signal produced by the scintillator into an amplified electrical signal. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 15 Z Z Figure 2.3: A photomultiplier position logic circuit connected to the back side of the scintillating detector. The amount of light produced in each quadrant ( X + , X - , Y + and Y~) are combined in the summing matrix circuit (SMC) to obtain the X and Y positions on the detector surface. The total light output is given by Z. Additionally, the photomultiplier array serves to determine the position of the de-tected event on the camera surface. That is, it localizes the position of each detected event into a discretized grid of (X,Y) camera locations (henceforth referred to as camera bins). This system works by combining the electrical signal from each photomultiplier in order to determine a resultant (X,Y) position of the detected event. When using a parallel hole collimator, only photons travelling perpendicular to the camera surface are detected by the scintillating crystal. As well, the location of this event on the camera surface is determined from the position circuitry. This type of information represents the projection of a three dimensional radiotracer source distribution into a two dimensional image much in the same way that an x-ray radiograph is a two dimensional representation of a three dimensional object. From this point on, information of this type Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 16 will be referred to as projection data. 4) The Acquisition Computer - Controlling the data collection and the motion of the gamma camera is an acquisition computer. Through the use of this computer, all data acquisition parameters such as the number of detector bins, the size of the detector bins, and the number and duration of projection angles are all specified. Additionally; this computer controls other camera specific details such as photomultiplier tube gain. 5) The Image Reconstruction/Display Computer - The image display computer takes the projection data that has been acquired by the acquisition computer and in the case of tomographic imaging, performs an image reconstruction procedure specified by the user in order to create an image set that can be interpreted by the viewer. In the case of SPECT image reconstruction, the final image set consists of a series of adjacent transaxial slices through the patient. These transaxial images can either be viewed independently as 2D images, or stacked and rendered using visualization tech-niques into a fully 3D volume data set. Various image manipulation tools can further be applied to this data. A common option performed at this stage is to apply various filtering kernels to the transaxial images in order to smooth out inherent noise from the reconstruction and data acquisition procedures. As well, displaying the images us-ing different color sequences allows the user to pick out details that may ordinarily be overlooked. 2 . 3 I M A G E R E C O N S T R U C T I O N P R O C E D U R E S When data is acquired by a SPECT camera, it consists of a collection of two dimensional projection measurements of the radionuclide source distribution at various angles. While this type of two dimensional imaging is useful for some studies, the purpose of tomo-graphic imaging is to produce a three dimensional representation of the interior of an Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 17 object. The method of determining the interior of an object by means of a collection of projection measurements at various angles around an object was first proposed by Radon in 1917. In this original work, an analytical formulation was derived that stated if an infinite number of projections were acquired around an object, one could determine what was in the interior of the object. This is known as a Radon Transform. In the simplest case of image reconstruction, the data for each transaxial slice through the object is reconstructed independently into a transaxial image. In order to do this, projection data corresponding only to the current slice under study must first be isolated from the array of projection data. This involves extracting the projection data for all the camera bins corresponding to one transaxial slice and reordering it in the form of a two dimensional array representing the number of counts collected in each camera bin at each projection angle (see Figure 2.4). When data is arranged in this manner it is referred to as a sinogram because a single position in space gets mapped into a sinosoid in the projection space. 2.3.1 A N A L Y T I C A L R E C O N S T R U C T I O N M E T H O D S In the original work by Radon, it was necessary to collect an infinite number of projection measurements around the object in order to accurately reconstruct the interior. However, once this is done, it is possible to formulate an analytical solution to the reconstruction problem. One method to determine the three dimensional distribution of activity is with the filtered backprojection method. In the filtered backprojection method of image reconstruction, projection data for each transaxial slice at each projection angle first gets ID transformed into the Fourier domain where it is filtered in order to reduce the effects of high frequency noise and Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 18 Figure 2.4: The reordering of projection data into a sinogram for use in tomographic reconstruction. to account for a change of coordinate space from Cartesian to polar coordinates. This modified data then gets inverse ID transformed back into the spatial domain and a process known as backprojection is applied to it. This step essentially consists of smearing the number of photons collected in each camera bin at each projection angle back into the positions in space that may have emitted the photons. The steps involved in the filtered backprojection method are shown in Figure 2.5. For more details on this and other analytical reconstruction methods, the interested reader is directed to [20] [21] [22] [23]. 2.3.2 I T E R A T I V E R E C O N S T R U C T I O N M E T H O D S An alternative to filtered backprojection reconstruction are the methods referred to as iterative types. With these methods, a series of estimates is made as to the transaxial radionuclide distribution based upon the projection data that has been collected. These estimates are constantly altered and reevaluated throughout the reconstruction process Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 19 Figure 2.5: The steps involved in a filtered backprojection reconstruction of tomographic data. The original projection data is first I D Fourier transformed and then multiplied by appropriate filters prior to being inverse transformed. This modified projection data is then backprojected at each angle in order to yield the transaxial image. and eventually an estimate is obtained which best matches the observed data and also takes into account any a priori information known about the object. As mentioned previously, measurements acquired with a S P E C T camera are simply equal to two dimensional projections of a three dimensional object. If the locations within the object are indexed consecutively (Figure 2.6), then the projection can be defined by the linear operation as the matrix-vector product, C x = d (2.1) where x is a vector representing the amount of activity present in each object location, d is another vector containing the measured data from the gamma camera, and the matrix C is the system matrix which maps the transaxial object into the projection space. The system matrix represents the contributions that each small volume in space (herein referred to as volume elements, or voxels) makes to each camera bin. That is, it can be thought of as a probability matrix which represents the probability of detecting a Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 20 Figure 2.6: The reordering of the object pixel indices into a vector. photon emitted from a selected object voxel by a specified camera bin, and in the simplest case, is equal to the volume of each object voxel intersected by each camera bin. Two effects are present in this matrix; i) geometrical effects due to the collimator and camera acquisition parameters, and ii) object dependent effects such as scatter and attenuation. For now, only the geometrical effects will be discussed. The purpose of the camera collimator is to limit the angle from which photons emitted from the object are detected by the scintillating detector. In the case of a perfect parallel hole collimator, it is assumed that only photons travelling perpendicular to the camera surface will be incident on the detector. Furthermore, values in the system matrix depend upon parameters used in the data acquisition process. For example, if data were collected with only a few large camera bins, then each camera bin would contain contributions from many object voxels, that is, the system matrix would map a large number of object voxels into each detector bin (Figure 2.7). Conversely, if the data is acquired with smaller camera bins, the system matrix would map fewer voxels into each of the camera bins. As a numerical example consider the system matrix shown below for a 8x8 voxel object (Figure 2.8). For the one camera bin depicted, the contribution from each object voxel is shown by the numerical value present in the matrix. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 21 a) b) -WLiIIF - I -• • • • I IB M M ' J l l i k I i MP" M I 1 J W i l l Figure 2.7: Examples of parallel hole system matrices acquired with different size camera bins. When all the probabilities are computed for a given geometry, the elements can be rearranged into the form of an array where the CV^ elements of this matrix correspond to the probability of detecting a photon emitted from the ith. voxel in the jth camera bin at projection angle k. As discussed, the simplest system matrix assumes perfect parallel hole collimation, however, more complicated acquisition models can be used as well. A true collimator does not have perfect parallel hole collimation but rather has a limited acceptance angle from which incoming photons originated. Such a model can be incorporated into the sytem matrix by modelling the true shape of the camera response to each object voxel. If one considers the case of independent transaxial slices, then each camera bin will contain data orginating from a trapezoid as shown in Figure 2.9a). The exact size and shape of the trapezoid depends on the camera acquisition parameters and on the size and type of collimator used (eg., high resolution, low resolution, etc). As an extension, one may Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 22 Figure 2.8: An example of the geometrical contribution of an 8x8 voxel object to a single camera bin. For each camera bin, the resultant matrix gets reordered into a single row of the system matrix C. (Values shown are approximate.) consider the fully 3D case as well. In such a model, it can be seen that the contents of each camera bin actually originate from a cone of object voxels (Figure 2.9b)). When these improvements are made in the modelling of the system matrix, they are referred to as collimator blurring correction or resolution recovery because of the fact that when used in the reconstruction process, they result in improved spatial resolution in the reconstructed image. Following the creation of the system matrix, it then falls upon an iterative reconstruc-tion algorithm to convert the camera projection measurements into a series of transaxial images representing the three dimensional distribution of radiotracer within the body. In order to do this, one of several algorithms can be used. P E N A L I Z E D W E I G H T E D L E A S T S Q U A R E S ( P W L S ) One iterative reconstruction technique that is used often in research environments is the penalized least squares algorithm [24]. This technique is normally based upon the Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 23 a) b) ,BIU 1 B H i m i i 15811 H • • • • I : m J B IBB' i f lB Figure 2.9: System matrices modelling the case of a) collimator blurring in the transaxial direction only (2D), and b) true camera response with blurring in all directions (3D). Gaussian likelihood function (see Appendix A). Here an attempt is made to minimize the difference between the measured projection data, d and an estimated projection of the source distribution within the patient, C x . That is, minimize | ] T (wjk£CijkXi - djk^J + aR(x.) (2.2) where Wjk is a data weighting term and a i?(x) is an additional regularizing term normally used to control one aspect of the reconstruction process (eg, limit noise, enforce image smoothness, etc). Table 2.4 summarizes a few common regularizers used in penalized weighted least squares reconstructions. Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 24 Name Functional Form Tikhonov i?(x) = Ei x\ Boltzmann-Shannon #(x) = EiXilog(xi) Table 2.4: Examples of some of the regularizing functions used in least squares image reconstruction methods. M A X I M U M L I K E L I H O O D - E X P E C T A T I O N M A X I M I Z A T I O N ( M L - E M ) One of the most popular iterative reconstruction techniques for SPECT imaging is the maximum likelihood - expectation maximization (ML-EM) algorithm [25] [26]. This al-gorithm is based on the fact that the number of photons emitted in a certain time interval during radioactive decay can be described by a Poisson distribution. The derivation of the EM algorithm is given is Appendix A and so the theory will not be discussed any further here. The implementation of the E M algorithm includes the following steps, 1) Model the camera acquisition process and create a forward projection of the current guess, xfd of the object. Djk = Y,Ci>jkxf (2.3) v 2) Find the ratio between the current estimated projection data, Djk, and the true collected projection data in each camera bin, djk. Rik = jj^ (2-4) 3) Backproject a normalized ratio along the projection rays. = y ' p J * - (2.5) Chapter 2. Introduction to Nuclear Medicine and Image Reconstruction 25 4) Scale the current estimated object activity value, xfd for each pixel by this new coefficient, a;. xnew = xolda, ( 2 g ) The EM algorithm has been shown to always maximize the Poisson likelihood function and while always maximizing the likelihood function, after several iterations, the EM algorithm will attempt to fit the image to the noise present in the projection data thereby deteriorating the reconstructed image quality [27]. Thus, one must be cautious when using the EM algorithm in order to stop the reconstruction prior to substantial image degradation. Typically the stopping point is chosen based upon the experience of the user. 2 . 4 C H A P T E R S U M M A R Y In this chapter, an overview of the field of nuclear medicine imaging has been presented. Included in this overview has been a discussion of nuclear medicine equipment and com-monly used radionuclides, as well as a brief introduction to various reconstruction meth-ods that produce a three dimensional image from a series of two dimensional projection measurements. Each of the reconstruction methods discussed relies on the fact that the tracer distribution remains fixed within the object over the imaging time of the study. This is required so that each projection is consistent from one view to the next (ie., each projection angle views the same activity distribution). In the next chapter, we will see what occurs when the tracer distribution does not remain fixed in the body over the imaging time of the study, but rather is allowed to change. C H A P T E R 3 D Y N A M I C S P E C T A P P L I C A T I O N S A N D C U R R E N T M E T H O D S The reconstruction techniques discussed previously in Chapter 2 all rely on the radio-tracer distribution within the patient remaining static during the data acquisition process. However, the human body is a dynamic organism and changes in radiotracer distribution will undoubtedly occur due to elimination of waste products, metabolism, etc. Further-more, the rate at which a radiotracer distribution changes may reflect the functional ability of an organ (eg. slow elimination may reflect abnormal function) [28]. For this reason, it may be useful clinically to determine these physiological rates for diagnostic purposes. In this chapter, an introduction to dynamic imaging will be presented with clinical applications and a review of currently used methods of dynamic imaging in nu-clear medicine. 3.1 A P P L I C A T I O N S O F D Y N A M I C F U N C T I O N A L I M A G I N G Dynamic nuclear medicine imaging can be used in a wide variety of clinical circumstances. A summary of some possible applications is presented in Table 3.1, and described in detail 26 Chapter 3. Dynamic SPECT Applications and Current Methods 27 below. Application Radiopharmaceutical Purpose Myocardial Imaging 99m Tc-teboroxime Myocardial blood flow 1 2 31 labelled fatty acids Fatty acid metabolism 2 U iTl-Thallous Chloride Myocardial blood flow Renal Imaging 9 9 m T c - D T P A Glomerular filtration rate 9 9 m Tc-MAG3 Renal plasma flow GI tract imaging 9 9 m T c compounds Venous and portal blood flow Table 3.1: Examples of possible dynamic processes in the body and current radiophar-maceuticals that may be used in dynamic imaging. 3.1.1 M Y O C A R D I A L M E T A B O L I S M O F F A T T Y A C I D S Through a metabolic process known as catabolism, the myocardium uses fatty acids and oxygen to create the energy necessary to fuel the heart. It has beem suggested that the rate of fatty acid metabolism in the myocardium is a good indication of myocardial func-tion [29] [30] [31]. In regions of poor cardiovascular blood flow, the regional concentration of oxygen is low and so the metabolism of fatty acids is low since not enough oxygen is present to meet the metabolic demands necessary for the breakdown of these com-pounds. Similarly, in regions with enhanced myocardial blood flow, the rate of fatty acid metabolism is high as there is sufficient oxygen present to take part in the breakdown of fatty acids into the metabolic by-products. The fatty acids metabolised by the heart are part of a group referred to as essential fatty acids and include such compounds as linoleic and linolenic acid. A similar compound has been labelled with radioactive 1 2 3 I to form the radiotracer iodo-phenyl-pentadecanoic acid, or IPPA [32]. When IPPA is introduced into the body, some of it is taken up by the myocardium. Chapter 3. Dynamic SPECT Applications and Current Methods 28 As the amount of radiotracer injected into the body is small, after a short time interval most of the compound will have either been taken up by the heart muscle or located at some other point in the body (eg, liver, bladder, etc.). The IPPA that does get taken up by the myocardium will then be metabolised and broken down into metabolic by-products to be subsequently washed out. The rate at which these by-products are eliminated from the heart is a direct reflection of the myocardial rate of metabolism and thus of myocardial functional ability [33]. Therefore, by measuring the rate of washout from the heart of these by-products, one is able to infer the rate of fatty acid metabolism and hence, determine the relative functional ability of the myocardium [29]. Studies have been carried out utilizing IPPA with the dynamic planar imaging tech-niques discussed later in this chapter. It has been found that the washout of the by-products of fatty metabolism can be modelled by a dual exponential function in time [34]. The first component of this washout function is interpreted as the fast washout of catabolite from the heart, while the second, slower component represents the turnover of the radioisotope into various complex lipids. Under ischemic conditions in the heart, these kinetic rates are lengthened demonstrating decreased fatty acid turnover rates as a result of decreased oxygen levels [35] [36]. 3.1.2 M Y O C A R D I A L P E R F U S I O N Studies using 9 9 m T c labelled hexakis-2-methoxyisobutyl isonitrile (MIBI) are routinely performed in order to assess myocardial blood perfusion. Upon injection of MIBI into the patient, some of the pharmaceutical is taken up by the heart and gets locked into the myocardial cells to remain there over the course of the imaging time (as discussed in Chapter 2). The amount of pharmaceutical taken up regionally in the myocardium can be related to the amount of blood flow to that region. In this way, by using a Chapter 3. Dynamic SPECT Applications and Current Methods 29 three dimensional imaging technique, regions with poor blood flow (ischemic) can be differentiated from regions with normal blood perfusion. However, there are problems with using MIBI for myocardial viability studies. One of the major difficulties is in the fact that an increase in myocardial blood flow does not result in a proportional increase of tracer uptake [37] [38]. The result of this is that there may be relatively poor contrast between regions with marginally different perfusion. Thus, it may prove difficult for viewers to diagnose small, but potentially relevant decreases in myocardial perfusion. Because of this problem with conventional MIBI perfusion imaging, it may be more appropriate to assess myocardial function based upon actual blood flow measurements rather than on pharmaceutical uptake. One pharmaceutical that may be used for this is 9 9 m T c labelled teboroxime [39], [40]. It has been shown that the amount of teboroxime uptake.in the myocardium can be directly related to myocardial blood flow, that is, an increase in tissue perfusion is reflected by a proportional increase in tracer uptake [41]. Teboroxime however remains locked in the myocardium for much less time than MIBI and in fact has a washout half-life on the same order of magnitude as the data acquisition time [42] [38]. As a result, there are two competing processes occuring; that of increased contrast between normally perfused and hypoperfused regions due to the direct blood flow rate dependence on uptake, and that of decreased contrast from the fast washout rate. Since the rate of teboroxime uptake and the rate of subsequent washout are both directly related to blood flow, it would therefore be apparent that if a dynamic imaging technique were available that could determine these rates, then perhaps a more accurate assessment of myocardial perfusion could be obtained. Chapter 3. Dynamic SPECT Applications and Current Methods 30 3.1.3 R E N A L F U N C T I O N One of the most obvious choices for dynamic imaging is the case of renal imaging. The kidneys are constantly filtering blood and so represent a dynamic process with the rate of filtration being a direct measure of the kidney's functional ability. Every minute, approximately 625 ml of plasma pass through the kidneys. Of this plasma volume, approximately 20% (125 ml), gets filtered through the glomerulus of the kidney. This filtrate consists mainly of water, but other components are also present in smaller quantities such as sodium, calcium, urea and other trace substances. Most of this glomerular filtrate gets subsequently reabsorbed by the renal tubules as it passes through the kidney, while those compounds not reabsorbed get passed into the bladder where they are excreted in the urine. Because much of the waste products excreted from the body is first filtered at the glomerulus, the functional state of the kidney is dependent upon the glomeruler filtration rate. In a normal, healthy person, the glomerular filtration rate ( G F R ) is approximately 125 ml/min. One method to determine this rate is through the injection of a substance into the bloodstream that gets completely filtered by the glomerulus. By measuring the venous blood concentration of the substance at various time intervals, one is able to determine the amount of compound removed from the blood per unit time [19]. Any measurement made in this manner however, will provide a quantitative value for the renal system as a whole rather than regionally for each kidney. In such a test if it is discovered that the G F R is not within normal values, follow-up testing must be performed in order to determine which kidney is malfunctioning. Because of this, a three dimensional imaging technique that is able to determine values for the glomerular filtration rate regionally within the kidneys would provide more useful diagnostic information as to what intervention steps are required in order to provide the most benefit to the patient. Chapter 3. Dynamic SPECT Applications and Current Methods 31 3 .1 .4 P H A R M A C O K I N E T I C ( P K ) M O D E L L I N G Prior to the clinical acceptance of new pharmacological agents, the kinetic behaviour of the drug within the body must first be demonstrated. That is, it must be known where the drug localizes in the body, as well as how long it remains there. This typically in-volves administration of the drug, either by intravenous or oral means, and subsequent measurements of plasma and/or urine drug concentrations at various times. These drug concentration measurements are made so that a (concentration vs. time) curve can be constructed in order to assess the global pharmacokinetic rates of distribution and elimi-nation of the drug. In the case of nuclear medicine imaging, it is more common to refer to (concentration vs. time) curves as (time activity) or TA curves, as the amount of activity per unit volume is plotted vs time. Depending on the interactions of the drug within the body, such curves are usually used with a technique known as compartmental modelling. A compartmental model is a simplified view of how a drug is taken up and eliminated from the body. The simplest compartmental model is that of a single compartment (Fig. 3.1). In this model, the assumption is made that all tissues in the body are in equilibrium with the blood at all times, that is, the drug concentration in the blood accurately reflects the equivalent drug concentration in the tissues. Therefore, by measuring the venous blood concentration, it is assumed that one is also measuring the same drug concentration in the tissues as well. Over time, the tissues will eliminate some of the drug through excretory means such as the renal system. The rate at which the drug is excreted can be determined by fitting the appropriate model to the corresponding time activity curve. In the aforementioned single compartment model, the appropriate function to fit the TA curve to is a monoexpo-nential function. Once the fit parameters are determined, other relevant pharmacokinetic parameters, such as the volume of distribution, V^, can also be determined. Chapter 3. Dynamic SPECT Applications and Current Methods 32 A(t) M I = - k A ( t ) dt A(t) = A e - k t k Figure 3.1: A simple one compartment model used in PK modelling. Following injection of the drug into the compartment, the amount of drug within the compartment (A(t)), is reduced in an exponentially decreasing manner with decay rate k. While a simple single compartment model may accurately describe the kinetics of some drugs, the majority of drugs do not obey such a model. For most drugs, it is necessary to subdivide the tissue space into more kinetic compartments that behave differently. These compartments typically would have very different drug concentrations, and hence, altered pharmacokinetic rates, yet all of these pharmacokinetic parameters have to be found from a single TA curve. The main problem with the determination of PK parameters from plasma drug con-centrations is as models become more complicated, it becomes increasingly difficult to fit the measured data in order to solve for the large number of kinetic parameters required. Additionally, it is very difficult to determine the drug concentrations at particular sites in the body. Rather, it must be assumed in the model that all regions of the same organ act homogeneously with respect to the administered drug. For these reasons, it is hoped that a time dependent 3D imaging technique will be able to provide accurate pharmacokinetic information on a regional basis in the body. With such a technique, it becomes possible to determine TA curves for each small volume in the body therefore leading to more accurate measurements of the pharmacokinetic parameters present. Chapter 3. Dynamic SPECT Applications and Current Methods 33 3.1.5 R E D U C T I O N O F D Y N A M I C A R T I F A C T S I N S P E C T I M A G I N G In many circumstances, certain regions of the body may undergo dynamic behaviour of a radiopharmaceutical while other regions of the body have a static activity distribution. Consider, for example, the case of pelvic imaging with a bone imaging agent. Under ideal circumstances, only the bones will take up the injected radiopharmaceutical, but realistically, only a small percentage of the injected dose is actually extracted by the bones. Much of the injected dose instead gets filtered in the kidneys and ends up stored in the bladder as urine. As the kidneys continually filter the blood, the amount of radio-pharmaceutical in the bladder will continue to increase over time. However, the amount of activity present in the bones remains essentially static over the imaging time and so projection data combines the bladder, which is constantly increasing in activity, with bones which contain a static activity distribution. As discussed below, projection data of this type leads to image artifacts when reconstructed with conventional reconstruction routines. As seen in Figure 3.2, when a changing activity distribution is reconstructed using conventional image reconstruction methods, severe streaking artifacts are seen in the resultant images. These artifacts appear because of the fact that the radiotracer dis-tribution within the patient changes from one projection to the next (referred to as inconsistent projections) and the reconstruction routine in unable to account for these changes [43]. In fact, it has been found that in order to avoid severe artifacts when using the filtered backprojection method, the time of rotation of the SPECT camera should be at least twice as fast as the rate of activity change in the object (eg., to reconstruct a 2 min washout half-life, a full SPECT camera rotation should be performed in 1 min) [44]. If however, a reconstruction procedure were available that could account for these temporal inconsistencies, then it would become possible to reconstruct a changing activity Chapter 3. Dynamic SPECT Applications and Current Methods 34 distribution from inconsistent projection data. Once the dynamic regions are accounted for in a reconstruction, dynamically induced artifacts can be reduced since it is these dynamic regions that give rise to the inconsistent projection data in the first place. Figure 3.2: Shown on the right is a reconstructed image using the standard filtered backprojection routine of an object undergoing a changing activity distribution. On the left is the true object. Notice the lack of definition in the lung regions and the streaks emanating from the leftmost heart region in the reconstructed image. 3 . 2 D Y N A M I C F U N C T I O N A L I M A G I N G M E T H O D S For the applications mentioned above, dynamic imaging offers the potential for a profound change in the way nuclear medicine diagnostic procedures are performed. To this end, several techniques currently exist in nuclear medicine to obtain kinetic parameters from these studies. 3.2.1 P L A N A R D Y N A M I C I M A G I N G The most well established dynamic imaging technique is known as planar dynamic imag-ing. For this technique, the gamma camera is positioned over the region of interest to be imaged (eg., kidney, heart, etc) and remains in this position during the entire scan Chapter 3. Dynamic SPECT Applications and Current Methods 35 time. A pharmaceutical is injected into the patient and projection images are acquired at this angle in fixed time intervals, for example, every 5 seconds over a total time of 20 minutes. The resultant images then depict "snapshots" of the tracer distribution at the time each projection was acquired (see Figure 3.3). In order to analyze such a dynamic sequence of images and to obtain the desired kinetic information, further analysis is carried out. To do this, a region of interest is selected either by user interaction, or via a computerized region drawing program. This region of interest is then applied to each frame in the sequence and the total number of counts determined within each region at each time interval. This data is then plotted against time as a time activity curve (TAC). Once the TAC curve is plotted it is then possible to fit an appropriate time dependent function (eg., mono exponential decay) to this data in order to obtain estimates of kinetic parameters. The time dependent function to fit to the data depends on the compartmental model chosen to represent the dynamics of the study. 4500 4000 3500 3000 2500 , r 2 0 0 0 1500 I 1000 0 5 10 15 20 25 Figure 3.3: Examples of a dynamic planar image sequence of a renogram acquired at 10 second intervals. Images shown correspond to the times 1 min, 5 min, 10 min and 15 min after radiotracer injection. On the right is a plot depicting the amount of activity present in each kidney as a function of time (a time activity curve) Chapter 3. Dynamic SPECT Applications and Current Methods 36 This method of dynamic imaging has the drawback that images are two dimensional representations of a three dimensional object and as such, lack any information regarding the three dimensional spatial distribution of radiotracer. This is not a large concern when imaging well separated organs such as the kidneys, but becomes increasingly problematic for organs such as the heart or brain where the determination of local three dimensional variations in the kinetic parameters is desirable. 3 . 2 . 2 P O S I T R O N E M I S S I O N T O M O G R A P H Y In comparison to the previously described dynamic imaging technique which uses single photon emitting isotopes in order to produce two dimensional dynamic information, positron emission tomography (PET) can be used to gather three dimensional dynamic information. The way in which this data is gathered is different from planar nuclear medicine and utilizes special equipment. While not directly comparable to the case of SPECT imaging, dynamic PET will be discussed in order to present the importance of dynamic information to clinical situations and to provide a gold standard in order to compare the methods of dynamic SPECT imaging. A PET camera consists of a circular array of detectors surrounding the patient. Usu-ally these detectors consist of several adjacent rings of a high density scintillating material (eg., bismuth germanate, BGO). The entire detector array may be about 40 to 60 cm in diameter and 10 cm long for the case of brain imaging, or up 100 cm in diameter and 15 cm long for whole body imaging. Surrounding the detector ring are a series of pho-tomultiplier tubes in order to amplify the small amount of light produced in the photon detection process. When a positron emitting isotope decays, the emitted positron travels a short distance before annihilating with an electron and producing two 511 keV photons. These photons Chapter 3. Dynamic SPECT Applications and Current Methods 37 are emitted from the annihilation centre in opposing directons to each other. The location of this annihilation event occurs in the vicinity of the radiopharmaceutical, and thus, if the location of the annihilation is known, then it is possible to infer the location of the radiopharmaceutical. In order to determine the location of the annihilation event, one must detect both emitted photons. When an incoming photon is detected in a detector bin, the opposing photon must also be detected by another detector within a specific short time interval (a few ns). If such a dual event occurs (called a coincidence event), it is known that the positron annihilation happened somewhere along the line connecting these two detectors. Since each detector bin is in coincidence with several other detector bins, it is possible to collect photons from all angles simultaneously, and hence there is no requirement for a camera rotation as in the case of SPECT imaging in order to collect angular projection data. Since PET cameras don't require a collimator to limit the angles of the incoming photons, and since most isotopes used in PET have a very short half-life, it is possible to collect a large number of statistical events within a very short time interval. As a result, enough data can be obtained in just a few seconds in order to reconstruct a reasonable transaxial image. Furthermore, many sets of projection data can be obtained sequentially and reconstructed into a series of transaxial slices depicting the radiotracer distribution at a variety of times. Analysis of these sequential images is carried out in a manner similar to that of dynamic planar imaging with the exception that usually 3D regions of interest are used instead of only 2D regions. Figure 3.4 depicts some typical dynamic compartmental models used to fit reconstructed 3D TAC's from PET data. Dynamic PET imaging is a very useful 3D dynamic imaging technique that is able to provide an abundance of data related to many physiological conditions, however, the Chapter 3. Dynamic SPECT Applications and Current Methods 38 FDG Model Blood Space FDG FDG FDG 6-phosphate Tissue Space Ammonia Model Blood Space Fast pool Slow pool Tissue Space Carbon Dioxide Model Blood Space F-Dopa Model Blood Space F-Dopa F-Dopa F-Dopa & Metabolites Tissue Space Figure 3.4: Examples of kinetic models used in dynamic positron emission tomography (PET) imaging. additional requirements needed for PET (eg., high cost, radionuclide availability, etc.) often put it out of the reach of most clinical departments. 3.3 C U R R E N T D Y N A M I C S P E C T M E T H O D S Because of the relatively high cost associated with positron emission tomography and with the limited information garnered in planar imaging, there is a need for an imaging technique that combines the advantages of both these methods, that is, one that is commonly available like planar imaging yet is able to provide three dimensional dynamic information similar to PET. One suggestion for this imaging method is to use existing nuclear medicine gamma cameras in order to perform dynamic single photon emission Chapter 3. Dynamic SPECT Applications and Current Methods 39 computed tomography (SPECT). As mentioned previously in Chapter 2, three dimensional information of a static tracer distribution can be obtained through the rotation of the gamma camera around the object of interest and subsequent data processing. In this form of SPECT imaging, no dynamic information is available as the reconstructed images represent the static distribution of tracer over the imaging time of the experiment. In order to obtain dynamic information from such a study, several techniques have been proposed involving altered acquisition procedures or advanced data handling and reconstruction methods. 3.3.1 F A S T R O T A T I O N I M A G E B A S E D S P E C T If a gamma camera rotates fast enough around an object, then it can be assumed that the tracer distribution within the object will change very little over the rotation time. The projection measurements can then be considered as temporally consistent and a con-ventional reconstruction can be performed with minimal streaking artifact. If successive rotations are performed and reconstructed in a similar manner the result is a series of three dimensional images of the same anatomical location, but with each image depict-ing a slightly different radiotracer distribution as a result of the time shift in the data collection sequence. Time activity curve analysis of this four dimensional data set then results in information related to the dynamics of the injected tracer. The critical assumption made in this type of dynamic imaging method is that the tracer distribution within the body remains essentially static during each camera rota-tion. In order to reconstruct an image sequence with high temporal resolution, one must rotate the camera very quickly around the patient with minimal delay between successive rotations. However, by rotating the camera quickly around the patient, the number of detected events in the camera is small, thus resulting in poor projection statistics and a Chapter 3. Dynamic SPECT Applications and Current Methods 40 very low signal to noise ratio. This low signal to noise ratio produces poor reconstructed images and results in a decreased spatial resolution due to the requirement for more heavily smoothing filters in the reconstruction. Slowing the camera acquisition proce-dure increases the spatial resolution, but at the expense of a lower temporal resolution and with possible streaking artifacts present. In order to increase the number of detected events and therefore improve the spatial resolution, multiple detector cameras are usually used for this type of dynamic imaging. As a triple head camera collects three times the number of counts as a single head camera, the resultant increase in the signal to noise ratio will be a factor of In a typical SPECT data acquisition procedure, the gamma camera usually rotates to a specified angle, stops at that position to acquire data for a few seconds, and then moves on to the next angle where the process is repeated. This is called a step-and-shoot acquisition because the camera only acquires projection data while it is stopped at each projection angle. In a conventional SPECT acquisition the time needed for the rotation motion is not usually a problem as the time spent at each projection is much longer than the time needed for the rotation (« 1 sec for camera movement). However, if the camera is rotated quickly, this movement time sets a limit on the speed of acquisition. Indeed, it is entirely possible that when using a step-and-shoot acquisiton protocol with a fast rotation, the camera spends more time moving and not collecting data than actually making measurements. In answer to this data collection problem, most camera manufacturers have allowed for a continuous collection-rotation protocol. In this protocol, the camera rotates contin-uously around the patient, collecting projection measurements as it proceeds. All data collected oyer a small rotation angle, typically, 3 ° or 6 ° , is then summed together into a projection measurement assumed to be taken at the midpoint of the rotation. With this Chapter 3. Dynamic SPECT Applications and Current Methods 41 acquisition method, there is no point in the rotation where the camera is not collecting data. To collect data that is the most temporally consistent from one rotation to the next, it would be desirable to minimize the amount of time required for the camera to complete one rotation and to start the next. The amount of data collected from a single 120 degree rotation of a triple head SPECT camera is of the order of 16 MB. In some cameras, this data is transferred from the acquisition computer to the analysis computer at the end of the current rotation, a procedure that may take several seconds. As well, most cameras are not capable of continuously rotating in the same direction but instead must reverse direction at the end of each rotation. This reversal may take a few seconds in order to occur and must be considered when performing a data acquisition procedure of this type. Despite the concerns with this type of dynamic SPECT imaging, the fast rotation method is the most common method of dynamic SPECT imaging in use today [45] [46] [47] [48] [35]. W E I G H T E D I N T E R P O L A T I O N Another reconstruction method that has been used with SPECT imaging in order to obtain dynamic information of a three dimensional tracer distribution is the weighted interpolation reconstruction method. Originally this method was used in order to reduce the dynamic artifact present in SPECT imaging of the gallbladder [49], but has since been extended to provide dynamic information as well [50]. This method involves a change in the conventional SPECT acquisition protocol sim-ilar to that of the multiple fast rotations, that is, rather than a single slow rotation of a SPECT camera, the acquisition is performed as several faster sequential rotations. Chapter 3. Dynamic SPECT Applications and Current Methods 4 2 Once this projection data is acquired, it becomes possible to create projection data corre-sponding to a selected angle at an arbitrary time by means of an interpolation of the data collected at the same projection angle but at different times. Consider two sequential projections, A and B, acquired at the same angle and at times ta and tb. Another pro-jection at the same angle, C, at time tc (ta < tc < tb), can be generated from projection data A and B by the interpolation given by, c = h ^ t 1 . A + t c ^ t a B ( 3 1 ) tb ta tb ta If multiple rotations are acquired, then it becomes possible to create interpolated projection data at any time in the data acquistion sequence. Since many sets of projection data can be generated, then it is possible to reconstruct this interpolated projection data sequence with conventional reconstruction methods in order to obtain images at each interpolated time frame. While the method appears to work well with slow activity changes, it suffers from the fact that each reconstructed image is produced from interpolated data and does not necessarily represent a true measured quantity. For this reason, only limited work has been done with this method. 3.3.2 S L O W R O T A T I O N D Y N A M I C S P E C T M E T H O D S Because of the complexity of the data acquisition process and poor image quality resulting from the fast rotation dynamic SPECT, other methods are being investigated which are able to provide three dimensional dynamic information -while utilizing an acquisition protocol that maintains the compromise between high statistical data and high temporal resolution. These methods usually involve only a few slow camera rotations in order to collect high count projections yet are able to account for the changing radiotracer Chapter 3. Dynamic SPECT Applications and Current Methods 43 distribution from one projection to the next. D I R E C T P A R A M E T E R E S T I M A T I O N M E T H O D S In the cases of dynamic PET and fast rotation dynamic SPECT, a sequence of images is reconstructed corresponding to a sequence of time frames. With this information, dynamic analysis is performed in order to determine the kinetic parameters of the ra-diotracer distribution within the object according to the selected compartmental model. One method of dynamic SPECT reconstruction however, bypasses the image reconstruc-tion step to provide estimates of dynamic parameters directly from the projection data. This method is referred to as the direct parameter estimation technique [51] [52] [53] [54]. With this method, one must first have some prior dynamic model that describes the expected dynamic behaviour for a selected region. Once this time dependent function is established, the measured data in each camera bin and at each camera angle can be described in terms of the projection of the time dependent object, x(t) at each projection time. Thus in a manner similar to the projection described in Chapter 2, the dynamic projection data collected can be modelled as, Cx(t) = d(t) • (3.2) where C represents the system matrix describing the projection operation, x(t) is the time dependent object and d(t) is the measured projection data at a variety of times. For example, in the case of fatty acid metabolism discussed previously, one might assume that the activity in each voxel can be described by a dual exponential function in time. Thus, the activity in each voxel, i, behaves according to the function, Xi{t) = one-Xit + pie~Uit + 7i (3.3) Chapter 3. Dynamic SPECT Applications and Current Methods 44 where C K ; , fa, A; and Vi correspond to the dynamic parameters, and 7, represents a static component present in each voxel. It can be noticed that in the absence of dynamic behaviour (ie., cxi = fa = 0), equation (3.2) reduces to the conventional case of static activity distribution. The problem in this image reconstruction procedure is to determine the functional parameters of each object voxel based upon the collected data. In order to determine these parameters, typically a nonlinear optimization process is performed whereby the difference between a current estimated object activity and the true collected data is minimized. In the case of fatty acid metabolism discussed, the objective function to minimize is: With this method of dynamic SPECT image reconstruction, rather than the recon-structed images representing the activity distribution within the object at specified times, the dynamic parameters a;, fa, A;, and 7$ are obtained instead. The inherent shortfall in this type of image reconstruction however, is in the fact that the objective function is nonlinear and is comprised of a large sum of multiple exponential functions. Such a sum of exponentials is extremely difficult to fit mathematically and does not necessarily result in a unique solution. 3.3 .3 L I N E A R M E T H O D S Because of the difficulty in fitting the nonlinear objective function from equation (3.4), other methods have been developed which attempt to linearize the objective function, yet still provide the a priori modelling information present in the direct parameter approaches (3.4) [55] [56]. Chapter 3. Dynamic SPECT Applications and Current Methods 45 F A C T O R A N A L Y S I S O F D Y N A M I C S T R U C T U R E S (FADS) In the factor analysis of dynamic structures (FADS) method, it is assumed that the dynamic behaviour of activity within each pixel can be described in terms of a linear combination of known discretized time dependent functions. For example, the time dependent functions may be taken to be exponential functions with known discretized kinetic parameters over a finite range. Thus, at each kth time interval, the amount of activity in the ith voxel can be described as, N Xi(tk) = AaF^tk) + Ai2F2(tk) + Ai3F3(tk) + ... = £ AinFn{tk) (3.5) n=l where Ain is the coefficient of the nth dynamic factor, Fn. In the case where the factors are given by a range of monoexponential functions, then the activity present within the zth voxel at time tk will be given as, N Xi(tk) = Aine-Xntk (3.6) n=l Using the projection operator, C, as presented previously, we can produce a projection of this time dependent object. Thus, the projection of the object at the kth projection time will simply be equal to, E ( E E AinFn(tk)) = djk (3.7) j,k \ i n ) Prior to the reconstruction, the functional form of the dynamic factors Fn(t) must already be established. Thus it falls upon the reconstruction process to determine the values of the coefficients, Ain for each object voxel. This method has been successfully applied to dynamic PET imaging in order to extract kinetic parameter information from regions of interest given a reconstructed dy-namic three dimensional dataset [57] as well as from the data acquisition space where one Chapter 3. Dynamic SPECT Applications and Current Methods 46 determines the kinetic parameters directly from projection measurements [58]. However, attempts to extend the FADS method into dynamic SPECT imaging have been less suc-cessful [59]. While the FADS method converts the nonlinear problem of equation (3.4) into a linear problem, the difficulty of nonuniqueness still exists. It has been shown that using a priori knowledge aids in this difficulty, but in clinical images, such information does not usually exist [60]. An additional concern in the FADS approach is that one must be selective in the appropriate dynamic factors used in the image reconstruction process. That is, enough factors must be selected in order to accurately describe the dynamic behaviour present, but not too many factors such that the computation time is exceedingly long [61]. As an alternative to specifying factors explicitly, some work has been done in attempting to determine the factors from the dynamic projection data itself [59]. The main drawback with this is in the computation time required as the reconstruction process must not only determine the factor coefficients, the Aj„'s, but also must determine the most likely factors, the F„'s. Therefore, in order to speed up the reconstruction process, typically a limited number of factors is chosen. Such a procedure limits the usefulness of FADS as a clinical dynamic imaging technique since it is often small variations in dynamic parameters that are informative from a diagnostic point and may not be distinguished with this method. 3 . 4 C H A P T E R S U M M A R Y In this chapter, potential applications of three dimensional dynamic imaging have been presented along with various methods that have been proposed to determine this dynamic behaviour. Most of these methods however either require special acquisition hardware (fast rotation SPECT), or make assumptions related to the expected kinetic behaviour Chapter 3. Dynamic SPECT Applications and Current Methods 4 7 (direct parameter estimation SPECT, factor analysis). Because of these limitations, dynamic SPECT imaging is not routinely used for clinical work, but rather is limited to research environments for the time being. In the following chapters, a new image reconstruction algorithm is presented which is able to determine the changing three dimensional radionuclide distribution within an object while using conventional SPECT imaging practices. C H A P T E R 4 T H E D S P E C T METHOD The various shortcomings of the dynamic imaging techniques discussed in Chapter 3 indicate the need for a new method of dynamic SPECT imaging. The method that is being proposed in this work uses a conventional SPECT camera and a typical SPECT acquisition protocol (ie, a single, slow camera rotation), and thus overcomes the problems associated with the difficult data acquisition process required with some of the other imaging methods. From this point on, the term dSPECT will be used to refer to this method of imaging where a single slow camera rotation is used. 4.1 D S P E C T R E Q U I R E M E N T S The dSPECT reconstruction method has been designed as a clinically versatile technique that is accessible to most hospitals. Because of this, the method requires no special cam-era equipment but rather is capable of using a wide variety of available SPECT cameras including single, dual and triple head units. From a data collection standpoint, again the dSPECT method was designed to work with a variety of acquisition techniques, in-cluding the most conventional method (a single 180° rotation). Because of the generality of the method and the fact that reconstruction is performed using iterative techniques, 48 Chapter 4. The dSPECT Method 49 there is no limitation on the complexity of data acquisition. For example, acquisitions performed using multiple camera rotations, gated acquisitions, different collimators can all be included in the dSPECT reconstruction method proposed. An additional feature of the dSPECT reconstruction method is in the correction for patient specific attenuation. This is a necessity since changes in tracer distribution due to dynamic processes cannot otherwise be distinguished from changes in activity measurements due to attenuation effects. This can be explained with the aid of Figure 4.1. a) b) Figure 4.1: Diagram depicting similar effects yielded in projection data due to a changing tracer distribution or because of attenuation effects with a static distribution. In a) a dynamic object is present in the absence of attenuation, while b) represents a static object with attenuation. In both cases, 180° projection measurements were acquired in a counterclockwise rotation with the camera starting at the 6 o'clock position. In the upper portion of Figure 4.1 are projection measurements of a dynamic object in the absence of attenuation. In the lower half is a static object acquired with the same acquisition rotation, but this time the object is placed off center in the medium such that Chapter 4. The dSPECT Method 50 the amount of attenuation increases at each projection stop. Notice that the resultant projection measurements appear similar in both cases. In order to obtain patient specific attenuation maps, many approaches have been used [62] [63] [64] [65]. Most of the methods are similar, however, in the fact that they tend to use a second radioactive source that is external to the patient to act as a transmission source. This additional radionuclide usually emits photons of a different energy than those arising from the emission agent in order that the detected photons from both the transmission and the emission sources can be differentiated on the basis of their energies. All experimental studies involving object attenuation presented in this work have used the Siemens Profile attenuation device [66]. This device consists of a series of metal rods filled with 1 5 3 G d . The sources in this apparatus are placed in a parallel line array directly opposite to the detector face (see Figure 4.2). The sources are further collimated with a parallel hole collimator in order to ensure that the angle of emission is small. In the absence of an attenuating medium, the number of photons emitted from the transmission source array and detected by the opposing detector is IQ. Once an object is placed between the transmission source array and the detector however, the number of counts detected by the camera will decrease in accordance with the amount of attenuating material. The new number of counts detected, / , is related to the original, non-attenuated counts, lo, by the equation, I = I0 exp (-Y^Pih) (4.1) i where pi is the linear attenuation coefficient for each ith object voxel, and U is the length of each object voxel that the emitted photons pass through. By making a series of transmission measurements around the patient, one can arrive at a series of equations of the form (4.1) for a variety of angles and camera bins which can then be reconstructed Chapter 4. The dSPECT Method 51 Figure 4.2: Schematic representation of the Siemens Profile transmission imaging geom-etry. Transmission imaging is performed simultaneously with emission imaging through the use of radioisotopes with different emission energies. into a three dimensional image of the patient specific attenuation coefficients, LI. These attenuation values are then used in the reconstruction process in order to account for the decrease in detected counts due to attenuation. Further discussion of attenuation correction in SPECT reconstruction can be found in references [67], [68], [69]. 4 . 2 T H E O R Y O F D S P E C T I M A G E R E C O N S T R U C T I O N Because in SPECT it is impossible to acquire all projection angles simultaneously as in PET, projection measurements must be acquired sequentially resulting in a time differ-ence between each measurement. When reconstructed using standard procedures, it is Chapter 4. The dSPECT Method 52 assumed that all projections were acquired at exactly the same time and that no change in radiotracer distribution has occured between projections. This is important as multi-ple views of the same object (ie., a static radiotracer distribution) are required in order to accurately reconstruct a three dimensional object. The dSPECT method is different in this regard as it relies on the fact that all pro-jections were not acquired at the same moment in time but rather that each projection is acquired at a slightly different time and that the tracer distribution does not need to remain stationary over the entire scanning interval. The end result of a dSPECT recon-struction is a series of three dimensional images each representing the tracer distribution at the midpoint time of each projection measurement. It is possible to achieve this by us-ing the projection data and the time information of each projection measurement, along with a set of linear inequality constraints that represent the dynamic behaviour of the tracer distribution within each object voxel. 4.2.1 L I N E A R I N E Q U A L I T Y C O N S T R A I N T S Initially, the dSPECT method was designed for the case of myocardial metabolic imaging using radiolabeled fatty acids. As mentioned in Chapter 3, the washout of a fatty acid tracer from the heart has been modelled by a dual exponential function in time. Thus, if such a study is performed, one would inject a radiolabeled fatty acid into the patient, wait until the tracer is taken up by the heart, and then start the imaging procedure at the point where the tracer starts to wash out of the heart. Using such an acquisition technique ensures that information is obtained relating to the decreasing activity within the heart pertaining to myocardial fatty acid metabolism. As an example of such a technique, consider the situation presented in Figure 4.3. In this example, the true activity contained within the myocardium can be represented Chapter 4. The dSPECT Method 53 by the time activity curve shown on the right. This curve consists of the ini t ial uptake of tracer (region A ) , the point where the activity is a maximum (region B ) , and the subsequent dual exponential washout (region C ) . A B C Time Figure 4.3: A simplified representation of an expected myocardial time activity curve for fatty acids. Consider the case of washout only from the myocardium, that is, region C in the T A curve in Figure 4.3. If projection measurements are acquired at equally spaced time intervals, then for each point in time, the amount of activity wi thin the myocardium wi l l always be less than at the previous time. If we now look at a single voxel of the myocardium, we see that again it follows the same decrease in activity over time. This can be rewritten such that for the ith voxel, the amount of activity present at projection time t\ is greater than the activity contained in the same voxel at the next projection time, t2, which again is greater than the activity at t3, and so forth up to the final projection measurement acquired at time ts. XiH > xit2 > ... > xit, (4.2) As the amount of activity contained within each voxel cannot be negative, we also Chapter 4. The dSPECT Method 54 must have a positivity condition on the activity. If we consider as well that not all voxels will contain activity, then any voxel with zero activity will remain zero over the entire scan time. Modifying equation (4.2) to reflect these changes we have, Xity > xit2 > ... > xiu > 0 (4.3) We now have formed a linear inequality condition that relates the amount of activity present in any given voxel at any time to the amount of activity present in the voxel at any other time for the case of a decreasing tracer distribution over time. If we now turn our attention to region A of the curve shown in Figure 4.3, we see that over this time interval the activity in each voxel increases from one time to the next. In a manner similar to that described above, we can formulate an inequality relating to this increase in voxel activity over time. In this case, the inequality becomes, 0 < xitl < xit2 < ... < xiu (4.4) In the final scenario, we will consider the case over the region B in Figure 4.3. This region contains both increasing and decreasing portions over the scan time, and so we can develop arevised inequality condition that depicts this. As can be seen, the activity in a voxel will increase up to the time tp, followed by a decrease after this point. Thus, by making use of both equations (4.3) and (4.4), we arrive at the inequality, 0 <xitl < ... < < xitp > xitp+1 > ... > xits > 0 (4.5) Equations (4.3) through (4.5), present inequalities which relate the amount of activity within the zth voxel at any given time, to the activity in the same voxel at a different time. Note that it is possible that each voxel in the object may have a different inequality condition such that some voxels may contain only an increasing activity, while others may Chapter 4. The dSPECT Method 55 contain only a decreasing activity, while still others may undergo combined increasing and decreasing activity distributions. With this concept, we now turn our attention to the problem of determining the unknown voxel activity at the time of each projection measurement. This problem can be posed as an extension of the general image reconstruction principle as discussed in Chapter 2. As mentioned previously, the image reconstruction can be regarded as a matrix-vector product, C x = d (4.6) where x is a vector representing the amount of activity present in each object voxel, d is the vector containing the measured projection data, and the matrix C is the system matrix which maps the object space into the projection data space. In the problem of tomographic -imaging, the question is to determine the unknown values of x such that they best match the measured data d when multiplied by the matrix C. When the situation involving a changing radionuclide distribution is included, the problem remains similar except for the fact that the voxel activity changes at each projection stop. We have seen that a variety of different dynamic behaviours can be represented by the inequalities given in (4.3) - (4.5). By incorporating these inequalities into the recon-struction process, it is possible to reconstruct the activity distribution within the object over time. That is, the reconstruction process attempts to determine the unknown ac-tivity within each voxel at each moment in time by solving equation (4.6) subject to the dynamic activity constraints given by (4.3) - (4.5). Previous image reconstruction routines do not allow the inclusion of such complex Chapter 4. The dSPECT Method 56 constraints such as those representing the dynamic behaviour of the radiotracer. There-fore, the implementation of these constraints in the image reconstruction process is the subject of the following section. 4.3 C O M P O N E N T S O F T H E D S P E C T M E T H O D In order to perform a dSPECT reconstruction, several items must first be defined. These consist of: i) the voxel activity vector, ii) the probability system matrix, iii) the projection data vector, and iv) the linear inequality difference matrix. These will now be presented in detail. 4.3.1 T H E A C T I V I T Y V E C T O R As mentioned in Chapter 2, each transaxial slice of an object is made up of many small volume elements (voxels). In the case of static SPECT, for each transaxial slice, voxel indices are numbered sequentially and thus each two dimensional slice can be represented as a vector of unknown activity values. In the case of dynamic imaging, each time frame corresponds to a different radiotracer distribution within the object. However, there is a link between successive images in time because of the inequality constraints already discussed. In a dSPECT reconstruction, each voxel will consist of a total of (# of proj) different activity values (one for each projection time). If voxels are numbered along the time axis first followed by spatially, then the unknown activity in each voxel at each time can be ordered as a single unknown activity vector of length (# of object voxels)^ of proj) (see Figure 4.4). Chapter 4. The dSPECT Method 57 i 2 3 16 17 I 12 4081 4082 4096 Figure 4.4: The reordering of object voxels from a single transaxial slice acquired at 16 time frames into a spatio-temporal vector for use in the dSPECT image reconstruction process. Notice voxel indices proceed over time first, followed by spatially. 4.3.2 T H E P R O B A B I L I T Y S Y S T E M M A T R I X Similar to the case of static SPECT reconstructions, in dSPECT, there is a probability system matrix that maps the activity contained within each voxel into the appropriate camera bin. In the case of a static SPECT reconstruction, the mapping operation can be represented by the matrix operator, C. The dimensionality of the static probability matrix is of size, [(# of bins) (# of proj) (# of heads) x (# of object voxels)} If all elements are stored for this matrix, a typical size for a single head camera might be [(64) (64) (1) x (4096)] w 17 million elements, however we are aided by the fact that the matrix is quite sparse (?» 2 million non-zero elements) and thus storage requirements are lower than first anticipated. In the case of dSPECT, at each projection angle a different transaxial object is present and so the mapping operation must take this into account. Thus, at each angle the Chapter 4. The dSPECT Method 58 dSPECT system matrix maps the current radionuclide distribution into the appropriate camera bins. Since the unknown activity vector for the dSPECT problem consists of [(# of object voxels) (# of proj)] variables, the new dSPECT system matrix must also be of this dimension. Thus the size of this matrix is then, [(# of bins) (# of proj) (# of heads) x (# of object voxels) (# of proj)} At first this size seems impossible to work with as a typical dynamic transaxial slice system matrix acquired with a single head camera will be of size [(64) (64) (1) x (4096) (64) ] « 1 billion elements. However, the camera acquisition geometry has not changed but rather only the object activity vector is different at each projection stop. As a result, the number of non-zero elements in the dSPECT system matrix has not increased compared to the static system matrix. The only difference is that the dynamic system matrix has been stretched out compared to the static matrix (see 4.7). For each ith voxel in the object, the components are no longer distributed along just a single column, but rather now are stretched out into the number of projections that were acquired. That is, Ci a n 12 ijk C(i, bins, proj) I d = a ill Ci ill C, ijk \ (4.7) y Cri,bins,proj) These individual system matrices are then further combined into a large system tensor that is comprised of all the object voxel elements, Chapter 4. The dSPECT Method 59 Ci C2 ... Cvoxeis (4.8) where C has a dimensionality of, [(# of proj) (# of bins) (# of heads) x (# 0 / voxels) (# 0 / proj)] 4.3.3 T H E D A T A V E C T O R Because the data collection process for dSPECT is the same as for a static SPECT reconstruction, it should be apparent that the measured projections consist of the same information for dSPECT as for static SPECT. That is, for each projection measurement, the data first gets divided into transaxial slice data comprised of transaxial sinograms and then further reordered into a data vector of length (# of bins) (# of proj). 4.3.4 T H E L I N E A R I N E Q U A L I T Y D I F F E R E N C E T E N S O R As discussed earlier, the dSPECT method reconstructs the unknown radionuclide distri-bution in each voxel subject to specific inequality constraints over time. In the case of a decreasing only activity, this constraint consists of a continually decreasing amount of activity over time. This means that the activity present in a given voxel at a time tx must be greater than the activity present in the same voxel at a later time t2. That is, Similarly for the case of an increasing activity, we can write the inequality constraint to be, %ih > x i t 2 or x i h - x i t 2 > 0 (4.9) x i h < x i t 2 o r x i t 2 ~ x i h > 0 (4.10) Chapter 4. The dSPECT Method 60 We see that in both cases, we have rewritten the inequality constraint such that the difference between two successive voxel activities is always greater than or equal to zero. Since the unknown voxel activities at each time are represented by a vector, we can write a matrix, A,, that is able to operate on the voxel activities to provide the difference in activity between two successive voxels. Thus between two successive time intervals, the difference in the ith voxel can be written as, / A i A t l ^ A{x.i — At2 Xi,ti X (4.11) Written explicitly, the form of the difference matrix for an individual voxel for the case of a decreasing activity over time will be, Ai = \ (4.12) 1 - 1 0 0 1 - 1 ••. '•• ••• 0 '•• 1 - 1 0 1 In the case of an increasing activity in each voxel over time, a modification exists in the difference matrix such that, V ) Chapter 4. The dSPECT Method 61 1 0 -1 1 o ••• 0 V (4.13) / 0 -1 1 In the final scenerio, the case of combined increasing and decreasing activity for each voxel, the two types of difference matrices are combined and result in another matrix. The determination of the point at which the change occurs between increasing and decreasing activity is the subject of a discussion in a following section. However, once this point is established, it is a straightforward combination of the above matrices in order to allow both increasing and decreasing behaviour for each voxel. Thus, for a combined increasing and decreasing activity, the difference matrix will follow the form, 1 -1 0 0 1 \ • 1 1 0 0 0 0 V -1 0 1 -1 0 (4.14) Here the row at which the upper and lower diagonal elements switch from 0 to -1 and -1 to 0, respectively, represents the position of the peak activity for the ith voxel. Notice that with a difference matrix of this type, once established, the position of peak activity is fixed during the reconstruction process, although this position may be different for each voxel. Without a priori knowledge, there is no way to know this point exactly and Chapter 4. The dSPECT Method 62 so the choice of the peak activity position is only an estimate. A method to estimate this position is described in a later section. Each voxel in the reconstructed image will have its own difference matrix which describes its dynamic behaviour. Since voxels cannot be separated in the image recon-struction process, we must formulate a method to deal with all voxels simultaneously in order to perform the image reconstruction. Thus, we combine each voxel difference matrix together into a tensor given by, A i 0 0 A2 \ V (4.15) ••. '•• 0 0 Ayoxels where the Aj's represent difference matrices for every ith object voxel. Notice that the new matrix A is of dimension, [(# of voxels) (# of proj) x (# of voxels) (# of proj)] Methods to reduce the size of this matrix will be discussed in later sections. Now that all the elements in a dSPECT reconstruction process have been described, it is possible to formulate the tomographic problem in the dSPECT case. Recall from equation 4.6 that the projection of a two dimensional static activity distribution into projection space can be written in terms of a matrix-vector product. In the case of a dynamic object, the activity in each voxel will change over time as the camera moves, but the projection operation can still be described in terms of the same matrix-vector product with one difference. In dynamic imaging, the sizes of the matrix C, and the vector x are now larger to reflect the fact that a different activity is present in each voxel Chapter 4. The dSPECT Method 63 at each projection time. As well, we have seen that the dynamic change in activity from one projection to another can be described in terms of a difference equation with the difference tensor A and the activity vector, x, where, Ax = x (4.16) With this, the projection operation in the dynamic case can be written as, Cx = CA~1 x = d (4.17) Thus, we see that the difficult constrained optimization problem of solving for x subject to the inequalities given in (4.3-4.5) can easily formulated by using the activity differences, x instead as we only need to include the positivity constraint x > 0 into the reconstruction routine. In the next section, further details of the dSPECT reconstruction procedure will be given. 4.4 T H E D S P E C T R E C O N S T R U C T I O N A flowchart of the various steps in a dSPECT reconstruction is given in Figure 4.5 with each section explained in detail in the text. 4.4.1 D E T E R M I N E O B J E C T B O U N D S A typical SPECT reconstruction involves determining the radiotracer distribution in typically 643 or 1283 spatial locations. As a dSPECT reconstruction determines the radionuclide distribution at each point in time in the projection sequence, it follows that the number of unknown variables is multiplied by the number of temporal stops. This increase in the number of variables means that the reconstruction process could become Chapter 4. The dSPECT Method 64 Determine Object Bounds Static Reconstruction eg., OS-EM, WLS, FBP Determination of Peak Activity Position and Dynamic Regions Reconstruction with Increasing Constraint Reconstruction with Decreasing Constraint Specify Shape Constraints Specify Regularizer R(x) -o Reconstruction 0 no Determine model . from aposterior fit V J Reconstruction with Shape Constraints and Regularizer Choose appropriate K, a, P, Y Reconstruction with Linear Difference Method Shape Constraints, and Regularizer Figure 4.5: A flowchart depicting the various steps involved in the dSPECT reconstruc-tion method. very time consuming. Given this, it would be advantageous to only reconstruct those voxels in the object that are known to or assumed to contain activity. It has been seen that when conventional filtered backprojection reconstruction is per-formed on dynamic data, images with poor spatial qualities result. When reconstructed with iterative type methods (eg., expectation maximization, least squares, etc.), it is seen that the amount of streaking artifact is reduced and for the most part, a reasonable spa-tial image is obtained. However, the reconstructed image is not quantitatively correct as Chapter 4. The dSPECT Method 65 a static reconstruction produces an image with voxel values corresponding to the average activity present over the entire acquisition time. If an initial reconstruction is performed with one of these iterative methods, it is possible to apply a thresholding operation in order to eliminate voxels that do not contain a significant amount of activity. Those voxels below the threshold level will then be deemed to contain a negligible amount of activity and therefore can be eliminated from further consideration. Those voxels above the threshold however will be regarded as containing a potentially dynamic behaviour. Using this method to determine the object bounds typically reduces the number of variables by a factor of two thus resulting in reconstruction speed increases on this order. At present the threshold level is operator dependent and is determined by viewing the effect that different threshold levels have on the reconstructed images. It is deemed that a reasonable threshold has been chosen when most of the voxels within the object are still present, while those voxels outside the object are eliminated. As an alternative to the static reconstruction method to determine the object bounds discussed, it is possible to use patient specific attenuation maps instead. By using these maps, only voxels within the attenuating medium are reconstructed in the dSPECT process, again resulting in a decrease of about half the number of variables. This method has the advantage that it is not as operator dependent as the previous method. However, one must be aware of potential artifacts in the reconstructed attenuation maps that could potentially affect the selection of the object bounds. 4.4.2 D E T E R M I N A T I O N O F P E A K A C T I V I T Y P O S I T I O N In the case of only decreasing or increasing dynamic behaviour, the inequality constraint is easily formulated prior to image reconstruction since it is known a priori that all voxels Chapter 4. The dSPECT Method 66 in the object will undergo similar dynamic behaviour. In the case of a combined increasing and decreasing activity however, the situation is not as simple since it is not known prior to image reconstruction at which point the inequality changes from increasing only to decreasing only. One method to determine this point is presented below [70]. Consider the combined increasing and decreasing activity time activity curve for a single voxel as shown in Figure 4.6a). For this voxel, the true position of peak activity occurs at projection number 15 out of 64, therefore, the modified inequality constraint for this voxel should be, 0 < xitl < ... < xitu < xitl5 > xitie > •••> Xi64 > 0 (4-18) Though this position is not known prior to the image reconstruction process, it can be estimated by performing two initial reconstructions on the same dynamic data; one assuming the voxel activity to decrease over the scan time, and the other assuming the activity to increase over the same scan time. After performing these reconstructions, time activity curves can be produced as shown in Figure 4.6b). Notice that when re-constructed using the decreasing only constraint, the time activity curve reconstructs correctly after projection 30, but remains flat over the region 1-29. Similarly the increas-ing only reconstruction yields time frames 1-7 reconstructing correctly as increasing, but with the times 8-64 depicted as static. If the correct portions of these two reconstruc-tions are combined (ie, projections 1-7 and 30-64), a flat region extends over projections 8-29, thus indicating that the true peak activity is situated somewhere in this interval. Therefore, as a first guess, we simply use the point midway on this plateau (projection 18) as the point titP in the reconstruction algorithm. Chapter 4. The dSPECT Method 67 0 10 20 30 40 50 60 0 10 20 30 40 50 60 Projection Number Projection Number Figure 4.6: a) The true time activity curve (TAC) of a single voxel for a dynamic in-creasing/decreasing simulation, b) The reconstructed TAC for the simulation assuming only increasing, and only decreasing activity. 4.4.3 D E T E R M I N A T I O N O F D Y N A M I C A N D S T A T I C V O X E L S In a typical dynamic SPECT study, many of the object voxels will either not contain an appreciable amount of radiotracer, or the amount of tracer that within each voxel will remain approximately constant over the imaging time of the study. As a result, it makes little sense to reconstruct these static voxels assuming dynamic behaviour. Based upon the apparent dynamic information determined from both the increasing only and decreasing only reconstructions, the following selection criteria are used to differentiate dynamic or static voxels. Since the approximate peak activity position has been located for each voxel using the method presented in the previous section, the maximum activity present within each voxel will occur at this position. Thus, we will denote this maximum activity as x i : m a x . It is also possible to determine the minimum activity in the voxel temporal sequence. This will be denoted as XitTJlin. Finally we will consider the average activity present in the voxel over the total number of projections, S, found through the equation, Chapter 4. The dSPECT Method 68 x. t,avg = E -c# (4-19) Due to the inherent Poisson noise present in a voxel, the expected standard deviation in the number of counts within a voxel is in the range ±y/xitavg. We consider a voxel to be static when the difference between the maximum activity and the minimum activity is less than twice the standard deviation. That is, in order for a voxel to be static, the absolute difference in activity within a voxel must obey the condition, Ax{ = Xi}Tnax Xi^min ^ 2iy/Xi^aVg (4.20) Once a voxel is deemed to be static, further modifications are carried out on the difference and system matrices to account for this. For static voxels, the activity present within these voxels at each projection measurement is the same and so rather than solving for a different activity value at each projection stop, it is possible to solve for only one value. In order to solve for this single variable, changes are made to both the difference matrix, Aj, and the dynamic system matrix, Ci, in order to account for this reduction. As the voxel activity at each point in time for a static voxel remains the same, we have, Xitx = xit2 = xit3 = ... = xit, > 0 (4.21) Thus, we recognize that the only unknown variable that has to be solved for is xu3 in the sequence. As a result, the difference matrix for a static voxel reduces to unity (ie., Ai = 1). As well, the system matrix also becomes reduced in size, although the number of nonzero elements remains the same. Whereas the sytem matrix for a dynamic voxel was of size [(# of bins) (# of proj) (# of heads) x (# of proj)}, a static voxel system ma-trix collapses back to the static case of size [(# of bins) (# of proj) (# of heads) x 1]. Chapter 4. The dSPECT Method 69 It it interesting to notice what occurs when the entire object is considered as static. In such a case, alterations are made to the dynamic system matrix, C in order to convert it to the static case, C and the difference tensor, A simply becomes equal to the identity matrix. Once these modifications are carried out, then it is seen that the dSPECT reconstruction process reverts back to the case of a conventional static reconstruction. That is, C A'1 x = C x when A = I = A'1 (4.22) In the implementation of the dSPECT method, the three previous procedures can be combined into a dynamic mask that is used at each stage in the reconstruction process. This dynamic mask consists of a value for each object voxel that describes its dynamic behaviour. That is, if a voxel does not contain significant activity, then the corresponding position in the dynamic mask will contain a value of-1. This signifies to the reconstruc-tion algorithm that it should skip past this voxel when performing the reconstruction as it does not affect the result. If a voxel is determined to remain static, then the cor-responding mask value will contain a 0 at this location. Again, when encountered in the reconstruction process, the algorithm will perform the steps necessary to reconstruct this voxel as a single unknown activity present at each time frame, thus maintaining the static nature and thereby reducing the number of variables to solve. In the final scenario, the case of a dynamic voxel, the voxel mask value will correspond to the projection number where the peak activity position is found. For example, if the dynamic behaviour for a voxel is determined to be decreasing only, then the corresponding location in the mask will contain a value of 1 denoting that the peak activity occurs at projection number 1. For increasing only behaviour, the mask voxel will contain a value corresponding to the number of the last projection (ie., projection 5). When a voxel has Chapter 4. The dSPECT Method 70 2 0 Figure 4.7: An example of the dynamic mask used in dSPECT reconstructions. A value of -1 signifies a null voxel, values of 0 correspond to static voxels, while other values correspond to the positions of peak voxel activity. both increasing and decreasing behaviour, the correponding mask value will contain the position of this peak activity. As it would be expected that each object voxel will have a different dynamic behaviour, it is possible to combine such information in the dynamic mask and reconstruct both static and dynamic voxels simulteously. In Figure 4.7 the true dynamic mask is shown for a simulated object that will be used in further studies. 4.4.4 I N I T I A L O B J E C T E S T I M A T E In many iterative type reconstruction methods, it is important to start with a good estimate of the object activity distribution. By starting with a reasonable starting guess, the reconstruction process is steered away from any local minima in the solution space and this decreases reconstruction time. For the dSPECT reconstruction, a static EM reconstruction is first performed in order to obtain a reasonable spatial image of the object. Thus, this reconstruction contains information about the structure present in the object but without any information relating to the dynamic behaviour. In order to provide an initial estimate of the dynamic behaviour, it is assumed that the activity Chapter 4. The dSPECT Method 71 It4 o < 12 10 * % i \ — Heart 1 - - Heart 2 /•' V It h f 1 1 V \ V \ . . ' r • # XV x> 10 Time 15 20 Figure 4.8: An example of the initial object activity values used in a dSPECT recon-struction. On the left is the static EM reconstruction shown and on the right are the initial time activity curves for the two heart sections for the dynamic mask shown in Figure 4.7. changes within this object in a linear fashion in accordance with the inequality constraints determined for each voxel. That is, once the appropriate inequality constraint is obtained, the dSPECT algorithm assumes that the activity present in each voxel of the object starts initially at half the static EM result, increases in a linear fashion up to the maximum value at the peak position and then is followed by a linear decrease from this point on until again half the static EM result is present at the final projection time (see Figure 4.8). 4.4.5 A D D I T I O N A L S H A P E C O N S T R A I N T S The inequalities shown in equations (4.3)-(4.5) represent shape constraints on the recon-structed time activity curve for a given voxel by limiting the values of the differences in voxel activity over time. It is possible to incorporate further shape constraints in the TA curves through the second differences in the voxel activities. Consider the case of a decreasing activity in a given voxel over time. We have seen Chapter 4. The dSPECT Method 72 that such a dynamic behaviour can be obtained by enforcing a positive difference in the voxel activity over time. That is, a decreasing activity means that the first difference is always positive (ie., Xitl — x^2 > 0, etc). Similarly, for an increasing activity over time, the first difference must always be negative (ie., xitl — xit2 < 0,etc). Consider the case of a convexly increasing activity over time. As the activity is always increasing, then the first differences must always be less than zero (ie., Xitl = X{tl —xit2 < 0, Xit2 = xit2 — xit3 < 0, etc ). The additional convex behaviour can be described by the fact that second differences must always be less than zero as well. That is, x\ ~ x2 = {xitl - xit2) - (xit2 - xit3) < 0 Similarly for concavity, £i - x2 = (xih - xit2) - (xit2 - xit3) > 0 It is also possible to combine both portions of convexity and concavity for instance to maintain convexity on the increasing and decreasing portions, but to allow for concavity near the peak position (Figure 4.9). These additional shape constraints are incorporated into the difference matrix, Ai for each object voxel. 4.4.6 S P E C I F Y R E G U L A R I Z E R Depending on the dSPECT reconstruction algorithm used, it may be possible to incor-porate additional a priori knowledge about the expected object activity distribution into the reconstruction process. This is done by including an additional regularizing term in the reconstruction algorithm. As discussed in the next section, this term may incorpo-rate knowledge of camera response, image smoothness, or additional time activity curve constraints. Chapter 4. The dSPECT Method 73 Convex Increasing Convex Decreasing Figure 4.9: Examples of more stringent inequality constraints to produce increasing and decreasing activities. On the right is depicted an example of a combined increas-ing/decreasing activity with convexity and concavity constraints. 4.4 .7 D Y N A M I C I M A G E R E C O N S T R U C T I O N As mentioned in Chapter 2, image reconstruction involves the minimization of specific objective functions. For the dSPECT method, two different objective functions have been implemented and tested. Both techniques are based on probability theory, but are implemented in different ways and so it has been beneficial to test both algorithms at this stage. The two algorithms used are the constrained least squares algorithm (CLS) which is based on Gaussian probabilities, and the dynamic expectation maximization (dEM) algorithm which is based upon Poisson probabilities. C O N S T R A I N E D L E A S T S Q U A R E S O P T I M I Z A T I O N - C L S The use of mathematical optimization techniques for image reconstruction is common practice in nuclear medicine research. This type of reconstruction usually results in im-proved image quality compared to filtered backprojection as it is able to provide a better physical model of the photon transport through the object (eg., attenuation, scatter, Chapter 4. The dSPECT Method 74 etc), and more accurate models of the camera acquisition process (eg., detector response, etc). In order to perform this type of reconstruction, an objective function based upon a likelihood function is created and then maximized in order to determine the most probable object activity distribution that would give rise to the measured data. Examples of such objective functions include the Gaussian based least squares (4.23), or the Poisson based likelihood function (A.6)(see Appendix A for a derivation of these objective functions). With either objective function, mathematical optimization techniques must be used in order to minimize the objective function and thus to find the best solution for the unknown activity distribution which best matches the actual measured data. Typi-cal optimization methods include conjugate gradient, coordinate descent or the simplex method. The algorithm that has been chosen to perform the optimization process in the dSPECT case is the Limited Memory BFGS-B algorithm [71], which is a quasi-Newton optimization method. This algorithm was chosen as it was designed for use with large scale optimization problems with boundary constraints. In most optimization methods, numerical constraints can be implemented into the minimization procedure in order to limit the values for the estimated variables. The most common constraint used in image reconstruction is to enforce positivity of the reconstructed object voxel values since it is physically impossible for images to contain negative amounts of radioactivity. (4.23) or (4.24) Chapter 4. The dSPECT Method 75 We have already seen that it is possible to express the dynamic SPECT problem in terms of a set of inequality constraints for each object voxel. Therefore by formulating an objective function that relates the voxel activity to the projection data at each time and by implementing the appropriate inequality constraint, it is possible to perform the dSPECT reconstruction. The objective function that has been chosen to be minimized is the Gaussian derived least squares. The appropriate inequality constraints are for-mulated into the difference matrix, A, as described previously and used in the objective function. Thus, the objective function used for the constrained optimization dSPECT reconstruction is, F = 1 {CA-1* - d} 2 (4.25) subject to x > 0 (4.26) Notice that this objective is consistent with (4.23) under the additional assumption that the weighting term, Ojk = 1. In order to speed up the convergence of the algorithm, the gradient is explicitly specified as, VF(x) = (A-1)TCTCA-lk - (A-1)TCTd (4.27) Following the CLS optimization, the variable x represents the difference in activity between successive time frames for each voxel. In order to relate the differences back to the actual amount of activity at each time frame, x is simply multiplied by the inverse difference matrix A'1 (ie., x = A _ 1 x) . One of the main advantages of the CLS reconstruction method is in the fact that regularization terms can be used with the objective function in order to stabilize the reconstruction process. This is especially important in dynamic SPECT image recon-struction since the problem is underdetermined due to the large number of unknown Chapter 4. The dSPECT Method 76 voxel activities and the small number of measured data points. The regularizing term is added to the original objective function such that the new objective function becomes, where a is a parameter through which the user specifies the relative importance between the goodness of fit term (the left term) and the regularizing term, i?(x). As mentioned in Chapter 2, the regularizing term can take various forms depending on the a priori knowledge of the expected object activity distribution. In the case of dSPECT image reconstruction, the selected regularization term would provide not only a spatial smoothing of the transaxial images, but also a temporal smoothing within each object voxel. Such spatiotemporal regularization would tend to produce higher quality images as well as smoother, more realistic time activity curves for each voxel. Different regularization terms for use in dSPECT are being investigated and are the subject of ongoing research [72] and therefore will not be presented in any further detail here. For this work, the term CLS reconstruction will refer only to the Gaussian likelihood fit term without the additional regularization term. T H E DYNAMIC E X P E C T A T I O N MAXIMIZATION A L G O R I T H M - D E M As an alternative to the CLS reconstruction algorithm which relies on an external math-ematical optimization routine, other algorithms have been investigated to handle the linear inequality constraints. One reconstruction algorithm that is becoming popular in nuclear medicine image reconstruction is the expectation maximization (EM) algorithm or variations of it such as ordered subset-expectation maximization (OS-EM) [73]. The EM algorithm has advantages over other reconstruction methods because it has a relatively rapid convergence rate, allows for the inclusion of complex photon transport (4.28) Chapter 4. The dSPECT Method 77 functions such as nonuniform attenuation and scatter compensation, and is well suited for dealing with geometrical effects such as two or three dimensional collimator blurring. The static E M algorithm has been presented in Chapter 2 and is discussed in further detail in Appendix A, therefore it will not be discussed further here, but rather the dynamic EM (dEM) algorithm will be presented. In the static case, at each iteration of the EM algorithm, the updated estimate of the amount of activity in the ith voxel can be found through the equation, „new _ x i ^j,k^i]k"jr» (A OQ\ x°ld Ej,k Cijkdjk Ejk Cijk Ez' Ci'jkX? In the dynamic case, the dEM algorithm can be written in a similar form to the static case but rather than working with the voxel activity values X{ directly, the voxel differences, xu~ are used instead. In addition, the static projector matrix, C is replaced by the dynamic system matrix and difference matrix combination CA'1. As a result of these changes, the dEM algorithm has the form, -new _ ^ik ^3 ^^3^ HjkyjK ^ ik ^ fi.... A-lz.old ^ ' ' EjiCijkA^l) E i ' Ci'jk iftxtfji where xfkew represents the updated estimate for the difference in activity in the ith voxel in the A;th time interval. As in the case of the CLS reconstruction, in order to obtain the voxel activity values x, following the dEM reconstruction, the difference vector x is simply multiplied by A"1. In the static case, it can be shown that the EM algorithm maximizes the Poisson like-lihood function. It can further be shown that in the absence of noise, the EM algorithm converges to a unique solution [25]. However, when noise is present, the E M algorithm will still always maximize the likelihood function but with high levels of noise, this is not always a desirable property. In fact, after just a few iterations, it has been shown Chapter 4. The dSPECT Method 78 that the EM algorithm starts to fit the noise in the data and as a result, diverges from reconstructing the true activity distribution in the object [74] [75] [27]. Because of this property, it is unclear at which point the EM reconstruction procedure should be stopped [76]. In the dynamic EM algorithm, the same pitfalls are present, although in the dynamic case it would seem that the number of iterations is not as important as for the static case. Presumably this is because the number of variables is several times larger in the dynamic case compared to the static case and so the error produced by performing too many iterations tends to get spread out over all the variables in time, rather than just at one time frame as in the case of static reconstruction. Thus, the error in any one voxel in time remains relatively small. One of the main advantages of the dEM algorithm over the CLS algorithm is in the fact that the positivity constraint required in the CLS method is intrinsic in the dEM algorithm. This is because at each iteration of the dEM algorithm, the current guess of the voxel activity gets multipled by a correction factor. This multiplicative factor must always be greater than 0 since the rightmost denominator in the dEM algorithm (known as the projection step) is always greater than 0 provided that the initial guess of the voxel activities was greater than 0. 4.4.8 a posteriori F I T Once dynamic projection data has been reconstructed with either the CLS or dEM method, the result is a sequence of images which represent the radionuclide distribu-tion within the object at the midpoint time of each projection. Further analysis can then be carried out identically to that of dynamic PET studies where three dimensional re-gions of interest are drawn and time activity curves are produced in order to fit to specific Chapter 4. The dSPECT Method 79 compartmental models. As no compartmental modelling was performed in this research project, no further discussion of these a posteriori fitting procedures will be given here. V A R I A B L E P E A K P O S I T I O N As seen from the difference matrix for a given voxel, (4.14), the position of the peak activity is very important in the reconstruction process while the method used to deter-mine this point is only approximative. Therefore, in order to relax the constraint in the vicinity of the peak, we have chosen the following method in which the position of peak activity is free to move within a selected time frame to one of three possible positions. This is done through the use of the revised constraint, Xitv > (xitp-i + xitp+i) (4.31) which modifies the inequality difference matrix for each voxel into the form, 1 0 -1 1 \ - 1 2 -1 0 1 V 0 With this relaxed constraint, the position of peak activity is now allowed to be in any of the positions xitp_1, xitp or xitp+i so long as the activity at time p is greater than the mean of the surrounding points (see Figure 4.10). This modified difference matrix has been utilized in all reconstructions of temporal behaviour consisting of a combined increasing and decreasing activity. Chapter 4. The dSPECT Method 80 Figure 4.10: The three possibilites of the peak activity position when used with the modified difference matrix above. 4 .5 C H A P T E R S U M M A R Y In this chapter, the general framework of the dSPECT reconstruction method has been discussed. Two reconstruction algorithms have been presented that are able to include the inequality constraints necessary for the dSPECT reconstruction. As well, many specific implementation details that either speed up the reconstruction process or result in enhanced reconstruction accuracy have also been shown. For the remainder of the work, the constrained least squares (CLS) and the dynamic expectation maximization (dEM) algorithm will be tested for accuracy in dSPECT reconstructions. C H A P T E R 5 CO M P U T E R SIMULATIONS This chapter presents a series of validations of the two dSPECT reconstruction methods using simulated projection data of a dynamic source distribution. Two different objects have been simulated; one to test the feasability of the dSPECT method with different camera acquisition protocols and dynamic parameters, and the other to test the limits of accuracy of the reconstruction under near clinical circumstances. For each simulation, the accuracy of the resulting dSPECT reconstruction was assessed using three different criteria. The purpose of each of these criteria will be discussed and the results reported for each simulation. 5.1 S I M U L A T I O N M E T H O D S In order to test the dSPECT reconstruction method for accuracy, various computer simulations have been performed. These simulations can be classified as of one of two types: analytical or discretized. The first set of simulations were of the analytical type and modelled a 2D annulus subdivided into four equally sized regions (Figure 5.1). This particular shape was chosen as it models a typical short axis slice through the myocardium. This annulus phantom has 81 Chapter 5. Computer Simulations 82 been used extensively to investigate different acquisition and reconstruction parameters such as the optimal number of camera heads or the number of projection measurements needed for accurate dSPECT reconstruction. Figure 5.1: A diagram of the annulus phantom used in the first set of simulations. Projection data was collected into 64 transaxial camera bins over a 20 minute acquisition time with various camera geometries. To create data simulating dynamic projection measurements, it is necessary to per-form a projection of the annulus into discretized camera bins at various angles corre-sponding to different camera angles. As each region is simulated with a different dynamic behaviour, it is necessary to perform this projection on each region individually. After all the projection operations have been performed the final step is to sum together these individual projection data sets in order to create projection data for the entire object. For each annulus segment, a different dynamic behaviour was assigned in the form of a time varying function (eg., mono exponential decay, dual exponential decay, etc). At each angle, a measurement was simulated by projecting the respective annulus segment into the appropriate camera bins and then scaling it by a factor corresponding to the value of the time dependent function at this time. As mentioned, the projection measurements for each annulus segment were created Chapter 5. Computer Simulations 83 independently and then simply added together to yield projection data for the entire annulus. Additional Poisson distributed noise was further introduced to the projection data at this point. The projection operation has been performed assuming perfect parallel hole collimation and no provision has been made for object attenuation. The second type of simulation more realistically models data expected during clinical acquisition procedures. These simulations will be referred to as discretized simulations as they use an object discretized into small voxels and allow for modelling of such effects as statistical noise, photon attenuation and camera response. The model used for these simulations consisted of a cross sectional slice of a thorax containing a myocardial region, lung tissue, and spinal column (Figure 5.2). Tests were performed using this phantom with clinically relevant levels of activity present in all regions. Dynamic behaviour was included only in the myocardial region as this is the only region that would be expected to display significant dynamic behaviour. The inclusion of the lung tissue and spinal column created a nonuniform attenuating medium and a uniform, static background in each region. 0.02 0.018 0.016 0.014 0.012 0.01 0.008 0.006 0.004 0.002 Figure 5.2: Left) Activity distribution map shown at t=5 min. Right) Nonuniform attenuation map phantom used in the attenuated simulations. Attenuation coefficients shown are in units of mm" 1. Chapter 5. Computer Simulations 84 Projection measurements were created by modelling the forward projection operator for a parallel hole collimator with nonuniform attenuation at each time interval and projection angle in the acquisition sequence. 5 . 2 D A T A A N A L Y S I S M E T H O D S In the case of computer simulations, the true object activity distribution is known at each point in time for each object voxel. This enables one to assess the accuracy of image reconstruction quantitatively. For dynamic S P E C T image reconstruction, there are two different criteria that must be evaluated for every reconstruction. The first criterion measures the quantitative accuracy in the reconstructed activity within each object voxel (or region) while the second measures the accuracy of the reconstructed time activity curve for a selected region (ie., how well the reconstructed dynamic function parameters compare with the true dynamic function parameters). In order to evaluate the quantitative accuracy of the reconstructions, it is possible to measure the deviation of the reconstructed activity values within each reconstructed region from the true activity values. This measurement will then refer to how well the reconstructed images compare to the true images over all time frames. To evaluate the second criterion, it is sufficient to measure the difference between the true dynamic functional parameters used in the simulations, and those obtained after dSPECT reconstruction and a subsequent time dependent function fit. These assessment techniques will be discussed further in the next sections. 5.2.1 R E L A T I V E S T A N D A R D D E V I A T I O N Since a dSPECT reconstruction produces a series of images with a different activity in each voxel at each time, in order to measure quantitative accuracy, it has been chosen to Chapter 5. Computer Simulations 85 perform an evaluation of each time image separately and then to combine these individual errors into a final reconstruction accuracy value. The accuracy value used is a measure of the relative standard deviation. The simu-lated amount of activity within each region is distributed uniformly and thus for a given dynamic region in an object, all voxels should ideally reconstruct to the same value. How-ever, due to factors such as statistical noise, in reality each voxel will typically reconstruct a value that will be slightly different than neighbouring voxels. One can measure this variation of reconstructed voxel activities by measuring the standard deviation of the reconstructed voxel activities to the true voxel activities for each region. The equation used to evaluate this is, where x^ is the true activity at each time frame and ok is the standard deviation within the region calculated as, where Xik is the reconstructed voxel activity within the ith voxel at the kth. time interval, total number of voxels within each region (208 for each region in the annulus phantom.) With this accuracy measurement, it should be evident that values closer to 0 represent reconstructions with less deviation from the true activity while larger values indicate a large degree of deviation within the reconstructed region. In order to evaluate the reconstruction as a whole, the scalled L2 — norm of the regional accuracy values over the four dynamic regions in the phantom is reported. That is, eR = (5.1) (5.2) Xik is the true simulated activity in the same voxel at the same time interval and NR is the Chapter 5. Computer Simulations 86 eT - (5.3) 5 . 2 . 2 K I N E T I C P A R A M E T E R D A T A ANALYSIS In all simulations performed, known time dependent functions were modelled. These functions consisted of either mono exponential decay, dual exponential decay, or a dual exponential function comprised of both an increasing and a decreasing portion. Following reconstruction, it is possible to perform an accuracy measurement by obtaining time activity curves over each dynamic region and fitting these time activity curves to the modelled time dependent function (eg., mono exponential washout, etc). Once fit in this way, dynamic parameters of the reconstructed T A curves can be compared to the true dynamic parameters modelled in the simulation. Such parameter fits were carried out by minimizing the Y 2 value between the estimated function value and reconstructed TA data for a given region of interest. With more complicated multiple exponential models, difficulties arise in possible nonuniqueness of solutions when using this type of accuracy assessment. It will be seen in some reconstructed T A curve parameter fits, that the estimated fit parameters do not agree well with the true dynamic parameters specified in the simulations. In such cases, it falls upon the reader to subjectively evaluate the accuracy of reconstruction by means of a comparison between the reconstructed TA curve and the true simulated T A curve and by means of the SR value previously described. 5 . 2 . 3 P E A K A C T I V I T Y POSITION E R R O R ANALYSIS In simulations involving both activity increase and decrease for a given voxel, a further analysis was used in order to quantify the accuracy of the peak finding algorithm. Since Chapter 5. Computer Simulations 87 all voxels within a given dynamic region were simulated with identical time activity curves, it follows that the position of peak activity should also be identical for each voxel within the region. In the reconstructed images however, each voxel has a slightly different peak activity position due to the inexact peak determination method. In order to assess this accuracy, the standard deviation between the estimated peak activity position, pi, and the true position of peak activity, pi, is calculated for each region R as, where NR is the number of voxels in the dynamic region. 5.3 A N A L Y T I C A L S I M U L A T I O N S In the annulus phantom simulations, acquisition protocols were simulated which modelled the various camera geometries possible with the Siemens line of nuclear medicine cameras (ie, single, dual, and triple head systems). In order to maintain a consistency among all simulations, identical angular sampling intervals of « 3° were used for all studies. As a result, for studies involving a rotation of 180° for each head, 64 projections have been acquired per head, thus resulting in 64 reconstructed dSPECT images, while those with less rotation (ie, triple head -120° and two head - 90°) contain fewer reconstructed images (40 and 32 frames respectively). All simulations were acquired over a total rotation time of 20 minutes resulting in the identical number of counts per head for all simulations of the same dynamic object. Acquisition parameters are summarized in Table 5.1. In order to obtain a measurement of the ideal reconstruction accuracy in the absence of any complicating factors, for each camera geometry, two simulations were performed using the same kinetic parameters; one without statistical noise in the projection data (5.4) Chapter 5. Computer Simulations 88 Acquisition Head Rotation Projections Bins Time Type Geometry (per head) (per head) (minutes) 1) Single head n/a 180° 64 64 20 2) Dual head 180° 180° 64 64 20 3) Dual head 90° 90° 32 64 20 4) Dual head 90° 180° 64 64 20 5) Triple head 120° 120° 40 64 20 6) Triple head 120° 180° 64 64 20 Table 5.1: Summary of simulation acquisition parameters used for all simulations. and the other adding a statisical noise to the simulated projection data in accordance with the level of activity present within each camera bin. By comparing accuracy values for both sets of reconstructed images, it is possible to estimate the extent to which additional statistical noise in the projection data corrupts the image reconstruction process. 5.3.1 M O N O E X P O N E N T I A L W A S H O U T As the first test of the dSPECT 'reconstruction method, simulations were performed using the annulus phantom with a different mono exponentially decreasing radioactivity distribution in each dynamic segment. For each region, the activity followed a time dependent function of the form, A(t) = A0 e~x 1 (5.5) Different kinetic parameters were present for each section, but initial activities were the same, corresponding to a total of 1.85 MBq in each dynamic section. The camera efficiency was modelled as 0.01%. A summary of the parametric values used in these simulations is shown in Figure 5.3 and Table 5.2. Note that T i / 2 = j^p. Chapter 5. Computer Simulations 89 o 3 "0 5 10 15 20 Time (min) Figure 5.3: True annulus object shown at t = 5 min for mono exponential decay simula-tion. Section A0 (cnts/voxel) A (min - 1) T\/2 (min) 1 57.8 0.35 2.0 2 57.8 0.17 4.0 3 57.8 0.087 8.0 4 57.8 0.043 16.0 Table 5.2: True annulus phantom dynamic parameters used in simulations of mono exponential decay. Chapter 5. Computer Simulations 90 RESULTS Simulations were reconstructed with both the dynamic constrained least squares method (CLS) and the dynamic expectation maximization (dEM) algorithm. A summary of accuracy values is presented in Table 5.3. While a total of six different camera acquisition protocols were simulated, because of space constraints, only three of the reconstructions are shown in expanded form in Appendix B. These consist of a single head 180° rotation, a dual head 90° configuration - 180° rotation per head, and a triple head 120° configuration - 180° rotation per head. Reconstruction accuracy measurements are given for both the noiseless simulations as well as for those involving statistical noise in the projection data. Acquisition Type Overall Error ET Noiseless With Noise CLS Method dEM Method CLS Method dEM Method 1) Single Head 0.41 0.46 0.46 0.46 2) Dual Head 0.41 •0.46 0.44 0.46 3) Dual Head 0.32 0.33 .0.35 0.34 4) Dual Head 0.22 0.26 0.30 0.27 5) Triple Head 0.21 0.24 0.28 0.25 6) Triple Head 0.19 0.23 0.28 0.24 Table 5.3: Summary of overall accuracy values for the mono exponential decreasing activity simulations. In dSPECT reconstructions performed using a single SPECT camera head and a slow rotation, it can be seen that rather significant streaking artifacts appear, primarily at early time frames and this is reflected in the various accuracy measurements. As well, it is apparent that regions undergoing a faster activity washout reconstruct with more error than regions which contain a slower varying activity distribution. Though not apparent from the images reproduced here, when viewing a movie of all the reconstructed frames, streaking artifacts that appear in each time frame are consistent with the location of the Chapter 5. Computer Simulations 91 camera at that projection time. Differences between the CLS and dEM methods can be seen in both time activity curves and in image quality where the dEM result appears to provide smoother images (less noise) as well as smoother TA curves for each dynamic region. Additionally, it is seen in the time activity curves of the CLS reconstruction method, that a sharp drop off occurs at the tail of the TA curves. It is not known at this point why this effect is present. It can be seen both in image quality and with the various accuracy measurement techniques, that acquisitions employing a 90° dual head SPECT protocol provide much more accurate reconstructions than single head acquisitions. Additionally, it appears that the dEM algorithm improves upon the CLS method by providing both smoother spatial images and resultant time activity curves as well as slightly lower error values when the effects of noise are included in the projection data. Reconstructions involving a triple head camera provide the best accuracy values and superior image quality. However, the triple head system appears to produce only marginally better accuracy values compared to the dual head system when acquiring data over a 180° rotation per head. In all cases, it is apparent that the inclusion of noise in the projection data results in a degraded reconstructed image sequence as evidenced by the smaller error value in the noisefree simulations. 5 . 3 . 2 D U A L E X P O N E N T I A L W A S H O U T As an extension of the mono exponential dynamic behaviour presented above, a modifi-cation was made to simulate dual exponential decreasing activities within each annulus segment. This was done in order to evaluate the effectiveness of the dSPECT method in differentiating multiple decay curves simultaneously. Again, camera acquisitions used the same geometrical configurations as in the previous case. Chapter 5. Computer Simulations 92 70 60 50 140 c = 30 < 20 10 0 Section 1 Section 2 — Section 3 • • Section 4 • \ -OO-r.-\ * ^ * \ \ - -• --10 Time (min) 15 20 Figure 5.4: True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation. Decay curves in this set of simulations consisted of both fast and slow components simulating the expected type of curves present in a myocardial fatty acid study. Again , the annulus was divided into four different segments, but this time wi th a decay function for each segment of the form: A(t) = A0 {e-Xlt + e~X2t) The decay parameters for each section are shown in Table 5.4. (5.6) Chapter 5. Computer Simulations 93 Section A0 Ai A 2 Tl/2,2 (cnts/voxel) (min - 1) (min) (min - 1) (min) 1 28.9 0.35 2.0 0.035 20.0 2 28.9 0.17 4.0 0.035 20.0 3 28.9 0.087 8.0 0.035 20.0 4 28.9 0.043 16.0 0.035 20.0 Table 5.4: True dynamic parameters used in simulations involving dual exponential decay in the annulus phantom. R E S U L T S Similar to the case for mono exponential decreasing activity, both CLS and dEM re-constructions were performed and assessed with a summary presented in Table 5.5 and. extended results in Appendix B. As in the case of mono exponential washout, reconstruc-tions with the dynamic behaviour modelling a dual exponential function show that triple head acquisitions yield superior accuracy results compared to both single and dual head acquisition methods. Again however, the differences between the dual and triple head acquisitions is marginal. Acquisition Type Overall Error eT Noise Free With Noise CLS Method dEM Method CLS Method dEM Method 1) Single Head 0.20 0.20 0.25 0.2.1 2) Dual Head 0.20 0.20 0.23 0.20 3) Dual Head 0.19 0.19 0.22 0.20 4) Dual Head 0.16 0.19 0.21 0.20 5) Triple Head 0.15 0.19 0.19 0.20 6) Triple Head 0.15 0.19 0.20 0.20 Table 5.5: Summary of overall accuracy values for dual exponential decreasing simulations using the annulus phantom. Chapter 5. Computer Simulations 94 It is seen in some circumstances (eg. annulus sections 1 and 4 of the dual head dEM reconstruction), that a dual exponential fit to the reconstructed time activity curve returns dynamic parameters that are very different from the true parameters. This is due to the fact that often when fitting multiple exponential functions, there may be a large number of widely different function parameters that are able to describe the observed data. In these circumstances, it is necessary to evaluate the accuracy based upon other criteria such as the e value, or qualitatively from the observed time activity curves. 5.3.3 D U A L E X P O N E N T I A L WITH S L O W WA S H I N AND W A S H O U T In the final stage of analytical simulations, the most clinically relevant case of increasing activity followed by decreasing activity was tested. With these simulations, the same data acquisition protocols were used as previous. For this case, the dual exponential equation simulated was of the form, A(t) = A0 {e~A l t - e-A2'} (5.7). Simulated dynamic parameters are given by Table 5.6. Section A0 Ai T V 2 , 1 A2 Ti/2,2 Peak Position (cnts/voxel) (min - 1) (min) (min - 1) (min) (frame no.) 1 28.9 1.39 0.5 0.035 20 9.0 2 28.9 0.69 1.0 0.035 20 15.0 3 28.9 0.35 2.0 0.035 20 24.0 4 28.9 0.17 4.0 0.035 20 38.0 Table 5.6: True dynamic parameters used in simulations involving a combined increasing and decreasing dual exponential function. Chapter 5. Computer Simulations 95 30 25 0) 20 2 115 > 3 10 — Section 1 Section 2 —— Section 3 • Section 4 1 ' ' j 1 / .... ' / ' / j / ' / ' ' . • ' / / / 111 .• 17 * : * 10 Time (min) 15 20 Figure 5.5: True image at t = 5 min along with time activity curves for each dynamic region for the slow dual exponential simulation. 5.3.4 D U A L E X P O N E N T I A L WITH F A S T W A S H I N AND W A S H O U T A second simulation was also performed with combined increasing/decreasing behaviour but this time wi th much faster dynamic rates. A summary of this second simulation is shown in Table 5.7. Section A0 (cnts/voxel) A i ( m i n - 1 ) Ti/2,1 (min) A 2 ( m i n - 1 ) Tl/2,2 (min) Peak Position (frame no.) 1 28.9 1.39 0.5 0.35 2.0 5.0 2 28.9 0.69 1.0 0.17 4.0 9.0 3 28.9 0.35 2.0 0.086 8.0 15.0 4 28.9 0.17 4.0 0.043 16.0 35.0 Table 5.7: True dynamic parameters used in simulations involving a combined increasing and decreasing dual exponential function wi th fast dynamic components. Chapter 5. Computer Simulations 96 | 1 2 110 & i 8 Section 1 — Section 2 . — Section 3 • Section 4 . j; f \ \ % . ; '/..-.. V \ ' / • \ V I'/- \ s x I • : \ ^ 7 X ^ ^ — ^ 10 15 Time (min) 20 Figure 5.6: True image at t = 5 min along with time activity curves for each dynamic region for the rapid dual exponential simulation. RESULTS The simulations presented here reflect the most difficult situation as they involve both increasing and decreasing activity and the point where the constraints change is not known a priori. Thus, these reconstructions test not only the more complex inequality constraint, but also the method used to determine the peak activity position for each object voxel. Results of these reconstructions are again summarized in Tables 5.8 and 5.9 and in extended form in Appendix B . It is seen in al l accuracy measurements that the d E M algorithm reconstructing the data acquired with a triple head acquisition comes out ahead of the other reconstruction and acquisition methods. In fact, the triple head d E M algorithm was even able to reconstruct accurate time activity curves for Section 1 of the fast uptake, fast washout model. This reconstruction also shows the best agreement with the true peak position and shows relatively litt le deviation in the peak activity position Chapter 5. Computer Simulations 97 (peak position = 5.6 ± 2.2), thus indicating that each voxel gets reconstructed with a very similar time activity curve. Acquisition Type Overall Error eT Noise Free With Noise CLS Method dEM Method CLS Method dEM Method 1) Single Head 0.21 0.32 0.43 0.35 2) Dual Head 0.21 0.32 0.31 0.33 3) Dual Head 0.22 0.29 0.31 0.31 4) Dual Head 0.35 0.25 0.35 0.28 5) Triple Head 0.19 0.24 0.32 0.27 6) Triple Head 0.17 0.23 0.34 0.26 Table 5.8: Summary of overall accuracy values for slow washin / washout simulations. Acquisition Type Overall Error ST Noise Free With Noise CLS Method dEM Method CLS Method dEM Method 1) Single Head 0.81 0.53 0.73 0.56 2) Dual Head 0.63 0.53 0.71 0.53 3) Dual Head 0.38 0.37 0.48 0.39 4) Dual Head 0.61 0.31 0.72 0.37 5) Triple Head 0.35 0.30 0.48 0.34 6) Triple Head 0.81 0.31 0.92 0.37 Table 5.9: Summary of overall accuracy values for fast washin / washout simulations. 5.4 D I S C R E T I Z E D S I M U L A T I O N S Based on results of the annulus phantom simulations, discretized thorax simulations were performed using the same camera geometries as for the dynamic annulus simulations. A summary of thorax acquisition parameters is given in Table 5.10. Again, a variety of Chapter 5. Computer Simulations 98 different dynamic parameters were simulated ranging from a mono exponential decreasing activity, to a dual exponential combined increasing / decreasing function. In this set of experiments, only simulations including Poisson noise are presented as they correspond to a more realistic scenario when object attenuation is included. Acquisition Head Rotation Projections Bins Time Type Geometry (per head) (per head) (minutes) 1) Single head n/a 180° 64 64 20 2) Dual head 90° 180° 64 64 20 3) Triple head 120° 180° 64 64 20 Table 5.10: Summary of the acquisition protocols used in discretized thorax simulations. 5.4.1 M O N O E X P O N E N T I A L W A S H O U T In the first dynamic thorax simulation, two dynamic regions of the heart were modelled with a dynamic behaviour corresponding to a different mono exponential washout for each section. All surrounding tissue (lung, thorax, spine) was modelled as a static background with different amounts of activity concentration. Table 5.11 summarizes the parameters simulated in this case. Thorax Section A0 (cnts/voxel) Ti/2 (min) heart 1 101.3 4.0 heart 2 101.3 8.0 lungs 0.3 oo thorax 6.1 oo Table 5.11: Summary of the true dynamic parameters used in mono exponential decay simulations of the thorax phantom. Chapter 5. Computer Simulations 99 0 5 10 15 20 Time (min) Figure 5.7: True image at t = 5 min along with time activity curves for each dynamic region for the mono exponential decay simulation. R E S U L T S Reconstructed results and images are presented in Appendix B and summarized in Table 5.12. As in the case of the annulus simulation, dSPECT reconstruction using data acquired with a single head camera system does not perform well enough to provide good spatial images or accurate time activity curves. It is also clear that when fast washout rates are present, the reconstruction routine tends to underestimate the amount of activity within each voxel and prolong the time of washout. Again severe streak artifacts are seen in the direction of the camera at each time frame as in the case of the annulus simulations. Results are similar for both the CLS and dEM methods in simulations of dual and triple head acquisitions. It is worthwhile noting that while the dynamic parameter fits more closely resemble the true dynamic parameters in the CLS reconstructions, the ER values are lower in the dEM result, thus indicating that the CLS method is more Chapter 5. Computer Simulations 100 susceptible to large variations in the voxel activities (an effect that is seen in the speckled nature of the CLS reconstructed images). Acquisition Overall Error e r Type CLS Method dEM Method 1) Single Head 2.78 2.16 4) Dual Head 2.05 1.38 6) Triple Head 1.78 1.14 Table 5.12: Summary of overall accuracy values for mono exponential decreasing thorax simulations. (Accuracy values appear large due to large errors in the low count lung regions.) 5.4.2 D U A L E X P O N E N T I A L W A S H O U T Similar to the annulus simulation, another set of simulations were performed modelling different dual exponentially decreasing activities in each of the heart sections. The dy-namic parameters consisted of both fast and slow components as expected for myocar-dial fatty acid imaging. Again, surrounding tissue (lungs, thorax, spine) contained static background activity and provided a nonuniform attenuating medium. The same acqui-sition protocols were used as previously. Chapter 5. Computer Simulations 101 Figure 5.8: True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation. Thorax A0 ^1/2,2 Section (cnts/voxel) (min) (min) heart 1 50.6 4.0 12.0 heart 2 50.6 8.0 20.0 lungs 0.2 oo oo thorax 2.2 oo oo Table 5.13: Summary of the true dynamic parameters used in dual exponential decay simulations of the thorax phantom. — Heart 1 H«art 0 — Thorax V \ \ • s V x \ N* V x • • x• *vi: < 10 Time (min) 15 20 Chapter 5. Computer Simulations 102 RESULTS Results of dual exponential decreasing activity appear very similar to those of the mono exponential decay and so will not be discussed in detail. A summary of the accuracy values is again shown in Table 5.14 as well as in Appendix B. Acquisition Overall Error ET Type CLS Method dEM Method 1) Single Head 3.06 2.10 4) Dual Head 2.25 1.67 6) Triple Head 1.90 1.20 Table 5.14: Summary of overall accuracy values for dual exponential decreasing thorax simulations. 5.4.3 D U A L E X P O N E N T I A L WITH WASHIN AND W A S H O U T The final set of simulations performed with the thorax model included the effects of attenuation with a dual exponential function modelling a rapid activity uptake in the myocardium, followed by a slow washout. Figure 5.9 and Table 5.15 summarize the parameters used in this final simulation. Section A0 Tl/2,1 Tl/2,2 Peak Position (cnts/voxel) (min) (min) (frame no.) heart 1 50.7 1.0 8.0 11.0 heart 2 50.7 2.0 10.0 19.0 lungs 0.3 oo oo 0 thorax 6.1 oo oo 0 Table 5.15: Summary of the true dynamic parameters used in dual exponential washin/washout simulations of the thorax phantom. Chapter 5. Computer Simulations 103 — 2 1 . "0 5 10 15 20 Time (min) Figure 5.9: True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential increase and decrease simulation. RE S U L T S The results of these final simulations confirm that the dEM algorithm with a triple head acquisition protocol outperforms the CLS reconstruction method and other acquisition protocols tested (see Table 5.16). Acquisition Overall Error e Type CLS Method dEM Method 1) Single Head 6.22 3.32 4) Dual Head 2.92 1.50 6) Triple Head 2.21 1.11 Table 5.16: Summary of overall accuracy values for dual exponential increasing and decreasing thorax simulations. As seen from these reconstructions, the peak determination algorithm also performs quite well for the triple head acquisition and dEM algorithm. In particular, it is able to Chapter 5. Computer Simulations 104 detect that the activity distributions in the lungs and thorax remain essentially static, while being able to find the peak heart positions for each section to within a few projec-tion times. Notice however that the time activity curves for the two heart sections in this reconstruction are quite ragged, owing to the fact that the peak position has been found to vary within each section. At this point it is not apparent why the CLS triple head reconstruction has undergone a marked decrease in accuracy compared to other simula-tions in which it performed very similar to the dEM algorithm. However the accuracy values do not seem to reflect the fact that the general shape of the reconstructed TA curve appears very similar to the true TA curve and is even able to reconstruct the point of crossover of the two heart sections quite accurately. 5.5 C H A P T E R S U M M A R Y In this chapter, validation of the dSPECT reconstruction method has been carried out through a collection of computer simulations. Simulations have progressed from first modelling only a mono exponential washout with four large segments of an annulus up to the final simulations which modelled a clinical situation involving nonuniform attenuation effects and similar dynamic behaviour in both heart sections. As well, more complicated simulations have modelled varying dynamic behaviour up to the point of a fast washin and washout without a priori knowledge of the peak activity position for each object voxel. Results of most dSPECT reconstructions indicate that a triple head acquisition with the dEM reconstruction routine yields the most accurate reconstructions. This is closely followed by dual head acquisitions with rotations of 180° per head. Based upon these sim-ulations, it appears that single head systems will not work adequately with the dSPECT method as substantial artifacts are present in reconstructed images. C H A P T E R 6 P H A N T O M A N D VO L U N T E E R EXPERIMENTS In conjunction with the previously described dynamic computer simulations, a series of experimental tests of the dSPECT reconstruction method was performed. In this chapter, a description of how these experiments were performed will be given and reconstruction results summarized. Also presented in this chapter will be a preliminary test of the dSPECT method applied to the clinical application of dynamic renal imaging. 6.1 T H E D Y N A M I C H E A R T - I N - T H O R A X P H A N T O M In order to test the dSPECT reconstruction method with projection data closely resem-bling that anticipated in a clinical dynamic application, we have developed a dynamic heart-in-thorax phantom [77] [78]. This phantom consists of the following dynamic com-ponents: • Central blood pool volume (32 ml) • Myocardial volume (350 ml) 105 Chapter 6. Phantom and Volunteer Experiments 106 • Up to 3 defect volumes spaced equally throughout the myocardial volume (17 ml each) The nonuniform attenuating thorax consists of the following static components: • Tissue equivalent medium (10 liters of water) • Spinal column (1 cm radius aluminum rod) • Two lung tissue medium (each 1 liter sponge) A schematic of the dynamic heart-in-thorax phantom is shown in Figure 6.1. Threaded Rod SpinaL Column „ Threaded Rod Heart Phantom Figure 6.1: The dynamic heart-in-thorax phantom shown with the surrounding hetero-geneous attenuating material. The dynamic heart insert can be seen located centrally in the thorax. Chapter 6. Phantom and Volunteer Experiments 1 0 7 6 . 1 . 1 D Y N A M I C H E A R T INSERT The heart insert consists of a total of five different potentially dynamic volumes. Within each volume, a radioactive compound is injected and gets subsequently diluted with a constant flow of pure water. In order to supply the inflowing water into each dynamic container, each volume contains both an inlet and outlet tube. The rate of change in the amount of radioactivity present in each volume can be described by the first order differential equation, dA(t) = Ajt) R dt V 1 j where A(t) corresponds to the amount of radioactivity present in the volume, V, at time t, and R is the rate of the incoming water supply. Solving for A(t) yields the familiar equation, A(t) = A0 e-v* (6.2) A common parameter used to describe the rate of change of the volume activity is the washout half-life, T i / 2 , which is calculated as, rp ln{2) V Ti/2 = — ^ — (6.3) In order to ensure uniform mixing of the radioactivity with the incoming water supply, within each dynamic volume is contained a small propeller (see Figure 6.2). The propeller mixing systems are further interconnected via a gearing system that is driven by an electric motor (Figure 6.3). In order to avoid any electrical interference between the electric motor and the sensitive SPECT camera electronics, the motor is situated two meters from the camera heads at the opposite end of the patient bed and connected via a flexible drive shaft to the central mixing unit. Chapter 6. Phantom and Volunteer Experiments 108 Figure 6.2: Cutaway diagram of the mixing apparatus used to ensure uniform mixing of the inflowing water with the radioactivity. In order to maintain a constant rate of inflowing water, a digitally controlled syringe infusion pump is used. This pump is capable of delivering accurately controlled flow rates of between 0.0001 and 4000 ml/hr from up to ten 60 ml syringes simultaneously. We have seen that if radioactivity is initially present in a volume at time t = 0 and diluted with the constant flow of pure water, the amount of radioactivity present at any time can be described by a mono exponential decreasing function. If however, the radioactivity is introduced at a time later than t = 0, the result will be an increase in activity initially, followed by a mono exponential decrease. If we consider that the radioactivity that is injected is a very short bolus in time, then we can approximate it as a delta function at time T, 8T- Since a constant flow of water in a volume creates a mono exponential washout over time, the introduction of a delta function of activity (the bolus) into a mono exponentially decreasing washout (the system response) results in a convolution operation. That is, the activity present at any Chapter 6. Phantom and Volunteer Experiments 109 Effluent Container Moto r \ 1 'hantom ~ Driveshaft • — J Patient B e d / Camera Head I I. I 1 f 1 f 1 f" T Syringe Pump Figure 6.3: A schematic of the experimental configuration of the phantom and the ex-ternally situated syringe pump and driving motor. Arrangement of equipment is shown from the top. time, t, can be found through the equation, A(t) = 6 T * e - £* (6.4) Notice that this type of resulting function has both an increasing and decreasing portion and that the position of the peak activity is dependent upon the time the bolus was injected. It is possible to recreate such dynamic behaviour experimentally using the dynamic phantom as each of the dynamic containers has a small device that allows for the injection of a radiolabeled compound. This is composed of a thin rubber membrane that allows a syringe to pass through in order to inject a substance into the incoming water supply, yet seals itself so that no leakage occurs when a syringe is not present. Chapter 6. Phantom and Volunteer Experiments 110 6 . 2 R E C O N S T R U C T I O N A C C U R A C Y A S S E S S M E N T M E T H O D S In computer simulations, the true radiotracer distribution is known at each point in time for each voxel and so it is possible to evaluate the accuracy of image reconstruction by comparing the reconstructed activity distribution to the true object distribution for each voxel position at each time frame. The case of phantom experiments presents a different matter however as the true object activity distribution cannot be known at each point in time and so other methods must be used in order to assess the accuracy of the dSPECT image reconstruction. The most useful assessment method for experimental phantom data is a comparison of the true and the reconstructed dynamic parameters found from fitting time activity curves to a desired time dependent function. Since it is possible to extract time activity curves for a 3D volume of interest by using the dSPECT reconstruction method, these curves can' be fit to the appropriate time dependent function in order to determine dynamic function parameters. These dSPECT reconstructed dynamic parameters can then be compared to either the preset experimental parameters or to the true dynamic parameters obtained in the following manner. If an experiment is performed with a single dynamic volume in the absence of atten-uation, then it is possible to determine the true dynamic parameters for the object as a whole directly from projection measurements. As each projection contains only informa-tion from one volume and it is assumed that the radioactivity is distributed uniformly throughout this volume, then it is possible to obtain a time activity curve for the entire object by simply summing the total number of counts collected at each projection and plotting these counts versus time. Once the dSPECT method has been proven, it makes little sense to reconstruct only single dynamic objects. In order to create more complex experiments, projection data is Chapter 6. Phantom and Volunteer Experiments 111 first acquired for each dynamic region separately so that the true time activity curve pa-rameters can be obtained from the method previously described, but when the dSPECT image reconstruction is performed, multiple object projection data sets are created by simply summing the projections from individual objects into one final projection set that contains all dynamic objects. 6.3 D Y N A M I C P H A N T O M E X P E R I M E N T S Phantom experiments have been performed with different kinetic parameters and utilizing a variety of camera acquisition protocols. As well, experiments with and without phantom attenuation have been performed. A selection of these experiments will be presented in the following section along with the subsequent data analysis. Based upon the results of the previous chapter involving simulated data, the dynamic expectation maximization (dEM) reconstruction method was chosen for all experiments and no further testing of the CLS method has been performed. 6.3.1 NON-ATTENUATED EXPERIMENTS In the first set of experiments, the dynamic heart phantom was used without the presence of any additional thorax attenuation. Dynamic behaviour consisted of a mono exponen-tially decreasing activity over time. Acquisition parameters for these experiments are presented in Table 6.1. As mentioned, in the absence of attenuation, true dynamic parameters for each dy-namic object can be obtained provided that each object is imaged independently. After obtaining true time activity curves, they were then fit to either mono or dual exponential functions by minimizing the Y 2 value in a manner identical to that performed in the simulation studies. True dynamic parameters obtained in this way are summarized in Chapter 6. Phantom and Volunteer Experiments 112 Acquisition Head Rotation Projections Bins Time Type Geometry (per head) (per head) (minutes) 1) Single head n/a 180° 64 64 12 2) Triple head 120° 180° 64 64 12 3) Dual head 90° 90° 64 64 12 Table 6.1: Summary of experimental acquisition protocols used for experiments of the dynamic heart phantom with no attenuating material. Tables 6.2 - 6.4. Experiment 1 - Single head, 180° rotation Region A 0 A Tl/2 (total counts) (min - 1) (min) i 70770 0.12 5.84 ii 86687 0.20 3.54 iii 28236 0.28 2.45 iv 16637 0.12 5.58 Table 6.2: True projection time activity curve parameters for a mono exponential washout of activity acquired with a single head camera. R E S U L T S Summarized in Tables 6.6 - 6.7 and in Appendix B are results of dEM reconstructions of the dynamic heart phantom data previously described. In these summaries, true time activity curves are presented for each dynamic region, along with time activity curves obtained from 3D regions of interest obtained from the dEM reconstructions. It is evident from all of these reconstructions that the dEM reconstruction algorithm is capable of producing accurate reconstructions of dynamic objects provided that no appreciable object attenuation is present and that the objects are well differentiated in space. In all cases, the shape of the reconstructed TA curves agree very well to the true TA Chapter 6. Phantom and Volunteer Experiments 113 Experiment 2 - Triple head, 180° rotation Region A 0 ; A Ti/2 (total counts) (min - 1) (min) i 25634 0.27 2.61 ii 81013 0.18 3.82 iii 15234 . 0.11 6.41 iv 65581 0.11 6.47 Table 6.3: True projection time activity curve parameters for a mono exponential washout of activity acquired with a triple head camera. Experiment 3 - Dual head, 90° rotation Region A i Ai T V 2 , 1 A 2 A 2 Tl/2,2 Peak (counts) (min - 1) (min) (counts) (min - 1) (min) (min) i 25482 0.34 2.02 n/a n/a n/a n/a ii 22322 0.11 6.59 n/a n/a n/a n/a iii 21816 0.18 3.85 n/a n/a n/a n/a iv -7387 0.78 0.89 23063 0.08 8.2 1.2 v -8002 0.68 1.02 21439 0.04 18.3 3.0 Table 6.4: True projection time activity curve parameters for a dual exponential change of activity acquired with a dual head camera. curves, even in the case of a single head 180° acquisition. However, it should be noted that for this acquisition protocol, some streaking artifact is again present mainly appearing in the vicinity of the fastest dynamic region. As well, the dynamic fit parameters for the fastest dynamic region fit with the lowest accuracy thus indicating a difficulty of the dEM method to fit a fast changing activity with reasonable accuracy. It should also be noted at this point that in every case, the reconstructed time activity curve for the 3D dynamic region is consistently lower in amplitude than the true TA curve determined from projection measurements directly. This is mainly due to the way in which the TA curves are obtained rather than any inherent problems with the dEM Chapter 6. Phantom and Volunteer Experiments 114 Experiment 1 - Single head, 180° rotation Region A 0 A Tl/2 (total counts) (min - 1) (min) i 66339 0.12 5.78 ii 80069 0.20 3.55 iii 15661 0.21 3.25 iv 18867 0.16 4.43 Table 6.5: Reconstructed dynamic paramters of non-attenuated, single head dynamic phantom acquisition. Experiment 2 - Triple head, 180° rotation Region Ac A Tl/2 (total counts) (min - 1) (min) i 59903 0.11 6.51 ii 75807 0.18 3.84 iii 20213 0.24 2.84 iv 13800 0.11 6.32 Table 6.6: Reconstructed dynamic paramters of non-attenuated, triple head dynamic phantom acquisition. algorithm. Recall that the true TA curves are obtained from the total number of counts over the entire detector at each projection. Because of this, it is possible that photons that get scattered from the object at small angles still get counted by the detector, thus increasing the apparent number of counts from each region. In contrast, the 3D regions of interest used for the dSPECT reconstructions include only each dynamic region and so do not include any of the scattered photons that may actually appear to be coming from the background regions. Chapter 6. Phantom and Volunteer Experiments 115 Experiment 3 - Dual head, 90° rotation Region Ai Ai Ti/2,1 A 2 A2 Tl/2,2 Peak (counts) (min - 1) (min) (counts) (min - 1) (min) (min) i 22812 0.33 2.11 n/a n/a n/a n/a ii 22283 0.11 6.59 n/a n/a n/a n/a iii 20890 0.18 3.83 n/a n/a n/a n/a iv -7132 0.87 0.80 21739 0.09 8.06 1.1 ± 0.3 V -8575 0.54 1.29 21796 0.04 15.63 3.8 ± 1.2 Table 6.7: Reconstructed dynamic paramters of non-attenuated, dual head dynamic phantom acquisition. Figure 6.4: An attenuation map of the dynamic phantom thorax acquired with the Siemens Profile transmission system. Attenuation coefficients are in units of c m - 1 6 . 3 . 2 A T T E N U A T E D E X P E R I M E N T S Further to the non-attenuated experiments discussed previously, a selection of experi-ments were performed utilizing the entire thorax phantom with non-uniform attenuation (an image of the attenuation map is shown in Figure 6.4). Experimental acquisition protocols are summarized in Table 6.8. The experiments have been designed with three bottles undergoing a dynamic be-haviour consisting of a quick increase in activity followed by a slow washout. To do this, we have made use of the concept mentioned in equation (6.4) whereby a bolus injection of Chapter 6. Phantom and Volunteer Experiments 116 Detector Head Rotation Projections Bins Time Heads Geometry (per head) (per head) (minutes) 1) Dual head 90° 180° 64 64 25 2) Triple head 120° 180° 64 64 25 Table 6.8: Summary of the acquisition protocols for experiments involving attenuation and a combined increase/decrease. activity in injected into a constant flowing water stream in order to create the necessary washin / washout curve. The true washout time was preset by the rate of the flowing water and are summarized for each experiment in Table 6.9. It is impossible to know the exact time that the bolus was injected into the bottle as there is a short distance of tubing between the injection point and the dynamic bottle. Therefore, no accurate assessment can be performed with this data, and therefore only a qualitative description will be given for these experiments. In the fourth phantom bottle, a static activity was present in order to test the dSPECT method at differentiating static and dynamic activ-ities. As well, the entire thorax phantom was filled with a static background activity in order to further make the experiments more realistic. As an additional test of the dSPECT method, projection data acquired using the triple head camera was reconstructed with a system matrix that modelled the two di-mensional point response of the camera. It was thought that this improvement in the modelling geometry of the camera would result in both increased spatial resolution as well as increased reconstruction accuracy. R E S U L T S As seen from the images and data tables in Appendix B, the dEM algorithm was able to successfully reconstruct the three different dynamic regions present in the experiments Chapter 6. Phantom and Volunteer Experiments 117 Experiment 1 - Dual head, 180° rotation Region Ac A Ti/2 Injection (total counts) (min - 1) (min) Time (min) i 22032 0.173 4.0 1.0 ii 19321 0.116 6.0 2.5 iii 24009 0.058 12.0 5.0 iv 16899 0.0 oo 0.0 Experiment 2 - Triple head, 180° rotation Region A 0 A Ti/2 Injection (total counts) (min - 1) (min) Time (min) i 13022 0.173 4.0 1.0 ii 12877 0.116 6.0 2.5 iii 12082 0.058 12.0 5.0 iv 10990 0.0 oo 0.0 Table 6.9: Preset washout times for experiments involving both increase/decrease be-haviour, and nonuniform attenuation. Half-life times have been calculated based upon the preset water flow rates. as having different dynamic parameters. The dEM algorithm was also able to determine that the background activity in the phantom contained a static radioactive distribution. Time activity curves for the small static region however do not reconstruct with the same level of accuracy as the dynamic regions do. It can be seen in both the dual head and the triple acquisition reconstructions that the small bottle that contained a static activity does not reconstruct as static, but rather with a small dynamic component. In both acquisitions, the absolute change in activity for this static region was about 50%, thus indicating potential regions in the vicinity of dynamic regions. Time activity curves for the other dynamic regions appear to reconstuct well however with a good differentiation between TAC's for neighbouring regions. Chapter 6. Phantom and Volunteer Experiments 118 6.4 D Y N A M I C R E N A L E X P E R I M E N T S 6.4.1 M A G 3 E X P E R I M E N T S The first renal dSPECT study that was performed used the renal pharmaceutical 9 9 m T c labelled mercaptoacetylglycylglycylglycine (MAG3). This agent has excellent renal imag-ing properties owing to the fact that it undergoes both glomerular nitration, as well as tubule secretion. As a result of these two methods of elimination, the contrast between kidney and background is quite high. As well, the removal of MAG3 from the body is fast and shows excellent potential as a dynamic SPECT renal imaging agent. However, because of the two pathways of elimination, any measurements of renal filtration will re-flect a global filtration rate (ie, renal plasma flow), rather than a more specific filtration quantity such as glomerular filtration. For this study, a single healthy male subject was imaged using a Siemens ecam dual head SPECT camera with the heads in the 90° configuration. The camera heads were initially located at the 12 o'clock and 9 o'clock positions with the subject situated prone. The acquisition protocol consisted of acquiring 64 projections per head over a 90° circular, clockwise rotation. Data was acquired into 128x128 pixel projections (pixel size 4.795 mm) with an acquisition time of 10 seconds per projection starting at approximately 1 minute following a 370 MBq injection of MAG3. The central 96 transaxial camera bins were then truncated in order to reduce the reconstruction space size and a total of 70 transaxial slices were reconstructed thus creating a reconstructed object consisting of 64 time frames of size 96x96x70 voxels. Prior to the dSPECT acquisition described above, a conventional dynamic planar renogram was performed in order to compare with the subsequent dSPECT reconstruc-tion. For this study, a series of posterior projection measurements were acquired at 10 second intervals over a total time of 25 minutes. Data was acquired into 128x128 pixel Chapter 6. Phantom and Volunteer Experiments 119 projections by a Siemens single head Diacam camera starting immediately following 74 MBq injection of MAG3. M A G 3 R E S U L T S Since MAG3 is not solely a glomerular filtration agent, no quantitiative measurements of GFR can be obtained from this experiment, but rather a measurement only of effective renal plasma flow (ERPF). However at this point, in order to evaluate the accuracy of the dSPECT reconstruction technique, only a qualitative comparison between the dynamic planar TA curves and the dSPECT TA curves is made. To do this, the reconstructed dSPECT data was reprojected into planar projections and resized in order to create identical pixel sizes to the planar renogram projections. Additionally, the resultant planar dSPECT images were translated in order to align the kidneys into the same pixel positions as the renogram images. Following these modifications, identical regions of interest could be drawn on both the reprojected dSPECT images and the planar images and TA curves compared for these regions. Shown in the left half of Figure 6.5 is an image of the projection data at one time frame acquired using the planar acquisition method along with the TA curves for each kidney. It is seen that the kidneys undergo an initial increase in activity as a result of blood flow, and then the activity decreases slowly in accordance with redistribution in the body and elimination from the renal system. A reprojected dSPECT image is shown in the right half of Figure 6.5. Again it can be seen from the TA curves that an initial uptake of activity is present followed by a slow decline in activity similar to that of the planar TA curves. One of the problems with dynamic planar imaging is that any region of interest drawn on a 2D image includes the contribution from the desired organ, but also as a necessary Chapter 6. Phantom and Volunteer Experiments 120 1400 1200 1000 800 600 400 200 0 x 10 Left Kidney Right Kidney v i fl s \ i . . . .\~~^>v . . . . \ ; 10 12 0 2 4 6 8 10 12 Figure 6.5: Top) Images acquired using a dynamic planar imaging technique. Bottom) Reprojected dSPECT images into planar images. consequence, includes data from the background regions as well. Because of the fact that the dSPECT reconstruction method produces a series of transaxial images, it is possible to create regions of interest that include only the organ of interest and so avoid all the overlapping background activity present in planar projections. In Figure 6.6, a transaxial dSPECT image is shown along with TA curves drawn over each kidney. 6.5 C H A P T E R S U M M A R Y In this chapter an experimental validation of the dynamic expectation maximization routine has been presented. Through a series of experiments using a dynamic phantom that consisted of heterogeneous attenuation and different dynamic rates, it has been Chapter 6. Phantom and Volunteer Experiments 121 Figure 6.6: Images of transaxial slices reconstructed using the dEM reconstruction algo-rithm on data acquired with a single 90° acquisition protocol. The plateau seen in the TA curve for the right kidney may represent a second point of increasing activity due perhaps to a second pass of radiopharmaceutical within the kidney. As the reconstruction algorithm can only reconstruct TA curves with a single peak position at this time, the second peak appears as a plateau. found that the dEM reconstruction algorithm is able to determine regional dynamic rates accurately provided that projection data is acquired with a multi-headed SPECT camera and accurate attenuation maps are available. C H A P T E R 7 CONCLUSION This chapter provides a brief overview of this work and discusses future improvements to implement into the dSPECT reconstruction method. 7 . 1 G E N E R A L C O M M E N T S In this work, a new method of determining the rate of change of a three dimensional radiotracer distribution within an object has been presented. The method is unique in that it uses a standard nuclear medicine SPECT camera and acquires projection data with a standard SPECT acquisition (ie., single slow camera rotation). Additionally it is able to correct for patient specific attenuation as the method requires the inclusion of accurate attenuation maps in order to produce accurate reconstruction results. The primary advantage of the dSPECT reconstruction method that it utilizes a con-ventional SPECT camera with a single slow camera rotation in order to acquire projection data for input into the reconstruction routine. This is important in that the speed of data acquisition directly relates to the quality of data acquired and hence to the accuracy of image reconstruction (ie, fast rotation results in noisy, less accurate reconstructions). 122 Chapter 7. Conclusion 1 2 3 Additionally, because of the less stringent hardware requirements and acquisition proto-col, it allows for routine clinical use in almost any clinical setting where SPECT imaging is currently performed. The dSPECT reconstruction method has been implemented by incorporating a linear inequality constraint into standard reconstruction algorithms (ie., least squares, expecta-tion maximization). The choice of these two reconstruction algorithms was based upon the results obtained from static imaging. However, the dSPECT inequality constraints can be easily implemented into any reconstruction routine, thus allowing the possibil-ity to use the dSPECT linear inequalities with many new and improved reconstruction algorithms. 7 . 2 F U T U R E W O R K Currently the dEM reconstruction has been shown to produce reasonable results based upon the simulations and experiments performed up to this point. With further im-provements to the reconstruction process and modelling of photon transport within the body, even better results should be obtained. Some possible future work on the method includes the following. 7 . 2 . 1 F U L L Y F O U R DIMENSIONAL R E C O N S T R U C T I O N In the current implementation of the dSPECT image reconstruction technique, the sys-tem matrix modelled data acquisition using either a perfect parallel hole collimator ( 1 dimensional), or a transaxial camera response function ( 2 dimensional). While improve-ments were noticeable between the I D and 2 D reconstructions ( 2 D providing clearer, more well defined images) true 3 D collimator blurring correction will benefit the re-constructions even more. However, work still needs to be done in order that such a Chapter 7. Conclusion 1 2 4 reconstruction become feasible as a typical reconstruction may involve on the order of 1 2 8 4 object voxels. 7 . 2 . 2 P H O T O N S C A T T E R COMPENSATION For some time, the problem of photon scatter has been investigated in nuclear medicine imaging in order to improve image quality and quantitation [79]. Several new methods are emerging that hold the promise of correcting for this degradation in the images [80] [81]. While these methods are currently being developed for the case of static imaging, the extension into the dSPECT algorithm is quite possible. Thus, when such methods are proven to provide accurate estimations of photon scatter within a static distribution, they will be incorporated into the dSPECT procedure as well. 7 . 2 . 3 M Y O C A R D I A L G A T I N G One nuclear medicine procedure that was not discussed in this work is that of gated myocardial studies. In such studies, data acquisition is controlled by the cycle of the heart. That is, a complete heartbeat cycle is divided into several time gates and projection data is only collected at these times. Thus, in the same time that a static imaging procedure can be performed, a complete gated study can be performed in order to produce a three dimensional image of the heart at a variety of times in the heartbeat cycle. With such images, it is possible to view abnormal myocardial wall motion. Using the dSPECT reconstruction method, it is possible to extend normal gated myocardial studies in order to obtain information related to the dynamics of tracer motion in addition to the motion of the myocardium. A P P E N D I X A L I K E L I H O O D F U N C T I O N S A . l G A U S S I A N L I K E L I H O O D O B J E C T I V E F U N C T I O N In nuclear medicine imaging, projection measurements are acquired at various angles of a three dimensional radiotracer distribution. Because of the random nature of radioactive decay, in addition to stochastic processes such as photon attenuation, the measured number of counts detected by the camera y does not represent the true amount of activity present within the object. On average the number of detected counts has an equal probability to over represent the true activity as it does to under represent it. As we have no prior information about the true activity present, then it has to be assumed that the mean number of counts collected over time must be equal to the true amount of activity present. So, given that the true number of counts that should be collected in the jth camera bin at projection k is Yjk, then the probability of measuring counts can be determined from the Gaussian distribution, P(yjk\Yjk) =—e 2°> (A.l) 20jfc where Oj\. is equal to the standard deviation in the measured number of counts. Because each measurement is an independent quantity, we can write the probability of collecting 125 Appendix A. Likelihood Functions 126 the measured projection data over all camera bins, y, given that the true projection should be equal to Y, is equal to, Taking the natural logarithm of this equation and dropping terms independent of the quantity Yjk, we arrive at, We will refer to this equation as the log likelihood function, InC. Thus, by maximizing the log likelihood function (or minimizing the negative log likeliuhood) in the variable Yjk, we will arrive at the most probable mean number of counts, Yjk, given the actual measured data yjk- Notice that this is the typical L2 norm, or least squares distance between a measured quantity yjk, and a true quantity Yjk. Recall that the activity within each object voxel is related to the number of counts collected on the camera face through the equation, Yjk = E CijkXi (AA) i where C{jk is the projection operator matrix and %{ is the unknown activity distribution within the object. Thus, we can rewrite equation (A.3) in terms of the measured data, yjk, and the unknown source distribution, Xi as, _/„ C = £ (fcCWxi2-V*)A\ ( A . 5 ) j,k \ 2°~jk J which is the familiar least squares image reconstruction objective function described in Chapter 2. In the next section, instead of working with Gaussian distributed variables, Appendix A. Likelihood Functions 127 we will use Poisson distributed variables. We will see that this leads to a different log likelihood objective function as a result. A . 2 P O I S S O N L I K E L I H O O D O B J E C T I V E F U N C T I O N Given a Poisson distributed variable, the probability of measuring the quantity m given that the population mean is LI can be found through the Poisson probability distribution function, P(m\u) = (A.6) For the case of nuclear medicine imaging, we will consider that the number of counts collected in each camera bin obeys Poisson statistics as the underlying process of radioac-tive decay can be described as a Poisson process. Thus, the probability of collecting the number of counts y in the jth camera bin at projection k can be found through equation (A.6). Of course, the true mean number of counts that should be collected is not known a priori. Rather, this quantity will be determined through a data fitting procedure based upon the actual measured values, y^. Note that the true mean number of counts will be equal to the quantity, Yjk-J2CiJkx.i (A-7) i where represents the probability system matrix and x^ corresponds to the unknown activity in the ith object voxel. Thus by adjusting the value of X{ for each voxel, we will eventually find the a quantity that will give rise to true mean value Yjk when summed according to equation (A.7). Since all camera bins are considered as indepedent measurements, the total probability of collecting the entire data set, y, given that the true collected data should be equal to Appendix A. Likelihood Functions 128 the mean value, Y , can be found through the Poisson probability, p(y\Y) = n — r ( A - 8 ) j,k y*-We will now define this probality to be the likelihood function, C as it represents the likelihood of detecting the data y from the population with a mean of Y. Taking the natural logarithm of this function we obtain, In C = J2(-Yjk + yjklnYjk-ln(yjk\)) (A.9) By maximizing this function over the unknown parameters Yjk we would hence find the most probable mean that matches the collected data y. Note that maximizing this log likelihood function is identical to minimizing —InC. In the minimization, the last term in the likelihood function can further be neglected since that term is a constant in the unknown parameters. Thus, the final objective function to be minimized is equal to, -lnC = 52(Yjk-yjklnYjk) (A.10) j,k A . 3 T H E S T A T I C E X P E C T A T I O N M A X I M I Z A T I O N ( E M ) A L G O R I T H M Given a particular radioactive source distribution in an object, consider the number of photons emitted from the ith pixel to be Xj. We can imagine that the number of counts actually collected in the jth detector bin from an emission from this pixel is given by the quantity, (A.l l ) Appendix A. Likelihood Functions 129 where represents the probability that pixel % contributes to camera bin j at projection k. Note that this value is not a real measureable quantity but rather is related to the measured number of counts by, Vjk = Ylziok ( A-12) i Recall the negative log likelihood equation (A. 10) that relates the expected number of counts, (here C^Xi) to the true measured number of counts, (zijk). Thus, the negative log likelihood objective function can be written as, Fi = Y. (CijkXi - zijkln(CijkXi)) (A.13) j,k We wish to minimize this expression so differentiating and equating to zero gives us, { £ - E C w » - 5 i 4 M _ o (A.14) Rearranging and solving for Xi gives us, X i = | M i ^ L (A.15) 2^ 7,/c ^ ijk which, after substituting in for Zijk becomes, + _ Xj Ej.fc CjjkVjk a 1 f i \ xi - ^ 7^ v r V~ I A . I O J Here we have the iterative form of the E M algorithm. With it, an estimate of the true mean counts per pixel is created, and at every stage of the iterative process this value is updated to form a new starting point for the next iteration, xf. A P P E N D I X B D SPECT RECONSTRUCTIONS B . l A N A L Y T I C A L S I M U L A T I O N R E S U L T S In the following pages is a summary of dSPECT reconstructions performed using the annulus phantom with different camera acquisition parameters for the two reconstruction methods. B . l . l M O N O E X P O N E N T I A L W A S H O U T In this set of simulations, mono exponential washout was simulated using the equation, ln(2) t A(t) = AQe T V 2 (B.l) with functional parameters given by, 130 Appendix B. dSPECT Reconstructions 131 70 60 140 £30 I 20 10 — Section 1 — Section 2 —— Section 3 . • Section 4 \ \ : \ s ' \ i v ^ s " v ^ \ : N \ x 7 ^ T*. ~ - _ _ _ 10 Time (min) 15 20 Figure B . l : True annulus object shown at t = 5 min for mono exponential decay simu-lation. True Dynamic Parameters Section A T l /2 ,1 (per voxel) (min) 1 57.8 2.0 2 57.8 4.0 3 57.8 8.0 4 57.8 16.0 Appendix B. dSPECT Reconstructions 132 CLS Single Head - 180° Rotation w/noise 70 60 ^ 5 0 140 t f"30 1 20 10 0 • Section 1 Section 2 — Section 3 • Section 4 \ V \ \ " ' 'S i ' V ' ' ' \ \ > \ \ \ \ ? v . - - _ 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section A Tl/2,1 (per voxel) (min) 1 17.1 4.8 0.73 2 48.6 4.3 0.25 3 58.5 6.9 0.19 4 65.4 10.3 0.22 Overall 0.41 Reconstructed Parameters - w/Noise Section A Tl/2,1 £R (per voxel) (min) 1 17.8 4.6 0.76 2 48.6 4.3 0.34 3 56.1 7.1 0.29 4 64.4 10.5 0.26 Overall 0.46 Appendix B. dSPECT Reconstructions 133 dEM Single Head - 180° Rotation w/noise 70 • 50 10 Section 1 Section 2 — Section 3 Section 4 \ v •• \ \ \ \ \ \ s N 15 Reconstructed Parameters - Noiseless Section Ai Ti/2,1 (per voxel) (min) 1 9.7 4.9 0.81 2 48.5 4.6 0.27 3 64.9 6.7 0.20 4 63.6 11.8 0.24 Overall 0.46 Reconstructed Parameters w /Noise Section Ai Tl/2,1 £R (per voxel) (min) 1 9.5 4.9 0.82 2 48.5 4.6 0.28 3 65.3 6.7 0.21 4 63.6 11.8 0.25 Overall 0.46 Appendix B. dSPECT Reconstructions 134 CLS Dual Head - 180° Rotation w/noise 70 60 "5 f 40 c . f"30 1 20 10 0 Section 1 Section 2 Section 3 • • Section 4 \ \ ^ \ \ > V \ V. N . \ \ V \ N N. v N ~~ - ^ ^ 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Tl/2,1 SR (per voxel) (min) 1 44.9 2.4 0.33 2 54.3 4.0 0.20 3 56.0 8.0 0.16 4 59.5 13.5 0.15 Overall 0.22 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 £R (per voxel) (min) 1 44.6 2.4 0.44 2 54.3 4.0 0.29 3 55.7 8.0 0.22 4 58.7 13.7 0.18 Overall 0.30 Appendix B. dSPECT Reconstructions 135 dEM Dual Head - 180° Rotation w/noise 70 60 — 5 0 f 40 g | " 3 0 I 20 10 0 \ \ A . V . \ . . Section 1 Section 2 —— Section 3 • Section 4 \ \ V V \ Y A - . \ \ ' S \ ^ ?v 5 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 £R (per voxel) (min) 1 44.1 2.4 0.37 2 54.8 4.0 0.24 3 55.3 8.1 0.19 4 58.3 14.7 0.19 Overall 0.26 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 £R (per voxel) (min) 1 43.1 2.5 0.38 2 55.8 4.1 0.25 3 53.7 8.2 0.22 4 61.7 14.8 0.20 Overall 0.27 Appendix B. dSPECT Reconstructions 136 CLS Triple Head - 180° Rotation w/noise 70 60 50 CD | 4 0 t | " 3 0 2 20 10 0 Section 1 Section 2 — Section 3 . Section 4 \ s \ "*' \ N N V \ \ \ \ \ \ \ N ~* — — — ^ 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Tl/2,1 £R (per voxel) (min) 1 52.7 2.2 0.27 2 55.7 4.1 0.17 3 56.7 7.8 0.16 4 58.6 14.1 0.15 Overall 0.19 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 eR (per voxel) (min) 1 52.5 2.2 0.38 2 55.5 4.1 0.29 3 56.6 7.9 0.21 4 58.2 14.4 0.19 Overall 0.28 Appendix B. dSPECT Reconstructions 137 dEM Triple Head - 180° Rotation w/noise 70 r f 40 I" 30 I 20 Section 1 Section 2 — Section 3 • Section 4 t \ \ • - - • \ x \ \ v \ \ \ > V \ \ \ \ N s v. 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 £R (per voxel) (min) 1 50.6 2.2 0.32 2 55.9 4.0 0.20 3 55.8 8.0 0.19 4 56.9 15.5 0.19 Overall 0.23 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 £R (per voxel) (min) 1 50.6 2.2 0.33 2 56.0 4.0 0.21 3 55.6 8.0 0.20 4 56.9 15.5 0.20 Overall 0.24 Appendix B. dSPECT Reconstructions 138 B . l . 2 D U A L E X P O N E N T I A L W A S H O U T Dual exponential simulations were modelled by the equation, / ln2 t ln.2 t \ A(t) = A0 (e T l / 2 , i + E ^ (B.2) and with true functional parameters given by, 70 60 — 5 0 140 I .|" 30 i 20 10 0 Section 1 Section 2 — Section 3 . Section 4 V s . • 10 Time (min) 15 20 Figure B .2 : True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation. True Dynamic Parameters Section A Tl/2,1 Tl/2,2 (per voxel) (min) (min) 1 28.9 2.0 20.0 2 28.9 4.0 20.0 3 28.9 8.0 20.0 4 28.9 16.0 20.0 Appendix B. dSPECT Reconstructions 139 CLS Single Head - 180° Rotation w/noise 70 60 _ 5 0 | 4 0 I .§"30 S 20 10 0 Section 1 Section 2 — Section 3 . - Section 4 \ \; \ - N -• • - > S i - • -s \ . V 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 17.2 5.8 51.2 0.23 2 27.9 3.7 17.0 0.19 3 27.8 11.2 11.2 0.22 4 31.6 7.0 18.1 0.23 Overall 0.20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 17.0 6.2 48.3 0.29 2 25.6 3.8 17.1 0.23 3 26.9 11.6 11.6 0.23 4 31.4 7.4 17.9 0.24 Overall 0.25 Appendix B. dSPECT Reconstructions 140 dEM Single Head - 180° Rotation w/noise 70 60 ^ 5 0 © f 40 c £ " 3 0 o < 20 10 0 -V : Section 1 Section 2 — Section 3 -• Section 4 x \ *"* . \ \ v X . \ ./SL.. . \ ^ "**.. X V v . \ . : \ " X . . , / .. ,"s . . . 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 19.3 5.2 36.8 0.21 2 27.6 3.9 19.5 0.21 3 28.0 6.7 24.3 0.17 4 31.1 6.4 37.4 0.20 Overall 0.20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 19.2 5.1 37.5 0.23 2 27.5 3.9 25.7 0.22 3 28.0 6.7 24.7 0.18 4 30.9 6.5 36.7 0.21 Overall 0.21 Appendix B. dSPECT Reconstructions 141 CLS Dual Head - 180° Rotation w/noise Section 1 Section 2 — Section 3 • Section 4 0 5 10 15 20 Time (min) Reconstructed Parameters - Noiseless Section Ai T l /2,1 Tl/2,2 £R (per voxel) (min) (min) 1 25.4 2.8 21.8 0.16 2 27.4 4.4 17.6 0.16 3 27.6 12.2 12.2 0.16 4 28.8 15.6 15.6 0.15 Overall 0.16 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 25.5 2.7 21.9 0.21 2 27.1 4.6 18.0 0.22 3 27.5 12.3 12.3 0.21 4 28.7 15.7 15.7 0.20 Overall 0.21 Appendix B. dSPECT Reconstructions dEM Dual Head - 180° Rotation w/noise 0 5 10 Time (min) Reconstructed Parameters - Noiseless Section Ai T l /2 ,1 Tl/2,2 £R (per voxel) (min) (min) 1 25.5 2.7 22.8 0.19 2 27.8 3.9 20.8 0.19 3 27.7 8.0 20.5 0.19 4 28.6 11.2 28.0 0.18 Overall 0.19 Reconstructed Parameters - w/Noise Section Ai T l /2 ,1 Tl/2,2 £R (per voxel) (min) (min) 1 25.3 2.7 22.8 0.20 2 27.9 3.8 20.8 0.20 3 27.6 8.3 20.0 0.19 4 28.6 11.0 28.7 0.19 Overall 0.20 Appendix B. dSPECT Reconstructions 143 CLS Triple Head - 180° Rotation w/noise 70 60 140 2 20 10 0 \vN**--, \ s *"**• Section 1 Section 2 — Section 3 -• Section 4 V • * V ' 1 *^SL ' '" ' **** \ ^ x ™—~ 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 27.6 2.27 19.92 0.15 2 27.9 4.52 18.09 0.15 3 28.1 12.16 12.16 0.16 4 29.0 16.15 16.15 0.15 Overall 0.15 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 27.3 2.3 20.2 0.21 2 27.7 4.5 18.3 0.21 3 27.9 12.2 12.2 0.20 4 28.9 16.0 16.0 0.19 Overall 0.20 Appendix B. dSPECT Reconstructions 144 dEM Triple Head - 180° Rotation w/noise 70 r Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 27.3 2.21 20.68 0.19 2 28.0 3.95 20.37 0.19 3 28.0 7.62 21.25 0.19 4 28.2 14.59 21.30 0.19 Overall 0.19 Reconstructed Parameters - w/Noise Section Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) 1 27.5 2.2 20.7 0.20 2 27.8 4.0 20.4 0.19 3 28.1 7.5 21.5 0.19 4 28.1 15.8 19.7 0.19 Overall 0.20 Appendix B. dSPECT Reconstructions 145 B . l . 3 D U A L E X P O N E N T I A L WITH S L O W W A S H I N A N D W A S H O U T Simulations involving both washin and washout were simulated with the dual exponential function, / tn2 t ln2 t \ A(t) = A0 le T i / 2 , i - E J (B.3) with true dynamic parameters consisting of, — Section 1 Section 2 — Section 3 • Section 4 0 5 10 15 20 Time (min) Figure B .3 : True image at t = 5 min along with time activity curves for each dynamic region for the slow dual exponential simulation. True Dynamic Parameters Section A, T l /2 ,1 Tl/2,2 Peak Position (per voxel) (min) (min) (frame no.) 1 28.9 0.5 20.0 9.0 2 28.9 1.0 20.0 16.0 3 28.9 2.0 20.0 24.0 4 28.9 4.0 20.0 28.0 Appendix B. dSPECT Reconstructions 146 CLS Single Head - 180° Rotation w/noise 30 r 25 •20 1« > <c 10 Section 1 I Section 2 — Section 3 • Section 4 |" j i •• f 5 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 23.7 1.0 30.9 0.23 26.6 ± 18.2 2 26.9 1.0 19.5 0.19 20.4 ± 7.8 3 25.4 1.4 22.8 0.19 26.9 ± 5.8 4 12.4 1.2 -294.6 0.23 33.2 ± 8.0 Overall 0.21 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 23.07 1.1 30.8 0.38 24.9 ± 21.6 2 27.06 1.0 17.5 0.37 24.8 ± 11.7 3 25.13 1.2 23.6 0.43 28.4 ± 7.4 4 10.98 1.0 -57.9 0.51 36.2 ± 7.2 Overall 0.43 Appendix B. dSPECT Reconstructions 147 dEM Single Head - 180° Rotation w/noise 30 r 25 f 20 x I S 15 £ 10 — Section 1 Section 2 — Section 3 • Section 4 / 1 / if If 1/ I-5 10 15 Time (min) 20 Reconstructed Parameters - Noiseless Section Ai T l / 2 , 1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 25.0 1.2 30.6 0.33 28.7 ± 23.0 2 29.1 1.2 18.0 0.27 21.4 ± 8.6 3 28.8 1.5 18.0 0.28 25.3 ± 5.9 4 11.9 1.1 177.5 0.38 37.8 ± 12.5 Overall 0.32 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 24.5 1.2 32.9 0.36 29.6 ± 23.7 2 29.5 1.2 17.6 0.29 21.4 ± 8.5 3 29.4 1.5 17.5 0.32 26.4 ± 6.0 4 11.7 1.1 199.0 0.42 38.4 ± 13.1 Overall 0.35 Appendix B. dSPECT Reconstructions CLS Dual Head - 180° Rotation w/noise Section 1 Section 2 — Section 3 • Section 4 0 5 10 15 20 Time (min) Reconstructed Parameters - Noiseless Section A, Tl/2,1 Tl/2,2 eR Peak Position (per voxel) (min) (min) (frame no.) 1 28.2 0.6 20.0 0.29 18.2 ± 11.1 2 26.5 0.9 22.0 0.33 21.7 ± 8.9 3 28.4 2.0 19.2 0.35 26.6 ± 5.8 4 16.5 2.5 78.0 0.42 34.9 ± 8.3 Overall 0.35 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 28.2 0.6 20.0 0.29 18.2 ± 11.1 2 26.5 0.9 22.0 0.33 21.7 ± 8.9 3 28.4 2.0 19.2 0.35 26.6 ± 5.8 4 16.5 2.5 78.0 0.42 34.9 ± 8.3 Overall 0.35 Appendix B. dSPECT Reconstructions dEM Dual Head - 180° Rotation w/noise 30 r 25 ^ 2 0 X a I 115 > < 10 Section 1 Section 2 — Section 3 • Section 4 / ' x ^ / ' / ' / • 1/ / : /' / .*•'"'* h.I /' / •'' /' / •' 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 29.9 0.7 17.7 0.25 16.2 ± 9.4 2 27.7 1.0 20.6 0.25 20.5 ± 8.2 3 33.8 2.5 15.3 0.23 28.0 ± 6.2 4 18.1 2.7 49.3 0.25 40.1 ± 6.4 Overall 0.25 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 29.9 0.7 17.8 0.27 16.7 ± 9.9 2 27.7 1.0 20.4 0.28 20.8 ± 8.5 3 33.6 2.5 15.4 0.28 29.1 ± 7.5 4 17.8 2.6 54.0 0.29 39.7 ± 6.4 Overall 0.28 Appendix B. dSPECT Reconstructions 150 CLS Triple Head - 180° Rotation w/noise 30r 25 f 20 x f I 15 > % 10 Section 1 Section 2 — Section 3 • Section 4 / / / ' / I / / ' / I' / i i . » v ! li / / !i / .•' ' / ' i 10 Time (min) 15 20 Reconstructed Parameters - Noiseless Section 4 Tl/2,1 Tl/2,2 eR Peak Position (per voxel) (min) (min) (frame no.) 1 28.2 0.5 20.2 0.16 15.0 ± 8.0 2 28.2 1.0 20.1 0.16 19.3 ± 6.4 3 29.2 2.1 19.0 0.16 26.8 ± 5.0 4 35.0 4.5 16.0 0.17 34.5 ± 7.1 Overall 0.17 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 28.4 0.5 19.6 0.28 17.6 ± 10.8 2 28.0 1.0 20.6 0.30 21.4 ± 9.2 3 29.5 2.1 18.5 0.35 27.9 ± 7.8 4 29.8 4.2 18.9 0.41 34.8 ± 8.5 Overall 0.34 Appendix B. dSPECT Reconstructions dEM Triple Head - 180° Rotation w/noise 30r 10 15 Time (min) Reconstructed Parameters - Noiseless Section Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 29.1 0.6 18.7 0.23 14.8 ± 7.6 2 29.6 1.1 18.2 0.23 19.2 ± 6.0 3 32.0 2.4 16.2 0.22 28.4 ± 6.0 4 32.3 4.4 17.4 0.24 39.3 ± 5.5 Overall 0.23 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 28.9 0.6 19.0 0.26 15.9 ± 8.8 2 29.6 1.1 18.1 0.25 26.4 ± 7.5 3 31.8 2.4 16.7 0.26 29.6 ± 7.6 4 34.0 4.5 16.6 0.28 39.2 ± 8.8 Overall 0.26 Appendix B. dSPECT Reconstructions 1 5 2 B . l . 4 D U A L E X P O N E N T I A L WITH F A S T W A S H I N A N D W A S H O U T Simulations involving fast washin and washout were modelled with the same dual expo-nential function in time as previously, ln2 t A(t) = A0 [e T ^ but with functional parameters given as, ( B . 4 ) •5.12 I.10 > 8 1 6 4 2 0 Section 1 — Section 2 — Section 3 • Section 4 . /' A -: \ ... .v.«.^ ,... 1' / \ \ \ Ir / . V lij • \ \ \ . . 1/. \ (/ • X «. _ 5 10 15 Time (min) 20 Figure B . 4 : True image at t = 5 min along with time activity curves for each dynamic region for the rapid dual exponential simulation. True Dynamic Parameters Section At Tl/2,1 Tl/2,2 Peak Position (per voxel) (min) (min) (frame no.) 1 2 8 . 9 0 . 5 2 . 0 5 . 0 2 2 8 . 9 1 . 0 4 . 0 9 . 0 3 2 8 . 9 2 . 0 8 . 0 1 5 . 0 4 2 8 . 9 4 . 0 1 6 . 0 3 5 . 0 Appendix B. dSPECT Reconstructions 153 CLS Single Head - 180° Rotation w/noise — Section 1 — Section 2 — Section 3 • Section 4 • t s» 0 5 10 15 20 Time (min) Reconstructed Parameters - Noiseless Section Ai T l /2 ,1 Tl/2,2 eR Peak Position (per voxel) (min) (min) (frame no.) 1 8.3 0.1 4.5 1.46 1.0 ± 4.8 2 31.0 0.9 3.6 0.53 9.6 ± 4.7 3 23.3 1.4 9.1 0.30 14.7 ± 7.5 4 29.0 3.6 14.5 0.31 28.8 ± 11.9 Overall 0.81 Reconstructed Parameters - w/N oise Section Ai Ti/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 2.0 0.6 9.9 1.15 1.3 ± 7.8 2 25.6 1.1 4.7 0.59 8.9 ± 3.4 3 22.2 0.6 8.6 0.53 8.9 ± 9.7 4 14.6 1.5 37.9 0.41 18.1 ± 18.5 Overall 0.73 Appendix B. dSPECT Reconstructions dEM Single Head - 180° Rotation w/noise 1 201 1 1 1 — Section 1 Time (min) Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 3.6 0.41 6.24 0.86 5.8 ± 10.8 2 22.4 1.03 5.04 0.34 10.7 ± 3.2 3 24.6 1.03 7.94 0.29 11.9 ± 6.8 4 11.9 1.03 85.51 0.40 30.7 ± 24.4 Overall 0.53 Reconstructed Parameters - w/N oise Section Ai Ti/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 3.7 0.40 6.07 0.88 7.8 ± 14.3 2 21.8 1.00 5.06 0.39 10.5 ± 3.4 3 25.1 1.11 7.82 0.34 12.2 ± 6.8 4 11.7 0.95 104.5 0.45 31.3 ± 24.4 Overall 0.56 Appendix B. dSPECT Reconstructions CLS Dual Head - 180° Rotation w/noise Section 1 Section 2 — Section 3 • Section 4 0 5 10 15 20 Time (min) Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 19.0 0.5 2.6 1.13 3.7 ± 4.4 2 27.6 0.9 3.9 0.47 8.8 ± 4.6 3 27.1 1.9 8.4 0.22 18.2 ± 4.1 4 17.3 2.5 31.4 0.21 33.4 ± 7.3 Overall 0.61 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 20.2 0.5 2.4 1.17 3.8 ± 4.9 2 28.8 0.9 3.8 0.60 9.3 ± 4.8 3 25.4 1.8 8.7 0.42 19.4 ± 5.5 4 17.4 2.6 31.5 0.41 30.4 ± 8.5 Overall 0.72 Appendix B. dSPECT Reconstructions dEM Dual Head - 180° Rotation w/noise — Section 1 18 — Section 2 1 — Section 3 16 I ' Section 4 Time (min) Reconstructed Parameters - Noiseless Section Ai Ti/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 22.8 0.57 2.38 0.42 5.3 ± 2.1 2 26.6 0.92 4.02 0.28 8.7 ± 2.9 3 34.6 2.31 7.07 0.26 17.4 ± 4.3 4 17.7 2.61 30.19 0.42 30.3 ± 13.4 Overall 0.31 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 22.4 0.56 2.41 0.49 5.6 ± 2.2 2 26.6 0.92 4.01 0.34 9.0 ± 3.4 3 34.5 2.30 7.05 0.29 17.2 ± 4.3 4 17.8 2.60 29.91 0.32 29.3 ± 14.2 Overall 0.37 Appendix B. dSPECT Reconstructions 157 CLS Triple Head - 180° Rotation w/noise Section 1 18 — Section 2 1 — Section 3 16 • Section 4 0 5 10 15 20 Time (min) Reconstructed Parameters - Noiseless Section Ai Tl/2,1 Tl/2,2 eR Peak Position (per voxel) (min) (min) (frame no.) 1 8.3 0.1 4.6 1.46 1.0 ± 4.8 2 31.1 0.9 3.6 0.53 9.5 ± 4.7 3 23.4 1.4 9.1 0.30 14.7 ± 7.5 4 29.0 3.7 14.5 0.31 28.8 ± 11.9 Overall 0.81 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 eR Peak Position (per voxel) (min) (min) (frame no.) 1 9.4 0.2 4.0 1.57 1.2 ± 5.2 2 32.8 0.9 3.6 0.73 9.9 ± 4.7 3 29.1 1.5 8.9 0.45 16.4 ± 6.9 4 29.4 3.8 14.4 0.46 28.8 ± 11.5 Overall 0.92 Appendix B. dSPECT Reconstructions 158 dEM Triple Head - 180° Rotation w/noise 20 1B 16 U •5.12 110 2-Section 1 Section 2 Section 3 Section 4 . / A / N /' Y s. e ^ r I V* • \ Ii / .' V -1 Time 0 (min) 15 20 Reconstructed Parameters - Noiseless Section A Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 22.1 0.48 2.39 0.43 5.5 ± 2.2 2 27.1 0.94 4.05 0.26 8.8 ± 2.9 3 33.8 2.31 7.14 0.24 16.6 ± 4.7 4 20.1 2.95 23.77 0.28 26.6 ± 17.2 Overall 0.31 Reconstructed Parameters - w/N oise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) 1 22.1 0.47 2.39 0.50 5.6 ± 2.2 2 27.2 0.95 4.06 0.32 9.0 ± 2.9 3 33.7 2.28 7.12 0.28 17.2 ± 4.6 4 20.8 3.05 22.16 0.32 26.1 ± 17.2 Overall 0.37 Appendix B. dSPECT Reconstructions 159 B . 2 D I S C R E T I Z E D S I M U L A T I O N R E S U L T S This section presents results of dSPECT reconstructions of the dynamic thorax phantom acquired with the same acquisition protocols as previously. B.2.1 M O N O E X P O N E N T I A L W A S H O U T Mono exponential decay was again modelled by the equation (B.l), with decay parameters for the different thorax regions given by, Heart 1 - - Heart 2 . — Thorax \ N \ v V \ V \ \ ' x \ \ s \ • X N. ~- - „ 0 5 10 15 20 Time (min) Figure B.5: True image at t = 5 min along with time activity curves for each dynamic region for the mono exponential decay simulation. True Dynamic Parameters Section A* Tl/2,1 (per voxel) (min) heart 1 101.3 4.0 heart 2 101.3 8.0 lungs 0.3 oo thorax 6.1 oo Appendix B. dSPECT Reconstructions 160 CLS Single Head - 180° Rotation w/noise 9 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section A Tl/2,1 £R (per voxel) (min) heart 1 35.1 7.3 0.58 heart 2 80.9 8.3 0.27 lungs 1.7 79.9 6.07 thorax 9.2 16.9 0.79 Overall 2.78 Appendix B. dSPECT Reconstructions 161 dEM Single Head - 180° Rotation w/noise 110r 100r 90 80 1 70 | 60 & 5 0 § 40 30 20 10 0 Heart 1 - - Heart 2 . Thorax • N. 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Ti/2,1 £R (per voxel) (min) heart 1 26.7 11.5 0.62 heart 2 64.3 13.5 0.29 lungs 1.7 45.8 4.76 thorax 9.3 20.7 0.75 Overall 2.16 Appendix B. dSPECT Reconstructions 162 CLS Dual Head - 180° Rotation w/noise 110 100 90 80 lixel; 70 $ 60 5 f 50 > 40 < 30 20 10 0 — Heart 1 — - Heart 2 . — Thorax N S. V . \ . \ s. S \ < >X N . . . V. \ \ , iTy, 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 £R (per voxel) (min) heart 1 65.5 5.5 0.41 heart 2 97.2 7.8 0.22 lungs 1.2 58.8 4.48 thorax 7.2 39.6 0.55 Overall 2.05 Appendix B. dSPECT Reconstructions 163 dEM Dual Head - 180° Rotation w/noise 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai Tl /2 ,1 £R (per voxel) (min) heart 1 55.1 6.2 0.38 heart 2 90.0 9.0 0.18 lungs 1.1 61.1 2.98 thorax 7.5 37.1 0.46 Overall 1.38 Appendix B. dSPECT Reconstructions 164 CLS Triple Head - 180° Rotation w/noise i f f Heart 1 - - Heart 2 . —— Thorax \ \ . \ X " ^ X x \ ^ \ N. N , • •' > s 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai Tl/2,1 £R (per voxel) (min) heart 1 75.2 4.9 0.34 heart 2 94.3 8.1 0.23 lungs 0.9 31.4 3.92 thorax 6.9 66.7 0.46 Overall 1.78 Appendix B. dSPECT Reconstructions 165 dEM Triple Head - 180° Rotation — Heart 1 — - Heart 2 — Thorax \ \ \ x • \ N ' \ x \ ^ \ s - ••" - - i - -• 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai T l / 2 , 1 £R (per voxel) (min) heart 1 72.9 5.3 0.25 heart 2 90.9 8.9 0.17 lungs 1.1 46.1 2.47 thorax 7.2 45.7 0.37 Overall 1.14 Appendix B. dSPECT Reconstructions 166 B .2 .2 D U A L E X P O N E N T I A L W A S H O U T Dual exponential simulations modelled the equation (B .2), with functional parameters of, Figure B .6: True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential decay simulation. True Dynamic Parameters Section Ai Tl/2,1 Tl/2,1 (per voxel) (min) (min) heart 1 50.6 4.0 12.0 heart 2 50.6 8.0 20.0 lungs 0.2 oo oo thorax 2.2 oo oo Appendix B. dSPECT Reconstructions 167 CLS Single Head - 180° Rotation w/noise 110 100 90 80 ixel) 70 1 60 50 > 40 < 30 20 10 0 — Heart 1 - - Heart 2 . —— Thorax ; N N. s \ ' — • • \ — VI 10 Time (min 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,1 £R (per voxel) (min) (min) heart 1 29.6 9.4 9.4 0.37 heart 2 45.7 11.0 11.0 0.25 lungs 0.9 62.5 62.5 6.75 thorax 4.2 23.6 23.6 0.74 Overall 3.06 Appendix B. dSPECT Reconstructions 168 dEM Single Head - 180° Rotation w/noise 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,1 £R (per voxel) (min) (min) heart 1 24.7 6.3 75.5 0.37 heart 2 38.2 9.0 58.7 0.22 lungs 0.8 48.3 48.3 4.65 thorax 4.4 78.6 78.6 0.69 Overall 2.10 Appendix B. dSPECT Reconstructions 169 CLS Dual Head - 180° Rotation w/noise 110 100 90 80 xel) 70 CL 60 I 50 > 40 < 30 20 10 0 Heart 1 - - Heart 2 -— Thorax x s x.....: s ^x V X . X X x. S ~~ N \ i : : ; -x 5 10 15 Time (min) 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,1 £R (per voxel) (min) (min) heart 1 38.9 8.5 8.5 0.28 heart 2 50.3 10.7 10.7 0.19 lungs 0.8 3.9 -35.5 4.95 thorax 3.5 58.4 58.4 0.54 Overall 2.25 Appendix B. dSPECT Reconstructions 170 dEM Dual Head - 180° Rotation w/noise 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) heart 1 36.2 5.0 17.5 0.25 heart 2 46.0 7.6 27.5 0.17 lungs 0.7 79.1 79.1 3.65 thorax 3.7 45.1 45.1 0.49 Overall 1.67 Appendix B. dSPECT Reconstructions 171 CLS Triple Head - 180° Rotation w/noise 110 100 90 80 (ixel) 70 60 50 > Tj 40 < 30 20 10 0 Heart 1 Heart 2 Thorax 5 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R (per voxel) (min) (min) heart 1 42.5 8.2 8.2 0.29 heart 2 48.7 11.8 11.8 0.21 lungs 0.5 5.4 -80.8 4.17 thorax 3.5 14.6 -44.6 0.49 Overall 1.90 Appendix B. dSPECT Reconstructions 172 dEM Triple Head - 180° Rotation w/noise 110 100 90 80 iixel) 70 i 60 ty (cr ty (cr 50 > 40 < 30 20 10 0 — Heart 1 — - Heart 2 — Thorax \ \ \ x \ x X x >. ; x 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,1 £R (per voxel) (min) (min) heart 1 43.0 3.9 17.2 0.21 heart 2 47.4 6.7 29.6 0.15 lungs 0.5 62.8 62.8 2.61 thorax 3.5 55.7 55.7 0.39 Overall 1.20 Appendix B. dSPECT Reconstructions 173 B . 2 . 3 D U A L E X P O N E N T I A L WITH W A S H I N A N D W A S H O U T The final simulation modelled a combined washin and washout and followed the equation, (B.3), with functional parameters specified as, — Heart 1 — - Heart 2 — Thorax Time (min) Figure B.7: True image at t = 5 min along with time activity curves for each dynamic region for the dual exponential increase and decrease simulation. True Dynamic Parameters Section Ai Tl/2,1 Tl/2,2 Peak Position (per voxel) (min) (min) (frame no.) heart 1 50.7 1.0 8.0 11.0 heart 2 50.7 2.0 10.0 19.0 lungs 0.2 n/a n/a 0 thorax 6.1 n/a n/a 0 Appendix B. dSPECT Reconstructions v CLS Single Head - 180° Rotation w/noise 40 r m « » 35 30 "55 $ 25 CL 1.20 I 15 < 10 5 0 — Heart 1 — - Heart 2 — Thorax Is \ ~ ~ 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) heart 1 13.5 0.1 32.7 0.58 11.9 ± 6.9 heart 2 14.3 0.1 111.5 0.51 19.5 ± 8.3 lungs 3.0 0 0.6 13.81 2.5 ± 9.7 thorax 6.4 0.1 86.3 1.00 21.9 ± 27.3 Overall 6.22 Appendix B. dSPECT Reconstructions 175 dEM Single Head - 180° Rotation w/noise 5.20 — Heart 1 — - Heart 2 — Thorax ' -0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Ai T l /2 ,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) heart 1 22.3 0.3 26.2 0.42 16.0 ± 8.0 heart 2 24.5 0.5 28.9 0.35 15.6 ± 11.3 lungs 2.1 0.0 oo 7.30 0.0 ± 0.0 thorax 7.1 0.1 65.9 0.65 9.3 ± 16.3 Overall 3.32 Appendix B. dSPECT Reconstructions 176 CLS Dual Head - 180° Rotation w/noise 1 * • * 0 5 10 15 20 Time (min) Reconstructed Dynamic Parameters - w/Noise Section Ai T l /2 ,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) heart 1 29.7 0.8 14.3 0.46 17.6 ± 8.2 heart 2 26.8 0.9 21.0 0.47 22.6 ± 7.9 lungs 1.3 0.0 743.3 6.35 -0.29 ± 27.4 thorax 6.4 0.1 144.6 0.93 23.3 ± 29.6 Overall 2.92 Appendix B. dSPECT Reconstructions 177 dEM Dual Head - 180° Rotation w/noise 40 35 30 © .5 25 a. 120 * 1 15 < 10 5 0 Heart 1 - - Heart 2 — Thorax " fl II \ / \ \ \ V li li .. . // i It 11 Time 15 20 'min) Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) heart 1 35.9 1.1 10.2 0.37 16.6 ± 7.5 heart 2 43.5 1.7 11.4 0.36 20.4 ± 6.3 lungs 1.0 0.0 oo 3.27 0.2 ± 2.8 thorax 6.2 0.0 224.8 0.42 0.8 ± 4.6 Overall 1.50 Appendix B. dSPECT Reconstructions 178 CLS Triple Head - 180° Rotation w/noise 40 35 30 S 25 f 120 | 15 10 5 0 — Heart 1 — - Heart 2 — Thorax ' • / ' ... y j 1 S \ " ^ 1 " 10 Time (min) 15 20 Reconstructed Parameters - w/Noise Section Ai Tl/2,1 Tl/2,2 £R Peak Position (per voxel) (min) (min) (frame no.) heart 1 35.0 0.5 10.2 0.45 14.8 ± 5.9 heart 2 29.1 0.7 20.1 0.54 23.2 ± 8.8 lungs 0.8 0.0 164.3 4.72 -0.7 ± 1.0 thorax 6.2 0.1 780.9 0.87 15.2 ± 24.1 Overall 2.21 Appendix B. dSPECT Reconstructions dEM Triple Head - 180° Rotation w/noise ft* 0 5 10 15 20 Time (min) Reconstructed Parameters - w/Noise Section Tl/2,1 T l / 2 ,2 eR Peak Position (per voxel) (min) (min) (frame no.) heart 1 46.4 1.2 8.8 0.32 14.4 ± 5.3 heart 2 46.7 1.9 10.6 0.35 17.9 ± 8.2 lungs 0.8 0.0 oo 2.40 0.0 ± 0.0 thorax 6.2 0.0 346.9 0.30 0.6 ± 4.0 Overall 1.11 Appendix B. dSPECT Reconstructions 180 B . 3 E X P E R I M E N T A L R E S U L T S B.3.1 N O N A T T E N U A T E D E X P E R I M E N T S In the case of washout only, time activity curves were fit to monoexponential functions of the form, ln(2) t A(t) = A0e T ^ (B.5) In cases involving an additional washin component as well the TA curves, the data was fit with a dual exponential function, ( ln2 t In2 t \ e T l / 2 , i _ e T i /2 ,2 j (B.6) Appendix B. dSPECT Reconstructions 181 Experiment 1 - Single head, 180° rotation 0 5 10 0 5 10 True Dynamic Parameters Region A 0 A Tl/2 (total counts) (min - 1) (min) i 70770 0.12 5.84 ii 86687 0.20 3.54 iii 28236 0.28 2.45 iv 16637 0.12 5.58 Reconstructed Parameters Region A 0 A Tl/2 (total counts) (min - 1) (min) i 66339 0.12 5.78 ii 80069 0.20 3.55 iii 15661 0.21 3.25 iv 18867 0.16 4.43 Appendix B. dSPECT Reconstructions 182 Experiment 2 - Triple head, 180° rotation 0 5 10 0 5 10 • 0 5 10 0 5 10 True Dynamic Parameters Region A 0 A Tl/2 (total counts) (min - 1) (min) i 65581 0.11 6.47 ii 81013 0.18 3.82 iii 25634 0.27 2.61 iv 15234 0.11 6.41 Reconstructed Parameters Region A 0 A Ti/2 (total counts) (min - 1) (min) i 59903 0.11 6.51 ii 75807 0.18 3.84 iii 20213 0.24 2.84 iv 13800 0.11 6.32 Appendix B. dSPECT Reconstructions Experiment 3 - Dual head, 90° rotation x104 X104 1 3 i I 3 | — 0 5 10 0 5 10 0 5 10 0 5 10 0 5 10 Time (min} True Dynamic Parameters Region A , Ai Tl/2,1 A 2 A 2 Tl/2,2 Peak (counts) (min - 1) (min) (counts) (min - 1) (min) (min) i 25482 0.34 2.02 n/a n/a n/a n/a ii 22322 0.11 6.59 n/a n/a n/a n/a iii 21816 0.18 3.85 n/a n/a n/a n/a iv -7387 0.78 0.89 23063 0.08 8.2 1.2 V -8002 0.68 1.02 21439 0.04 18.3 3.0 Reconstructed Parameters Region Ai Tl/2,1 A 2 A 2 Tl/2,2 Peak (counts) (min - 1) (min) (counts) (min - 1) (min) (min) i 22812 0.33 2.11 n/a n/a n/a n/a ii 22283 0.11 6.59 n/a n/a n/a n/a iii 20890 0.18 3.83 n/a n/a n/a n/a iv -7132 0.87 0.80 21739 0.09 8.06 1.1 ± 0.3 V -8575 0.54 1.29 21796 0.04 15.63 3.8 ± 1.2 Appendix B. dSPECT Reconstructions 184 B . 3 . 2 A T T E N U A T E D E X P E R I M E N T S In attenuated experiments, true dynamic parameters were obtained from incoming water flow rates by the use of the equation, A(t) = Aoe-v* (B.7) Values of A0 were obtained from well counter measurements performed a short time prior to injection into the dynamic phantom and rescaled by the measured efficiency of the SPECT camera. Appendix B. dSPECT Reconstructions Dual head, 180° rotation x 10* 3 2.5 2 1.5 1 0.5 0 Region 1 Region 2 —— Region 3 Region 4 s J . A / I i I / . | . . / . \ I / / ' ' / • " " — 10 15 20 25 Time True Dynamic Parameters Region A 0 A Tl/2 Injection (total counts) (min - 1) (min) Time (min) 1 22032 0.173 4.0 1.0 2 19321 0.116 6.0 2.5 3 14009 0.058 12.0 5.0 4 16899 0.0 static 0.0 Reconstructed Parameters Region Ac A Ti/2 Peak Position (total counts) (min - 1) (min) (min) 1 25444 n/a n/a 1.8 2 8449 n/a n/a 4.9 3 21812 n/a n/a 5.7 4 14420 n/a n/a 4.9 Appendix B. dSPECT Reconstructions 186 Experiment 2 - Triple head, 180° rotation 4 3 * 0 5 10 15 20 25 Time True Dynamic Parameters Region A 0 A Tl/2 Injection (mCi) (min - 1) (min) Time (min) 1 13022 0.173 4.0 1.0 2 12877 0.116 6.0 2.5 3 12082 0.058 12.0 5.0 4 10990 0.0 oo 0.0 Reconstructed Parameters Region A 0 A Tl/2 Peak Position (total counts) (min - 1) (min) (min) 1 11704 n/a n/a 1.4 2 16997 n/a n/a 2.5 3 14901 n/a n/a 6.4 4 12502 n/a n/a 4.5 B I B L I O G R A P H Y J. Keyes. Clinical applications of SPECT. Int. J. Card. Imaging, 5:25-32, 1989. B. Mandalapu, M. Amato, and H. Stratmann. Technetium tc99m sestamibi myocar-dial perfusion imaging: current role for evaluation of prognosis. Chest, 115(6):1684-1694, 1999. H. Gallowitsch, O. Unterweger, P. Mikosch, E. Kresnik, J. Sykora, and P. Lind. Attenuation correction improves the detection of viable myocardium by thallium-201 cardiac tomography in patients with previous myocardial infarction and left ventricular dysfunction. Eur. J. Nuc. Med., 26(5):459-466, 1999. E.G. Depuey, S. Parmett, M. Ghesani, A. Rozanski, K. Nichols, and H. Salensky. Comparison of Tc-99m sestamibi and tl-201 gated perfusion SPECT. J. Nucl. Car-diol, 6(3):278-285, 1999. M. Verani. 201T1 myocardial perfusion imaging. Curr. Opin. Radiol, 3(6):797-809, 1991. S. Braat. 99mTc myocardial perfusion imaging. Curr. Opin. Radiol, 3:810-816, 1991. H. Iida, H. Itoh, M. Nakazawa, J. Hatazawa, H. Nishimura, Y. Onish, and K. Ue-mura. Quantitative mapping of regional cerebral blood flow using Iodine-123-IMP and SPECT. J. Nuc. Med, 35:2019-2030, 1994. M.D. Devous Sr., J.K. Payne, J.L. Lowe, and R.F. Leroy. Comparison of technetium-99m-ECD to xenon-133 SPECT in normal controls in patients with mild to moderate regional cerebral blood flow abnormalities. J. Nuc. Med., 34:754-761, 1993. A. Gottschalk. V / Q imaging and the diagnosis of PE: Can we shift the gray to black and white? J. Nuc. Med., 35(12):1931-1932, 1994. J. Juni. Lung scanning in the diagnosis of pulmonary embolism: the emperor re-dressed. Semin. Nuc. Med., 21:281-296, 1991. J. Ridge, J. Bading, Gelbard A, R. Benua, and J. Daly. Perfusion of colorectal hepatic metastases. Relative distribution of flow from the hepatic artery and portal vein. Cancer, 59:1547-1553, 1987. S. Kuiken, M . Samsom, M. Camilleri, B. Mullan, D. Burton, L. Kost, T. Hardy-man, B. Brinkmann, and M. O'Connor. Development of a test to measure gastric accomodation in humans. Am. J. Physiology, 277:G1217-G1221, 1999. 187 Bibliography 188 [13] J. Lee, W. Shih, and C. Chuang. Tc-99m RBC SPECT showing colonic bleeding in traumatic pseudoaneuryism with arterllcolonic fistula. Clin. Nucl. MEd., 23:397-399, 1998. [14] S. Stevens, T. Male, and J. Turner. Pelvic fractures diagnosed by bone scintigraphy in patients with normal radiographs after a fall. Med. J. Aust., 171(9):476-478, 1999. [15] G. Cook and I. Fogelman. The role of positron emission tomography in the man-agement of bone metastases. Cancer, 88:2927-2933, 2000. [16] Y. Nakamoto, T. Higashi, H. Sakahara, N. Tamaki, M. Kogire, M Imamura, and J. Konishi. Contribution of PET in the detectionof liver metastases from pancreatic tumours. Clin. Radiol., 54:248-252, 1999. [17] R. Bar-Shalom, A. Valdivia, and M. Blaufox. PET imaging in oncology. Semin. Nucl. Med., 30:150-185, 2000. [18] A. Peters. Scinitgraphic imaging of renal function. Exp. Nephoi, 6(5):391-397, 1998. [19] F. Gaspari, N. Perico, and G. Remuzzi. Measurement of glomerular filtration rate. Kidney Int. Suppl., 63:S151-154, 1997. [20] S.R. Deans. The Radon transform and some of its applications. John Wiley & Sons, Toronto, 1983. [21] C. Nahmias, D. Kenyon, K. Kouris, and E. Garnett. Understanding convolution backprojection. Single Photon Emission Computed Tomography and Other Selected Computer Topics, pages 19-29, 1980. [22] T.F. Budinger, G.T. Gullberg, and R.H. Huesman. Emission computed tomography Image Reconstruction From Projections. Springier-Verlag, New York, 1979. [23] G.T. Herman. Image Reconstruction from Projections. Academic, New York, 1980. [24] J.A. Fessler. Penalized weighted least-squares image reconstruction for positron emission tomography: IEEE Trans. Med. Imag., 13(2):290-300, 1994. [25] L.A. Shepp and Y. Vardi. Maximum liklihood reconstruction for emission tomogra-phy. IEEE Trans. Med. Imag., MI-L113-122, 1982. [26] J Kay. The E M algorithm in medical imaging; Stat. Methods Med. Res., 6:55-75, 1997. Bibliography 189 [27] D.W. Wilson and B.M.W Tsui. Noise properties of filtered-backprojection and ML-EM reconstructed emission tomographic images. IEEE Trans. Nuc. Sci., 40(4):1198-1203, 1993. [28] J. Juni. SPECT of rapidly cleared tracers: Imaging a Cheshire cat. J. Nuc. Med., 33:1206-1208, 1992. [29] J. Taki, K. Nakajima, I. Matsunari, N. Tonami, H. Sumiya, K. Ichiyanagi, S. Takada, and K. Hisada. Impaired 1-123 beta-methyl-branched fatty acid (BMIPP) uptake in relation to wall motion in ischemic but viable myocardium. Eur. J. Nuc. Med., 22:744, 1995. Abstract. [30] W.S. Richter, S. Beckmann, G. Armbrecht, T. Schuppenhauer, H. Eichstaedt, and R. Felix. Dynamic SPECT with I-123-phenyl-oentaaadecanoic acid (IPPA) in pa-tients with coronary artery disease — feasibility and first results. Eur. J. Nuc. Med., 22:744, 1995. Abstract. [31] S. Nellis, A. Liedtke, and B. Renstrom. Fatty acid kinetics in aerobic my-ocardium:Characteristics of tracer carbon entry and washout and influence of metabolic demand. J. Nuc. Med., 33:1864-1874, 1992. [32] H. Machulla, M. Marsmann, and K. Dutschka. Biochemical concept and synthesis of a radioiodinated phenylfatty acid for in vivo metabolic studies of the myocardium. Eur. J. Nuc. Med., 5:171-173, 1980. [33] J. Corbett. Fatty acids for myocardial imaging. Semin. Nucl. Med., 29:237-258, 1999. [34] M.J. van Eenige, F.C. Visser, C.M.B. Duwel, P.D. Bezemer, G. Wester a, A.J . P. Karreman, and J.P. Roos. Analysis of myocardial time-activity curves of 1 2 3 I -heptadecanoic acid. I. curve fitting. Nuklearmedizin, 26:241-247, 1987. [35] I. Matsunari, T. Saga, J. Taki, and et al. Kinetics of iodine-123-BMIPP in patients with prior myocardial infarction: Assessment with dynamic rest and stress images compared with stress thallium-201 SPECT. J. Nuc. Med., 35:1279-1285, 1994. [36] M.P.J. Hudon, D.M. Lyster, W.R.E. Jamieson, A.K. Qayumi, C. Sartori, and H.A. Dougan. The metabolism of (123I)-iodophenyl pentadecanoic acid in a surgically induced canine model of regional ischemia. Eur. J. Nuc. Med., 16:199-204, 1990. [37] T. Chua, H. Kiat, G. Germano, K. Takemoto, G. Fernandez, Y. Biasio, J. Fried-man, and D. Berman. Rapid back to back adenosine stress/rest technetium-99m teboroxime myocardial perfusion SPECT using a triple-detector camera. J. Nuc. Med., 34:1485-1493, 1993. Bibliography 190 [38] L.L. Johnson. Myocardial perfusion imaging with technetium-99m-teboroxime. J. Nuc. Med., 35:689-692, 1994. [39] R. Narra, A. Nunn, B. Kuczynski, T. Feld, P. Wedeking, and W. Eckelman. A neutral technetium-99m complex for myocardial imaging. J. Nuc. Med., 30:1830-1837, 1989. [40] J. Nickles, A. Nunn, C. Stone, and B. Christian. Technetium-94m-teboroxime: Syn-thesis, dosimetry and initial PET imaging studies. J. Nuc. Med., 34:1058-1066, 1993. [41] R. Stewart, M. Schwaiger, G. Hutchins, P. Chiao, K. Gallagher, N. Nguyen, N. Petry, and W. Rogers. Myocardial clearance kinetics of technetium-99m-SQ30217: A marker of regional myocardial blood flow. J. Nuc. Med., 31(7):1183-1190, 1990. [42] A.M. Smith, G.T. Gullberg, P.E. Christian, and F.L. Datz. Kinetic modeling of teboroxime using dynamic SPECT imaging of a canine model. J. Nuc. Med., 35:484-495, 1994. [43] M.K. O'Connor and D.S. Cho. Rapid radiotracer washout from the heart: Effect on image quality in SPECT performed with a single-headed gamma camera system. J. Nuc. Med., 33(6):1146-1151, 1992. [44] B. Bok, A. Bice, M. Clausen, D. Wong, and H. Wagner. Artifacts in camera based single photon emission tomography due to time activity variation. Eur. J. Nucl. Med., 13:439-442, 1987. [45] K. Nakajima, J. Taki, H. Bunko, M. Matsudaira, A. Muramori, I. Matsunari, K. Hisada, and T. Ichihara. Dynamic acquisition with a three-headed SPECT sys-tem: application to technetium 99m-SQ30217 myocardial imaging. J. Nuc. Med., 32:1273-1277, 1991. [46] S. Ross, A. Welch, G. Gullberg, and R. Huesman. An investigation into the effect of input function shape and image acquisition interval on estimates of washin for dynamic cardiac SPECT. Phy. Med. Biol, 42:2193-2213, 1997. [47] G. Gullberg, R. Huesman, S. Ross, E. Di Bella, G. Zeng, B. Reutter, P. Christian, and S. Foresti. Nuclear Cardiology: State of the art and Future Directions, chapter Dynamic cardiac single photon emission computed tomography. Mosby-Year Book, Inc, Philadelphia, PA, 1999. [48] E.V.R. Di Bella, G.T. Gullberg, A.B. Barclay, and R.L. Eisner. Circumferential profiling for region-based analysis of dynamic SPECT data. In 1996 IEEE Nuclear Science Symposium Conference Record, pages 1608-1612, November 1996. Bibliography 191 [49] B.E. Oppenheim and J.D. Krepshaw. Dynamic hepatobiliary SPECT: a method for tomography of a changing radioactivity distribution. J. Nuc. Med., 29:98-102, 1988. [50] T. Farncombe, A. Celler, D. Noll, J. Maeght, and R. Harrop. An evaluation of dy-namic SPECT imaging methods. In 1998 IEEE Nuclear Science Symposium Con-ference Record, November 1998. [51] M.N. Limber, A. Celler, J.S. Barney, M.A. Limber, and J.M. Borwein. Direct re-construction of functional parameters for dynamic SPECT. IEEE Trans. Nuc. Sci., 42:1249-1256, 1995. [52] G.L. Zeng and G.T. Gullberg. Estimating kinetic parameters directly from projec-tion measurements. IEEE Trans. Nuc. Sci, 42:2339-2346, 1995. [53] B. Reutter, G. Gullberg, and R. Huesman. Kinetic parameter estimation from attenuated SPECT projection measurements. IEEE Trans. Nuc. Sci., 45(6):3007-3013, 1998. [54] R. Huesman, B. Ruetter, G. Zeng, and G. Gullberg. Kinetic parameter estimation from SPECT cone-beam projection measurements. Phys. Med. Biol., 43(4):973-982, 1998. [55] E. Hebber, D. Oldenburg, T. Farncombe, and A. Celler. Direct estimation of dy-namic parameters in SPECT tomography. IEEE Trans. Nuc. Sci, 44(6):2425-2430, 1997. [56] J. Maltz, E. Polak, and T. Budinger. Multistart optimisation algorithm for joint spatial and kinetic parameter estimation in dynamic ECT. In 1998 IEEE Medical Imaging Conference Record, 1998. [57] A. Houston and W. Sampson. A quantitative comparison of some FADS methods in renal dynamic studies using simulated and phantom data. Phys. Med. Biol, 42:199-217, 1997. [58] J. Matthews, D. Bailey, P. Price, and V. Cunningham. The direct calculation of parametric images from dynamic PET data using maximum-likelihood iterative re-construction. Phys. Med. Biol, 42:1155-1173, 1997. [59] A. Sitek, G. Gullberg, E. DiBella, and A. Celler. Reconstruction of tomographic dynamic studies using a factor model. IEEE Trans. Med. Imag., 1999. submitted. [60] A. Sitek, E. DiBella, and G. Gullberg. Factor analysis with a priori knowledge -Application in dynamic cardiac SPECT. Phys. Med. Biol, 1999. submitted. Bibliography 192 [61] J. Maltz, B. Reutter, R. Huesman, and T. Budinger. Direct kinetic parameter es-timation from dynamic ECT sinogram using dimension-reduced time-activity basis. In 1999 IEEE Medical Imaging Conference Record, 1999. [62] R.J. Jaszczak, D.R. Gilland, M.W. Hanson, S. Jang, K.L. Greer, and R.E. Cole-man. Fast transmission ct for determining attenuation maps using a collimated line source, rotatable air-copper-lead attenuators and fan-beam collimation. J. Nuc. Med., 34:1577-1586, 1993. [63] M. King, B. Tsui, T.S. Pan, S.J. Glick, and E.J. Soares. Attenuation compen-sation for cardiac single-photon emission computed tomographic imaging: Part 2. attenuation compensation algorithms. J. Nuc. Card., 3:55-64, 1996. [64] T. Morozumi, M. Nakajima, K. Ogawa, and S. Yuta. Attenuation correction meth-ods using the information of attenuation distribution for single photon emission ct. Med. Imag. Technoi, 2:20-28, 1988. [65] C.H. Tung, G.T. Gullberg, G.L. Zeng, P.E. Christian, F.L. Datz, and H.T. Morgan. Non-uniform attenuation correction using simultaneous transmission and emission converging tomography. IEEE Trans. Nuc. Sci., 39:1134-1143, 1992. [66] A. Celler, A. Sitek, E. Stoub, P. Hawman, R. Harrop, and D. Lyster. Investigation of an array of multiple line sources for SPECT transmission scans: Simulation, phantom and patient studies. J. Nuc. Med., 39(12):2183-2189, 1998. [67] E.P. Ficaro, J.A. Fessler, P.D. Shreve, J.N. Kritzman, P.A. Rose, and J.R. Corbett. Simultaneous transmission/emission myocardial perfusion tomography, diagnostic accuracy of attenuation-corrected 99mtc-sestamibi single-photon emission computed tomography. Circulation, 93:463-473, 1996. [68] M.A. King, B.M.W. Tsui, and T.P. Pan. Attenuation compensation for cardiac single-photon emission computed tomographic imaging: Part 1. impact of attenu-ation and methods of estimating attenuation maps. J. Nuc. Card., 2(6):513-524, 1995. [69] F.J.Th. Wackers. Attenuation correction, or the emperor's new clothes? J. Nucl. Med., 40:1310-1312, 1999. [70] T. Farncombe, A. Celler, D. Noll, J. Maeght, and R. Harrop. The incorporation of organ uptake into dynamic spect (dSPECT) image reconstruction. IEEE Trans. Nuc. Sci., 2000. submitted. [71] D. Liu and J. Nocedal. On the limited memory BFGS method for large scale opti-mization. Math. Prog. B, 45:503-528, 1989. Bibliography 193 [72] J. Maeght. Analyse et methodes pour un probleme inverse en tomographic dy-namique. PhD thesis, Universite Paul Sabatier, 2000. [73] H.M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag., 13:601-609, 1994. [74] H.H. Barrett, D. Wilson, and B. Tsui. Noise properties of the E M algorithm - I. theory. Phys. Med. Biol, 39:833-846, 1994. [75] D.L. Snyder, M.I. Miller, L.J. Thomas, and D.G. Polite. Noise and edge artifacts im maximum-likelihood reconstruction for emission tomography. IEEE Trans. Med. Imag., 6:228-238, 1987. [76] D.L. Snyder and M.I. Miller. The use of sieves to stabilize images produced with the E M algorithm for emission tomography. IEEE Trans. Med. Imag., 35:3864-3871, 1985. [77] A. Celler, T. Farncombe, R. Harrop, and D. Lyster. Dynamic heart-in-thorax phan-tom for functional SPECT. IEEE Trans. Nuc. Sci., 44(4):1600-1605, 1997. [78] T. Farncombe. The design and development of a phantom for use in dynamic SPECT imaging. Master's thesis, University of British Columbia, 1998. [79] F. Beekman, C. Kamphuis, and E. Frey. Scatter compensation methods in 3D iterative SPECT reconstruction: A siulation study. Phys. Med. Biol., 42:1619-1632, 1997. [80] E. Frey, Z.W. Ju, and B.M.W. Tsui. Modeling the scatter response function in inhomogeneous scattering media for SPECT. IEEE Trans. Nuc. Sci., NS-4L1585-1595, 1994. [81] D. Kadrmas, E. Frey, and B. Tsui, application of reconstruction-based scatter com-pensation to Thallium-201 SPECT: Implementations for reduced reconstructed im-age noise. IEEE Trans. Med. Imag., 17:325-333, 1998.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Functional dynamic spect imaging using a single slow...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Functional dynamic spect imaging using a single slow camera rotation Farncombe, Troy 2000
pdf
Page Metadata
Item Metadata
Title | Functional dynamic spect imaging using a single slow camera rotation |
Creator |
Farncombe, Troy |
Date Issued | 2000 |
Description | Dynamic single photon emission computed tomography (SPECT) is a relatively new imaging method that uses radioactive tracers and tomographic data acquisition techniques in order to quantify temporal changes in regional radiotracer concentrations within a patient. This is important as the rate of change in tracer concentration within an organ often can be related to the functional ability of that organ. In this work, a new method is presented that is able to determine these kinetic rates while using a conventional single or multiple detector SPECT camera, and more importantly, a single, slow camera rotation in the data collection process (herein this reconstruction method will be referred to as dSPECT). This reconstruction method is based on the fact that a temporal change in the activity concentration at a given location can be represented by a linear inequality constraint over time. Two iterative reconstruction algorithms, constrained least squares (CLS) and dynamic expectation maximization (dEM), have been tested using this approach with a variety of computer simulations and phantom experiments. In simulations involving a slow dynamic change of activity, results indicate that the dSPECT reconstruction procedure typically produces kinetic parameter estimates with a 7% error when using projection data acquired with a single 180° rotation of a triple headed SPECT camera system. This error increases to about 15% for data acquired with a dual head system and further increases to about 50% for single detector acquisitions. When simulated with faster dynamic parameters, errors increased slightly to about 8%, 12% and 55% for acquisitions involving triple, dual and single head systems respectively. |
Extent | 27330709 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085703 |
URI | http://hdl.handle.net/2429/11070 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-565394.pdf [ 26.06MB ]
- Metadata
- JSON: 831-1.0085703.json
- JSON-LD: 831-1.0085703-ld.json
- RDF/XML (Pretty): 831-1.0085703-rdf.xml
- RDF/JSON: 831-1.0085703-rdf.json
- Turtle: 831-1.0085703-turtle.txt
- N-Triples: 831-1.0085703-rdf-ntriples.txt
- Original Record: 831-1.0085703-source.json
- Full Text
- 831-1.0085703-fulltext.txt
- Citation
- 831-1.0085703.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085703/manifest