UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Record Information

Full Text 
A MULTIRESOLUTION RESTORATION METHOD FOR CARDIAC SPECT By JUAN M. FRANQUIZ A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1997 ACKNOWLEDGMENTS The author expresses special thanks to Dr. Mario Ariet for his support and for providing access to the computer facilities at the Division of Computer Sciences of the College of Medicine, University of Florida, for this research work. The author also thanks Mr. Stefan Lupkiewicz for providing technical help in installing and using part of the software employed in this research. The author is grateful to his academic advisor, Dr. Shailendra Shukla, for the support, criticisms, technical advice and assistance offered throughout this research. The author also would like to express his gratitude to the staff of the Department of Nuclear Medicine at the Veterans Administration Medical Center in Gainesville, for supporting the experimental part and providing access to their nuclear instruments and radioactive sources. The author is specially grateful to Dr. Michael Fagein for his comments and Mr. Morris Thompson for preparing the radioactive sources used in this research and his general help in the experimental work. The author thanks Dr. Nils Diaz and Dr. David Hintenlang for their support in the beginning of the graduate course work at the University of Florida. The author also thanks Dr. Andrew Laine of the Department of Computer and Information Science and Engineering for his initial comments and revision of the proposal of this research and for providing part of the software used in this work. Finally, the author wishes to thank ii supervisory committee members Dr. Wesley Bolch, Dr. Walter Drane, Dr. Janice Honeymann and Dr. William Properzio for their assistance and the time dedicated to the revision, comments, criticisms and recommendations throughout this research. iii TABLE OF CONTENTS Page ACKNOWLEDGMENTS .................................................... ii LIST OF TABLES ..................................................... vi LIST OF FIGURES .................................................... viii ABSTRACT ........................................................... xi CHAPTERS 1 INTRODUCTION ................................................. 1 The Problem ....................................... .............. 3 Hypothesis and Objectives ........................................ 8 Content of Chapters ............................................. 9 2 CARDIAC SPECT .................................... ............... 11 Radiopharmaceuticals ..... ...................................... 11 SPECT Data Acquisition ........................................ 14 SPECT Reconstruction .......................................... 19 Myocardial Tomographic Slices ................................. 25 Degrading Factors in SPECT .................................... 28 Artifacts in the Inferior Myocardial Wall ....................... 45 3 THE RESTORATION PROBLEM IN SPECT ................................ 48 Degradation Models in SPECT ..................................... 49 Restoration Filters .............................................. 64 Computer Implementation of the Metz Filter ..................... 72 4 MULTIRESOLUTION RESTORATION ALGORITHM ............................ 87 Wavelet Multiresolution Representation .......................... 89 Multiresolution Restoration ................................... 99 5 RESTORATION IN A REALISTIC CARDIAC CHEST PHANTOM .................. 111 Material and Methods ......................................... 112 Results and Discussion ......................................... 119 Conclusions .................................................... 156 iv 6 CONCLUSIONS AND FUTURE WORK ...................................... 159 Other Potential Applications ................................... 161 Future Developments .............................................. 163 APPENDIXES A COMPUTER IMPLEMENTATION OF A METZ FILTER ........................ 170 Al Calculation of MTF Parameters .............................. 170 A2 Simulation of a CardiacLiver Projection ..................... 175 A3 Metz Filtering ........................... ................ 177 B MULTIRESOLUTION DECOMPOSITION AND RESTORATION ALGORITHM .......... 179 Bl Multiresolution Decomposition .............................. 179 B2 Multiresolution Restoration ................................. 181 C MULTIRESOLUTION RESTORATION OF SPECT PROJECTIONS ................. 184 Cl Restoration of SPECT Projections ........................... 184 C2 Analysis Programs ...................................... 187 REFERENCES .................................................... 192 BIOGRAPHICAL SKETCH ........................................... 202 V LIST OF TABLES Table Page 2.1 Contribution of outofplane scattered photons per voxel as a function of the activity ratio of outofplane to emission voxels (ROE) for different neighboring planes ................. 42 3.1 Restoration Models in SPECT and their approximate solutions...... 53 3.2 Line spread function (LSF) and point spread function (PSF) parameters fitted to experimental data from line sources ........ 75 3.3 Analytical expressions of fitted modulation transfer functions .. 77 3.4 Normalized mean square error (NMSE) as a function of the power factor X in image restoration using the Metz filter ............. 83 5.1 Quantitative analysis of myocardial bull's eye polar maps ....... 123 5.2 Quantitative analysis of a midventricular shortaxis slice ...... 125 5.3 Normalized chisquare of bull's eye polar maps as a function of the attenuation coefficient and livertoheart activity ratio (LHAR) ............................................... 129 5.4 Uniformity of bull's eye polar maps as a function of the attenuation coefficient and livertoheart activity ratio (LHAR) for the cardiacliver phantom in water ................. 130 5.5 Normalized chisquare calculated with different reconstruction protocols as a function of the livertoheart activity ratio .... 135 5.6 Uniformity of bull's eye polar maps calculated using different reconstruction protocols ....................................... 137 5.7 Quantitative analysis of a midventricular shortaxis slice of the cardiac insert in water without liver activity .............. 141 5.8 Quantitative analysis of bull's eye polar maps of the cardiac liver insert in water ............................................ 142 vi 5.9 Quantitative analysis of the effect of multiresolution restoration on the normalized chisquare of bull's eye polar maps for different livertoheart activity ratio (LHAR) and reconstruction protocols ...................................... 143 5.10 Quantitative analysis of the effect of multiresolution restoration on uniformity of bull's eye polar maps as a function of the livertoheart activity ratio for different reconstruction protocols ...................................... 146 5.11 Contrast of myocardial defect measured in a shortaxis slice for different reconstruction protocols ........................ 148 5.12 Results of geometric calculation of true defect size on the bull's eye polar map ........................................ 150 5.13 Defect size in percent of the polar map area calculated from the histogram of counts of bull's eye polar maps as a function of the threshold level ..................................... 150 5.14 Threshold values for calculating defect size in cardiac SPECT.........................................................154 6.1 Data of the three patients included in a preliminary evaluation of the multiresolution restoration algorithm in attenuation corrected cardiac SPECT scans ................................... 165 vii LIST OF FIGURES Figure Page 21 Density function g(x,y) in the transaxial plane XY and its projection P,[r] at an angle 9 ................................ 20 22 Filtered backprojection technique .............................. 21 23 SPECT reconstruction convolution kernels and their frequency response ............................................... ......... 26 24 Myocardial tomographic slices .................................. 27 25 Detection geometry of unscaterred and outofplane scattered photons.........................................................38 26 Plot of the relative number of outofplane scattered photons per voxel as a function of scattering angle ..................... 41 31 Coordinates for source distribution and detection planeD .................. ....................................... 57 32 Line spread function (LSF) as a function of the distance from the face of the collimator (depth) in air and in a water tank of 20 cm in diameter .......................................... 63 33 Plots of modulation transfer function (MTF) at three different distances from the face of the collimator ....................... 65 34 Plots of Metz filter for power factors of 10, 4 and 2 ............. 70 35 Fitted line spread functions and experimental counts ............ 76 36 Plots of modulation transfer functions calculated from fitted line spread functions .......................................... 79 37 Simulated cardiacliver projection .............................. 81 38 Plots of normalized mean square error (NMSE) as a function of Metz power factor X for projections with different noise levels..........................................................84 viii 41 The "mexican hat" function and its translations and dilatations (wavelets) ......................................... 93 42 Frequency response of analyzing functions and MTF .............. 101 43 Multiresolution decomposition of a high resolution MRI brain image into 5 resolution levels ........................... 102 44 Multiresolution decomposition of a low clinical brain SPECT projection into 5 resolution levels ...................... 103 45 Square 12 norm expressed in percent as a function of the subband image m ................... ............................. 107 46 Count profiles through the simulated projection ................ 110 51 Decomposition of a projection into 5 images or resolution levels......................................................... 115 52 Bull's eye polar maps of the cardiac phantom without defects and LHAR = 0 .............................................. 120 53 Histograms of polar map counts for the cardiac insert without liveractivity .................................................. 122 54 Horizontal profile lines through the center of a midventricular shortaxisslice ................................................ 126 55 Graphs depicting the effect of variation in the attenuation coefficient value upon the uniformity of reconstructed bull's eye polar maps ......................................... 131 56 Graphs depicting the effect of variation in the attenuation coefficient value upon the normalized chisquare of polar maps ... 132 57 Normalized chisquare of bull's eye polar maps as a function of the livertoheart activity ratio (LHAR) and cutoff frequency for four different reconstruction protocols .................... 136 58 Bull's eye polar maps of the cardiac liver insert (LHAR = 3.5) inwater ......................................................... 138 59 Graphs depicting the effect of restoration on uniformity of bull's eye polar maps as a function of the livertoheart activity ratio (LHAR) ........................................ 144 510 Effect of restoration on bull's eye polar maps ................ 145 511 Myocardial tomographic slices for the cardiac insert in water with apical and inferoapical defects ......................... 149 ix 512 Graphs depicting the calculated defect size as a function of the threshold level ......................................... 151 61 Bull's eye polar maps of patient number 3 ...................... 166 x Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A MULTIRESOLUTION RESTORATION METHOD FOR CARDIAC SPECT By Juan M Franquiz December, 1997 Chairman: Shailendra Shukla, Ph.D. Major Department: Nuclear and Radiological Engineering Singlephoton emission computed tomography (SPECT) is affected by photon attenuation and image blurring due to Compton scatter and geometric detector response. Attenuation correction is important to increase diagnostic accuracy of cardiac SPECT. However, in attenuationcorrected scans, scattered photons from radioactivity in the liver could produce a spillover of counts into the inferior myocardial wall. In the clinical setting, blurring effects could be compensated by restoration with Wiener and Metz filters. Inconveniences of these procedures are that the Wiener filter depends upon the power spectra of the object image and noise, which are unknown, while Metz parameters have to be optimized by trial and error. This research develops an alternative restoration procedure based on a multiresolution denoising and regularization algorithm. It was hypothesized that this representation leads to a more straightforward and xi automatic restoration than conventional filters. The main objective of the research was the development and assessment of the multiresolution algorithm for compensating the liver spillover artifact. The multiresolution algorithm decomposes original SPECT projections into a set of subband frequency images. This allows a simple denoising and regularization procedure by discarding high frequency channels and performing inversion only in low and intermediate frequencies. The method was assessed in bull's eye polar maps and shortaxis attenuationcorrected reconstructions of a realistic cardiacchest phantom with a custommade liver insert and different 9"Tc livertoheart activity ratios. Inferior myocardial defects were simulated in some experiments. The cardiac phantom in free air was considered as the gold standard reference. Quantitative analysis was performed by calculating contrast of shortaxis slices and the normalized chisquare measure, defect size and mean and standard deviation of polar map counts. The performance of the multiresolution method was also assessed in 201Tl clinical SPECT studies. Phantom results demonstrated that the multiresolution algorithm compensated the liver spillover artifact, yielded uniform polar maps and improved significantly the accuracy in calculating myocardial defect size. The procedure does not require operator intervention and can be easily implemented in the clinical setting by using Fourier transform techniques. Finally, the extension of the multiresolution method to other SPECT procedures is discussed and recommended. xii CHAPTER 1 INTRODUCTION Singlephoton emission computed tomography (SPECT) is the clinical procedure most commonly used in nuclear medicine. The objective of SPECT is the quantitation of the threedimensional (3D) distribution of a radiopharmaceutical in an organ or region of interest. This is accomplished by calculating, from a set of images (projections) about the human body, the radioactivity concentration in twodimensional (2D) transaxial slices. Radiopharmaceuticals are usually labeled with gamma emitters with emissions between 80 keV and 300 keV. The projection images are acquired with one or more rotating scintillation or gamma cameras. Quantitation in SPECT studies includes two broad categories of experimental procedures. The first is absolute quantitation, which means the accurate and precise measurement of certain quantities in absolute units. Examples are measurements of organ volumes and concentration of the radiotracer expressed in units of activity per unit of volume or mass. Absolute quantitations are needed in absorbed dose calculations for target and normal tissues when radionuclides are administered with therapeutic purposes (Zanzonico et al., 1989; Qian and Clarke, 1996). The other category is relative quantitation, in which relative comparisons of the radiotracer uptake in two regions or at different times are needed. This is the type of quantitation used in diagnostic studies. Examples are the assessment of brain SPECT where the radiotracer uptake in regions of the 1 2 left and right hemispheres are compared, or myocardial SPECT where the radiotracer uptake in the cardiac muscle with the patient at rest is compared with the uptake during a physical or pharmacological stress (Berman et al., 1991; Leppo, 1996). Absolute quantitation could be an ideal goal in SPECT studies and has demonstrated great utility in physiological research (Maddahi and Czermin, 1996) and dosimetry calculations (Zanzonico et al., 1989). However, experience has demonstrated that absolute quantitation in SPECT is not required for clinical studies. The major use of SPECT in nuclear medicine is for diagnostic purposes. Thus, relative quantitation is a much more relevant problem in SPECT than absolute quantitation. In addition, the particular problem studied in this research is related to myocardial SPECT in which absolute measurements are not needed (Maddahi and Czermin, 1996). Therefore, absolute quantitation problems in SPECT will not be treated in this research. In all the material that follows, SPECT quantitation means relative quantitation. The 3D reconstruction problem in SPECT has been approached using iterative algorithms based on parameter estimation techniques (Shepp and Vardi, 1982; Tsui et al., 1994a; 1994b) and noniterative Fourier transform based methods (Larssen, 1980). Iterative algorithms have the potential advantage that accurate quantitative corrections for the three major degrading factors in SPECT (photon attenuation, Compton scatter and finite detector geometric resolution) can be included into the reconstruction process (Buvat et al., 1994; Tsui et al., 1994a). However, these corrections are done at the expense of such large hardware and computational requirements, that iterative algorithms have remained 3 unsuitable for routine clinical use. Recently, some commercial systems (Case et al., 1996; Cullom et al., 1996; McCartney et al., 1996) have included iterative algorithms but only for attenuation correction with results which are simpler than that obtained with scatter and geometric detector response compensations. Fourier transform based methods are computationally efficient but they have the limitation that accurate corrections can not be included into the reconstruction process (Larssen, 1980). Any correction for SPECT degrading factors has to be performed before or after reconstruction (Chang, 1978; Gullberg and Budinger, 1981; King et al., 1987; Buvat et al., 1994). Among these methods, the filtered backprojection algorithm (Larssen, 1980) is the most easily implemented and so far it is the only SPECT reconstruction technique used routinely in clinical studies. Therefore, because this research is focused on the analysis and solution of a particular clinical problem, only reconstruction and correction methods based on the filtered backprojection technique are considered. The Problem SPECT quantitation is significantly affected by photon attenuation and Compton scatter. Both effects can induce inaccuracies in clinical diagnosis, mostly in myocardial SPECT studies (Nuyts et al., 1995; O'Connor et al., 1995a; Ficaro et al., 1996; Maniawski et al., 1996). Because tissue attenuation is not uniform, the attenuation of photons in some regions produces a larger reduction in the number of counts than in other regions (Tsui et al., 1994a; 1994b). It is known that in 9"mTc SPECT studies, the attenuation effect in some regions can reduce the detected 4 counts to 25% of the unattenuated counts (Tsui et al., 1994a). This artifactual reduction of counts can be interpreted as heterogeneous or pathological uptake distributions. A low clinical specificity (high number of false positive studies) could be obtained if apparent uptake defects due to attenuation are attributed to a true reduced radiotracer uptake. Similarly, a reduced clinical sensitivity (low number of true positive studies) could be obtained if true radiotracer uptake defects are attributed to attenuation. Another important degrading effect is Compton scatter (O'Connor et al., 1995a; Welch et al., 1995) which produces blurred images and degradation of lesion contrast (Buvat et al., 1994). Consequently, radiotracer defect severity and extent can be obscured or can appear within normal limits because of the scattered radiation that reduces clinical sensitivity. An additional and very important blurring effect is due to the limited spatial resolution of the detector (Tsui et al., 1994a). Photon attenuation as well as Compton scattering depends on the photon energy, the body size of the patient and its composition. This is why they are considered patientdependent effects (Tsui et al., 1994a). In addition, the physiological distribution of the radiotracer also influences the magnitude of the scatter effect. This is because the scatter contribution depends on the source distribution (Buvat et al., 1994). Finally, Compton scattering is also dependent on the width and offset of the energy window, the energy resolution of the imaging system and the acceptance solid angle of the detector collimator (O'Connor et al., 1995a; 1995b). 5 It is a widely recognized fact that photon attenuation is the major factor that affects the accuracy of quantitations in SPECT (Tsui et al., 1994a; Ficaro et al., 1996). This effect has been an important source of false positive studies of coronary artery disease in myocardial SPECT studies (Ficaro et al., 1996). Although myocardial SPECT without attenuation and scatter corrections has demonstrated high sensitivity in detecting coronary heart disease (Leppo, 1996), its specificity in assessing the absence of disease has been reported as suboptimal in some groups of patients (Ficaro et al., 1996). This has been attributed to the photon attenuation in the thorax which produces false uptake defects in patients with normal myocardial perfusion (Tsui et al., 1994a; Ficaro et al., 1996; McCartney et al., 1996). As a result of intensive investigative efforts during the last decade, the new commercial SPECT systems includes hardware and software for correcting attenuation in myocardial studies (Case et al., 1996; Cullom et al., 1996; Ficaro et al., 1996; McCartney et al., 1996). Although these systems are in the process of clinical validation, preliminary data have demonstrated a significant improvement in specificity when myocardial SPECT is attenuationcorrected (Ficaro et al., 1996; McCartney et al., 1996). However, the spillover of activity from nearby organs into the inferior myocardial wall resulting from photon scatter has become a significant problem in attenuationcorrected cardiac SPECT (Ljunberg et al., 1994; King et al., 1995; Nuyts et al., 1995; Cullom et al., 1996; Ficaro et al., 1996; McCartney et al., 1996). This effect could obscure myocardial defects in the inferior wall yielding false negative results (lower sensitivity). 6 The main source of scattered photons in myocardial SPECT is the radiotracer uptake in liver and bowel (Case et al., 1996; King et al., 1995; Nuyts et al., 1995). The effect of photon scatter from the liver has been a recognized problem in the assessment of inferior myocardial wall perfusion in some patients, even without attenuation correction (Beller and Watson, 1991; Heo et al., 1992; Johnson, 1994). However, this undesirable effect is more severe in attenuationcorrected studies (Ljunberg et al., 1994; King et al., 1995; Nuyts et al., 1995). The reason is that attenuation correction methods can not distinguish between unscattered and scattered photons. Therefore, scattered photons are artificially amplified with attenuation correction. Although the spillover effect has been reported in 20T1 myocardial SPECT (Case et al., 1996; McCartney et al., 1996), the effect is particularly important in studies with "9eTclabeled compounds which have a higher liver and gall bladder uptake (Beller and Watson, 1991; Heo et al., 1992; Johnson et al., 1994). The general consensus is that some scatter compensation has to be used in conjunction with attenuation correction in order to solve the liver spillover artifact. Different methods have been proposed to compensate for scatter in SPECT (Buvat et al., 1994), but only the simplest approaches can be used in the clinical setting because of constraints in computer capacity and processing time. Although there is not a completely satisfactory solution to the scatter problem, even the simplest methods could provide adequate compensation, at least for the accuracy required in most clinical studies (King et al., 1987; Frey and Tsui, 1995). 7 One common approach for compensating blurring due to both scatter and detector geometric response is restoration by deconvolution of SPECT data with a representative point spread function of the imaging system (King et al., 1983; 1987; Frey and Tsui, 1995). Assuming a linear and shiftinvariant imaging system, the restoration operation is usually performed in the spatial frequency domain, where the deconvolution operator is the inverse modulation transfer function of the imaging system (King et al., 1987). However, due to the singularity of the inverse modulation transfer function and the illconditioned nature of the restoration problem, the presence of even small amounts of noise can lead to unacceptable restorations (Andrew and Hunt, 1977, p.114). Regularization approaches which overcome this drawback have led to the Wiener and Metz restoration filters in nuclear medicine (King et al., 1983; 1987; 1991a). The Wiener filter depends on the power spectra of the object image and noise, which are unknown (King et al., 1987; Penney et al., 1990). The Metz filter is easier to implement but there is not a general and consistent methodology for determining its optimal parameters (King et al., 1987; 1991). In addition, both types of filters are image dependent. Therefore, the problem of SPECT restoration can be summarized as the optimal design of a lowpass frequency filter for denoising and regularization operations. Optimal filter design is important because the tradeoff between blurring and noise suppression depends critically on filter parameters (King et al., 1987; 1991). This is a common problem in single scale representations where it is difficult to completely separate noise from significant image features (Laine et al., 1994). 8 An alternative and relatively new restoration procedure based on multiresolution analysis by wavelet expansion (Daubechies, 1992, p. 7  11) is proposed in this research. The advantage of multiresolution representation is that it provides a general spacefrequency framework in which significant image features are more easily identified and segmented from noise (Laine et al., 1994). Projections can be decomposed into a set of subband images corresponding to frequency channels of different bandwidth which allow a combined spacefrequency representation of images. In the case of SPECT projections, high frequency or noisy subbands can be eliminated while restoration is performed only in those subband images where significant image features are present. This approach has been previously proposed by Wang et al. (1995) for general image restoration and Qian and Clarke (1996) for restoration of bremsstrahlung images of beta emitters. In this research, a multiresolution restoration algorithm is developed for deblurring cardiac SPECT images. The practical contribution of this new approach is that denoising and regularization operations can be performed in a more straightforward and automatic way, without trial anderror optimization (Metz restoration) or power spectra estimation (Wiener filtering). An automatic and general restoration methodology, which does not require operator intervention, would greatly help its practical implementation and application in the clinical setting. Hypothesis and Objectives The ultimate goal of this research is the development and critical assessment of a multiresolution regularization and denoising algorithm for compensation of the liver spillover artifact in attenuationcorrected 9 cardiac SPECT. The investigation is based on the hypothesis that a multiresolution restoration algorithm can compensate scatter and geometric detector response in noniterative cardiac SPECT reconstructions and with sufficient accuracy to be used in clinical studies. The demonstration of this hypothesis and the assessment of the multiresolution restoration algorithm were performed using a realistic cardiacchest phantom in which a custommade liver insert was included. The methods, algorithms and experimental results are presented in the chapters that follows. All the source code of the software written specifically for this research is included in appendixes. Content of Chapters Chapter 2 describes in detail the methodology of cardiac SPECT reconstructions, correction methods usually employed in the clinical setting and the physical degrading effects that could produce clinical artifacts. Special emphasis is dedicated to those artifacts that are due to the scattering of photons emitted within the liver in attenuation corrected SPECT scans. Chapter 3 discusses the restoration problem in SPECT and the mathematical models in which restoration is based. The discussion starts with the most realistic and complete approach, which is the nonlinear shiftvariant model. Then, it is demonstrated how practical limitations in the clinical implementation of complex models have led to approximate solutions in which the imaging system is represented by a linear and shiftinvariant model. This chapter also includes a discussion on restoration filters, the implementation and validation of the Metz filter, and the methodology and results of the experimental determination of the modulation transfer functions used in this research. Chapter 4 10 briefly presents the theory behind multiresolution wavelet representations and discusses its application in restoration of SPECT projections. This Chapter also describes the algorithms used in this research for multiresolution decomposition and restoration, the criteria for denoising and restoration and their validation in a computer simulated cardiac projection. Chapter 5 presents the experimental results obtained with the application of multiresolution restoration to attenuationcorrected SPECT reconstructions of a realistic cardiacchest phantom with different *""Tc livertoheart activity ratios. Experiments simulated the liver spillover artifact with myocardial homogeneous distribution of the radiotracer (normal subject) and with the inclusion of inferior myocardial defects in the cardiac phantom (pathologic lesions). The performance of the multiresolution algorithm was compared with that of the conventional Metz restoration. Finally, Chapter 6 summarizes the conclusions of this research, discusses a preliminary application of the multiresolution restoration method to clinical studies and makes recommendations for future clinical and physical work. Special emphasis is dedicated to the extension of the multiresolution restoration method to PET and SPECT using a linear shiftvariant restoration model. CHAPTER 2 CARDIAC SPECT The objective of cardiac SPECT is the relative quantitation and 3D display of the distribution of a radiopharmaceutical in the cardiac muscle. Myocardial radiopharmaceuticals reach the viable myocardium by a blood flow dependent process and tomographic slices provide an accurate indication of the amount of tissue perfusion from the coronary arteries. This information plays an important role in the diagnosis, followup and risk stratification of patients with suspected or known coronary artery disease (CAD), which constitutes the single most common cause of death in the United States (Cullom, 1995). In the last decade, major improvements have been made in radiopharmaceuticals, instrumentation and processing methods in order to improve accuracy in cardiac SPECT (Galt, 1994). Nevertheless, some inaccuracies persist as a consequence of the photon absorption and scatter in tissue. This chapter reviews the basic principles of cardiac SPECT, the newest systems, the physical degrading factors and the correction methods mostly used in routine clinical settings. Specific artifacts encountered in the inferior myocardial wall are discussed in more detail. Radiopharmaceuticals Cardiac SPECT has been performed using 20T1 for almost 15 years. This is a monovalent cation which is taken up by the viable myocardium by 11 12 a blood flow and passive transport process (Beller, 1994). The myocardial uptake of 20T1 starts immediately after its intravenous administration and also begins to redistribute into other tissues and the myocardium itself. Significant change in distribution is observed after more than four hours (Garcia et al., 1985). An initial cardiac SPECT depicts the myocardial perfusion at the time of radiotracer injection. A later cardiac SPECT reflects the initial perfusion and the rate of regional radiotracer loss. This biokinetics behavior has been used in the study of transient changes in myocardial perfusion. Clinical protocols usually include two cardiac SPECT scans for each patient (Mahmarian et al., 1993; Matsunari et al., 1996a). In the first study, the radiotracer is administered at the peak of a pharmacologic or exercise stress. The second one is performed at rest and after redistribution has occurred. If both stress and rest cardiac SPECT scans show homogeneous distribution of the radiotracer, they are interpreted as normal studies. Regions with significant decreased activity (regional defects) are considered abnormal. If the regional defect is present at rest, this is interpreted as suggestive of myocardial infarction. If the regional defect is present after stress but is either not present or less apparent at rest, then it is considered as indicative of regional ischemia. Unfortunately, the physical properties of 201T1 are less than ideal for cardiac SPECT imaging. Thallium201 decays by electron capture to 201Hg emitting characteristic xrays in the 60 keV to 83 keV range. This low energy produces low quality images. Therefore, accurate quantitations are very difficult to make because of photon absorption and scatter in tissue. In addition, the relatively long halflife of 201T1 (73 hours) does not 13 allow the administration to the patient of high activities, resulting in noisy SPECT images. On the contrary, 99mTc, the most widely used radionuclide in nuclear medicine, has a monoenergetic gamma emission of 140 keV which is ideal for gamma camera imaging. The shorter halflife (6 hours) of 99mTc allows the administration of activities that are up to 10 times higher than those of 201T1. This yields higher quality images in a shorter time period. Another advantage is that 99mTc is available from a "9Mo""Tc generator 24 hours a day, while 201T1 is cyclotrongenerated and requires offsite delivery. These superior properties of 99mTc motivated in past decades the search for different compounds to be used as 99"Tc labeled cardiac perfusion agents. In December 1990, the Food and Drug Administration approved the clinical use of two of those compounds: 9mTc sestamibi (DuPont Nemours, Inc.) and 99mTcteboroxime (Squibb Diagnostic, now Bracco Diagnostics). Both have shown successful results in clinical trials (Beller and Watson, 1991; Berman et al., 1991; Johnson, 1991; Van Train et al., 1994). More recently, a new compound, ""mTctetrofosmin (Amersham International, plc), has also been approved and is used commercially (Braat et al., 1994; Matsunari et al., 1996a). The compound 99mTcsestamibi is a monovalent cation that is taken by the viable myocardium by a passive transport process, similar to the uptake of 201T1 (Beller and Watson, 1991). This compound has a long myocardial residence time (approximately 5 hours) and little redistribution (Berman et al., 1991). In contrast, S""Tcteboroxime is a neutral lipophilic compound that diffuses rapidly across phospholipid membranes in a manner similar to freely diffusible radiotracers such as 14 '33Xe (Johnson, 1991). Its myocardial washout is flow dependent and shows a rapid biexponential clearance with halftimes of 3 to 6 minutes and 60 minutes for the early and late components, respectively (Johnson, 1994). The myocardial agent 99mTctetrofosmin is a lipophilic cationic compound with a rapid heart uptake and relatively slow clearance (Braat et al., 1994). These different redistribution properties have required different clinical protocols from those used for 201T1 SPECT (Berman et al., 1991; Braat et al., 1994; Johnson, 1994; Matsunari et al., 1996a). One adverse aspect of 99mTclabeled myocardial agents is the hepatobiliary excretion of these compounds. Photon scatter from hepatic activity up into the heart can interfere with interpretation and quantitation of inferoapical segments, specially in obese patients or in individuals with high diaphragms (Johnson, 1991). Artifacts in the inferior myocardial wall due to liver and abdominal background have been reported in clinical studies with 99mTcsestamibi (Middleton and Williams, 1993; DePuey, 1994), 99sTcteboroxime (Johnson, 1991; 1994; Chua et al., 1993) and e9mTctetrofosmin (Braat et al., 1994; Matsunari et al., 1996b). In spite of their physical advantages, 99mTc compounds have not replaced 201T1 in clinical studies. Both radionuclides are currently used in cardiac SPECT. SPECT Data Acquisition The first step in SPECT studies is the acquisition of a set of planar images or projections about the human body using one or more rotating gamma cameras. Typical SPECT studies acquire 60 projections at intervals of 6 degrees. The gamma camera orbits the patient in a stepand 15 shoot motion mode. This means that there is no data acquisition during the time in which the detector is moving between two positions. Projections are converted and stored as computer files with formats of 64 x 64, 128 x 128 or 256 x 256 pixels. The most common format in cardiac SPECT is 128 x 128 pixels. Gamma Camera The gamma camera is an image formation device that uses a scintillator detector with a large NaI(Tl) crystal coupled to an array of photomultiplier tubes (Cho et al., 1993, p. 165). The NaI(Tl) crystal can be either circular, with diameter of 40 cm, or rectangular, with dimensions of 40 cm x 20 cm. The thickness of the crystal is approximately 1.25 cm. The number of photomultiplier tubes depends on the manufacturer, but presentday gamma cameras use from 37 to 101 tubes arranged in an hexagonal pattern. A multihole collimator, made of a high atomic number substance, is attached to the external surface of the crystal and confines the direction of the incident photons to an extremely small acceptance solid angle. Thus, the crystal receives gammarays only from sources directly in front of it, making a projection from the 3D source distribution onto the 2D face of the scintillation crystal. When a gamma photon is absorbed within the crystal, the photomultiplier tubes closest to the interaction position receive an amount of light proportional to their distance to the event. The output signals of the photomultiplier tubes are processed by an analog or digital circuitry to calculate a pair of coordinates (x,y) representing the 2D position of the interaction on the crystal. Finally, the output from each tube is summed to obtain a signal proportional to the total energy absorbed within the crystal. A 16 pulseheight analyzer retains only those events whose signal amplitude lies within a selected small window ( 10%) around the amplitude of the photopeak. This energy discrimination excludes a high percent of Compton scattering events that occur in the patient, the collimator or in the crystal itself. Spatial Resolution and Sensitivity Two important variables that characterize the performance of a gamma camera are spatial resolution and sensitivity. Spatial resolution indicates the ability of the system to identify the exact location at which a photon has been emitted. This variable is commonly measured in terms of the full width at halfmaximum (FWHM) of the line spread function (LSF) obtained from the count profile recorded across the image of a line source. The FWHM value critically depends on the energy resolution of the detector, the energy window, the sourcetodetector distance, the geometric response of the collimator and the presence of a scattering medium. The best spatial resolution is nearest the collimator. By increasing the distance from the detector to the source, the spatial resolution is significantly degraded. The FWHM in air at most typical organ depths (= 10 cm) for lowenergy high resolution (LEHR) and low energy general purpose (LEGP) collimators are approximately 7 mm and 9 mm, respectively (Tsui et al., 1994a). In SPECT studies, at the typical distance from the collimator to the center of rotation (= 18 cm), the FWHM is approximately 14.5 cm. The presence of a scattering medium (e.g., the patient) should increase the FWHM up to 16.5 mm. Sensitivity is defined for a point source and expressed in terms of counts per minute per unit of activity per meter from the center of the 17 camera detector (Nichols and Gait, 1995). Unlike spatial resolution, sensitivity decreases only slightly with distance. When distance increases, the number of photons per unit of area decreases by the inverse square of the distance, but then the area of the crystal exposed to the source increases by the square of distance. This is valid only in the absence of an attenuating medium. Conventional parallelhole collimators used in SPECT permit only about 0.015% of the emitted photons to interact with the scintillation crystal. For most SPECT studies, the average count rate is on the order of 104 counts per second. Because of the limited number of counts per pixel, projections can be severely corrupted by random noise. Sensitivity could be improved by increasing the size of collimator holes or by using more holes of smaller diameter but thinner septal thickness. Unfortunately, there is an inverse relationship between spatial resolution and sensitivity. Both solutions to increase sensitivity allow undesirable photons to penetrate the collimator with the subsequent degradation in spatial resolution. The possibility of increasing the number of counts by increasing the acquisition time is not technically viable. Long acquisition times increase the probability of patient involuntary motion which can give rise to image artifacts. On the other hand, the biokinetics of some radiotracers commonly used in cardiac SPECT do not allow long acquisition times. Increasing the amount of radiotracer also is limited by dosimetry considerations and the dead time of SPECT systems (< 1 us). Practical solutions to increase the number of counts in cardiac SPECT while preserving spatial resolution have been proposed by modifying 18 acquisition modes (Galt, 1994; Gait and Germano, 1995). The most common modalities have been the following: a) The introduction of elliptical orbits to minimize the distance from the camera to the body throughout all the acquisition is one modality. b) Acquisition of projections in an orbit of 180 degrees instead of the full 360 degrees rotation is another (this method excludes those projections which are severely distorted by photon attenuation and scatter. The acquisition time is then confined to count acquisitions for the most reliable projections), c) The third modality is rotation of the gamma camera and data acquisition in continuous mode rather than the conventional stepandshoot. In this mode, counts are not wasted during the time the gamma camera is moving. In the stepandshoot mode the waste of counts can be significant because detectors take 2 to 4 seconds to move between two positions (Galt and Germano, 1995). MultiDetector SPECT Systems The most important advance for improving sensitivity has been the development of multidetector SPECT systems: the higher the number of detectors, the larger the number of counts that are acquired in the same period of time. The higher sensitivity of these systems has also improved spatial resolution by allowing the use of high resolution collimators (Galt and Germano, 1995). Two opposed detector systems have been used for many years. They were developed for acquiring simultaneous anterior and posterior projections in body skeleton imaging. Threedetector rotating SPECT systems were introduced by Trionix Research Laboratory (Twinsburg, OH) and Picker International (Cleveland, OH) in the late 1980s (Galt, 1994). Today, almost every manufacturer produces threeheaded SPECT systems which are mostly used for cardiac and brain studies. Other systems 19 designed specifically for cardiac SPECT are those with two rectangular gamma cameras mounted at 90 degrees, for example, the Optima gamma camera from General Electric Medical Systems (Milwaukee, WI), the CardiaL from Elscint (Hackensack, NJ), and others (Galt, 1994; Case et al., 1996; Cullom et al., 1996). SPECT Reconstruction The Inverse Problem in SPECT The 3D reconstruction or inverse problem in SPECT can be stated as follows: Given a set of collimated projections {P#[r]} around a transaxial plane t of constant thickness and along an axis r at right angles to the axis of rotation, calculate the density function of the radioactivity concentration g(x,y) in the transaxial plane t that produces (Pe[r]} (Figure 21). This problem has been extensively investigated using iterative algorithms and Fourier transform based methods. However, currently only the filtered backprojection (FBP) technique has been routinely used in clinical studies and universally incorporated in commercial SPECT systems. The Filtered Backprojection Technique The FBP is a Fourier transform method based on the application of two mathematical operators to projection data. The first one filters projections using a reconstruction bandpass frequency filter. The second operator, which is the backprojector or backprojection operator, assigns each projection filtered pixel values to all pixels along a perpendicular corresponding line through the reconstructed transaxial image plane. This is done for all the pixels in the projection and over all projection 20 r x Pe li1 r =x cos + y sin e Figure 21. Density function g(x,y) in the transaxial plane XY and its projection Po[r] at an angle 8. angles. The reconstruction is achieved by the superposition of the filtered projection values for all angles (Figure 22). Briefly, the method can be explained as follows: For a narrow ray, assuming ideal detector response, no photon attenuation and no photon scatter, the relationship between g(x,y) and the set (Po[r]} (Figure 21) is given by the Radon transform (Tretiak and Metz, 1980): Pg[r] = R[g(x,y)] = fg(x,y) 6 (xcos6 +ysin r) dxdy, where R is the Radon operator. The xy coordinate system is stationary and defines the transaxial image plane. The rs system is rotating and gives 21 0.3 SOriginal Point Profile 0.2 Reconstruction (Ramp) 0.1 t eredProfile 0 0.1 Single Backprojected Profile 0.2 10 5 0 5 10 SWoitionr SBckprojected Profiles from pixel position fib Ang"es (a) (b) Figure 22. Filtered backprojection technique. (a) Convolution kernel or reconstruction filter in the spatial domain; (b) Filtering, backprojection operation and reconstructed transaxial image plane of a central point source. 22 the reference for the projection Pe[r] (Figure 21). Equations for coordinate transformation between the two systems are given by r = xcosO + ysin8 s = xsin8 + ycosO. From the Radon transform and equations for coordinate transformation, it can be demonstrated (Rogers and Clinthorne, 1987) that the onedimensional (1D) Fourier transform (FT) of the projection P,[r] is equal to the FT of g(x,y) along a radial line in the frequency domain at the same angle as the projection. This is the "Central Slice" theorem, which constitutes the basis of Fourier reconstruction methods. Mathematically, it can be expressed as FT[Pe[r]] = Ig(x,y) exp[2j(wxx + wyy) ] dxdy , where wx = w.cos8 and Wy = w.sinS. Then, the above equation can be expressed as FT[Pe[r]] = FT[g(x,y)]e = G(w,O) , where this equation indicates that g(x,y) in the frequency domain G(Wx,wy) can be synthesized from the radial samples G(w,0), so that: g(x,y) = FT1[G(wx, wy)] 23 Expressing Wx and Wy in polar coordinates (w,0), the above equation becomes g(x,y) = ffG(w, ) Iw.exp[2xnjw(x.cosO +y.sin) ] .dwd , where Iwl is the Jacobian of the coordinate transformation. The above integral can be arranged as g(x,y) = Pfp x.coso +y.sin0] d (2.1) where P [r] = FT1[FT[P[r]] .lw ] (2.2) Equations (2.1) and (2.2) indicate that g(x,y) is calculated by applying two operators to the projection data. Equation (2.2) corresponds to a filter operator, where PFo[r] is the filtered projection and Iwl is the ramp filter (Larssen, 1980). Equation (2.1) represents the backprojection operator and indicates the backprojection of PFo[r] onto the plane xy. Because of the divergent nature of Iwi, the integral (2.2) can not be evaluated. Ramachandrian and Lakshminaryanan (1971) solved this problem by considering that FT[P,[r]] is band limited and then Iwl can be truncated at some value IWMIx. Therefore, the ramp filter was modified according to L(w) = IwI if IwI < 2nfN L(w) = 0 Otherwise, 24 where fN is the Nyquist frequency (Gonzalez and Wintz, 1977, p. 70) defined by 1 N 2AX ' and Ax is the onedimensional pixel size. The filtering operation in the frequency domain is equivalent to the convolution of P,[r] with the ramp filter in the spatial domain (Gonzalez and Wintz, 1977, p. 61). The convolution kernel L(r) = FT'[L(w)] is a symmetrical function with negative sidelobes (Figure 22a) that introduces negative or null values in the filtered projection. These negative values compensate the blurring and star effects produced by the backprojection operation (Figure 22b). Smoothing Window Functions There are two practical limitations in the use of the ramp filter. First, L(w) is a highpass band filter. This is an ideal function to reconstruct noisefree data, but not real SPECT projections which are degraded by statistical noise. Second, the truncated ramp filter introduces rim and overshoot artifacts (Gibbs phenomenon) into the image. This is because the hard cutoff of the ramp filter produces ripples in the filter response to image edges (Algazi et al., 1995). Both drawbacks are overcome by multiplying L(w) by a lowpass frequency filter or smoothing window function. Thus, the highpass frequency ramp filter becomes a band pass frequency filter which gradually rolls off the high frequencies rather than cutting off sharply at fN. The most commonly used window functions in clinical SPECT are the Butterworth, Hann, Hamming and Parzen 25 lowpass frequency filters (Larssen, 1980). Figure 23 shows examples of some convolution kernels used in SPECT and their frequency response. The exact form of the convolution kernel depends upon how much lowpass filtering is acceptable in order to reduce the noise in the images while preserving the resolution. The width of the central positive peak determines the amount of smoothing and the negative sidelobes determines the amount that is subtracted from the neighboring backprojected rays. In practical reconstruction algorithms, the window function is a 20 radially symmetric lowpass frequency filter that is applied to projection images before reconstruction. Then, the ramp filter is applied one dimensionally to the transaxial projection data, row by row. The objective of the first operation is to avoid streaking artifacts in sagittal and coronal slices, which are calculated from transaxial data (Larssen, 1980). These artifacts appear when projection data are smoothed in just one direction. The 2D prereconstruction smoothing avoids streaking artifacts but increases blurring in both transaxial and axial directions. Myocardial Tomographic Slices Coronal and sagittal slices are calculated from transaxial slices for further evaluation of the 3D distribution of the radiotracer (Larssen, 1980). They are displayed in three perpendicular axes parallel to the natural vertical long, horizontal long and short axes of the heart (Figure 24). Axes reorientation is performed by two successive image rotations through angles selected by the operator (Borello et al., 1981; Garcia et al., 1985). The most important slices from the clinical point of view, are those in the short axis. These slices are normal to the left ventricular long axis and provide a set of crosssectional slices of the 26 0.1 8 0.08 0.06 6 0.04 0.02 4 0 0.02 0.04 0 0.1 8 0.08 0.06 6 0.04 0.02 4 0 0.02 0.04 0 0.1 0.08 8 0.06 6 0.04 0.02 4 0.02 0.04 0 8 4 0 4 8 0 5 10 15 20 PIXEL POSITION FREQUENCY Figure 23. SPECT reconstruction convolution kernels (left) and their frequency response (right). From top to bottom: Hamming, Hann and Parzen filters. 27 Right RV LV Left Short Axis View Axis View Base Apex (a) Plane for vertical long asi Apex $R / Horizontal Long Axis View le Right Left Plane 11ftor honzontal a long axis Mid Caviy Mid Cavity Short Axis Slice Vercl Long Ai Slie Sep LVL A rl / I2i4 At 415 (b) \2 Septal Lateral Intferior Polar Map Figure 24. Myocardial tomographic slices. (a) Standardized display according to the natural axes of the heart (Berman et al., 1991); (b) Computergenerated bull's eye polar map display of circumferential maximal count profiles of short axis slices. Apex is derived from both short and vertical long axis slices (Matsunari et al., 1996b). 28 cardiac muscle. This reorientation standardizes the cardiac image display and facilitates qualitative and quantitative analysis by exploiting the symmetry of normal left ventricles in the short axis (Berman et al., 1991). Quantitative Analysis Quantitative estimates of the radiotracer distribution are commonly performed using circumferential maximal count profiles of shortaxis slices (Berman et al., 1991). An alternative semiquantitative method is the bull'seye plot or display of SPECT data in the form of a polar map (Garcia et al., 1985). This procedure converts the circumferential maximal count profiles of short axis slices into polar coordinates profiles. Slices are displayed as a series of concentric circles, with the apex at the center, the inferior wall at the bottom, the anterior wall at the top, the septum at the left and the lateral wall at the right (Figure 24). Pixel counts inside all rings are normalized to the maximum pixel counts and expressed in percent of this value. This representation results in a simple 2D display of the 3D data and allows the comparison of patients's bull'seye with reference normal standards (Ficaro et al., 1996; Matsunari et al., 1996b). Degrading Factors in SPECT The FBP method provides accurate 3D reconstruction for idealized projections free of the degrading effects due to the geometric response of the detector, photon attenuation and Compton scatter. When the method is applied to real projections the reconstruction is significantly limited in terms of accurate quantitation, spatial resolution and contrast (Tsui et al., 1994a; 1994b). Therefore, a number of pre and postreconstruction 29 semiquantitative techniques have been proposed for compensating distortions when the FBP technique is used in SPECT reconstructions (Chang, 1978; Larsson, 1980; King et al., 1991; Buvat et al., 1994; Tsui et al., 1994a; Case et al., 1996). The discussion that follows, describes the three main degrading factors and the correction methods mostly used in experimental and clinical cardiac SPECT studies. Geometric Response of the Detector The geometric response of the detector is the main cause of blurring and loss of spatial resolution in SPECT studies (Tsui et al., 1994a). This degrading factor critically depends on the geometric response of the collimator and the intrinsic resolution of the scintillator crystal. However, the major contributor to the detector response is the collimator. The number of collimator holes and their acceptance solid angle are the determining factors in the spatial resolution of a SPECT system. The geometric response of the detector is experimentally characterized by the point spread function (PSF), which fully represents the blurring characteristics of the SPECT system. This is a spatially variant function, which rapidly broadens as the distance from the collimator increases. Photon Attenuation Photon attenuation is the most important source of inaccuracy in SPECT quantitations (Tsui et al., 1994a; 1994b) and is considered the major cause of falsepositive cardiac SPECT studies (King et al., 1996). The attenuation effect can reduce the detected counts to 25 % and 40 % of the unattenuated counts for ""Tc and 201T1, respectively (Tsui et al, 1994a). Because photon attenuation is a function of the thickness and composition of the absorbing medium, the attenuation effect is not uniform 30 for an imaged organ. This is specially important in cardiac SPECT due to the geometric position of the heart and because surrounding tissues (lung, muscle, adipose tissue, bone) have very different attenuation coefficients. When cardiac SPECT is performed in supine position, there is a significant decrease in myocardial inferior wall counts as a result of attenuation by the left hemidiaphragm (Segall and Davis, 1989; Wallis et al., 1995) and the overlaying right ventricle and right ventricular blood pool (DePuey and Garcia, 1989). However, the problem is not limited to the reduction of counts in the inferior wall. Significant artifacts are induced in the anterior wall due to attenuation from the breast tissue in females (DePuey and Garcia, 1989). Also, in obese patients the accumulation of adipose tissue in the lateral chest wall may result in an apparent perfusion defect in the lateral wall (Wallis et al., 1995). Attenuation Correction. The attenuation problem in SPECT can be mathematically described by the inclusion of a negative exponential term in the Radon transform (Tretiak and Metz, 1980). The resultant projection is known as the attenuated Radon transform (Gullberg and Budinger, 1981). So far, an analytical solution to the attenuated Radon transform has not been found. Different attenuation correction numerical algorithms have been proposed, but a complete solution to the problem has not as yet been proposed (Gullberg and Budinger, 1981; Tanaka et al., 1984; Tsui et al., 1994a; King et al., 1996). However, considerable progress has been achieved in the last decade due to improvements in correction algorithms, hardware of SPECT systems and increased computing power available in current SPECT computers. 31 Depending on their application before, during or after reconstruction, there are three broad approaches for attenuation correction: a) preprocessing correction; b) intrinsic correction; and c) postprocessing correction. Preprocessing methods are based on the assumptions of a uniform concentration of activity and a constant attenuation coefficient in the absorbing medium. One method corrects by hyperbolic functions the arithmetic or geometric mean of conjugate opposing projections (P,[r] and P,+,[r]) before filtering and backprojection (Keyes, 1976; Larssen, 1980; Tsui et al, 1994a). A practical inconvenience is that the body contour is needed. The simplest and easiest prereconstruction correction makes the reconstruction for 180" with the arithmetic or geometric mean of the conjugate opposing projections which were taken over 360*. This correction is employed in commercially available SPECT systems. Intrinsic compensation methods include attenuation correction directly into the reconstruction process. They can be divided into two categories: a) analytical and b) numerical. The first group is based on an approximate solution of the attenuated Radon transform with uniform attenuation coefficient (Gullberg and Budinger, 1981; Tanaka et al, 1984). Analytical intrinsic methods use a first correction for the body outline. This is accomplished by multiplying each projection pixel by the inverse of the attenuation along a ray from the projection pixel to the center of rotation. After that, a second correction compensates for the actual attenuation distance. This can be performed in various ways, but the best known method is the Gullberg and Budinger (1981) attenuationweighted backprojection technique. This method uses window functions which depend 32 on a uniform attenuation coefficient (Gullberg and Budinger, 1981). Filtering is performed with a modified ramp filter in which its value is zero in the frequency range below p/2w (Gullberg and Budinger, 1981). Then, filtered projections are backprojected and multiplied by an attenuation weighting factor which depends on the position of the pixel in the reconstruction matrix. One important limitation is that these algorithms amplify noise and introduce artifacts in lowcount images. Bellini et al. (1979) have proposed an analytical solution for the attenuated Radon transform, which relates the nonattenuated and attenuated projection sinograms in the frequency domain. The sinogram is an image which represents the projection data as a function of projection angles. The sinogram is formed by taken the same row of data from all projection images and creating a new image in which the xcoordinate values are those of the projection and the yvalues are determined by the projection angle. The image is called "the sinogram" because a single point should follow a sine curve in this type of representation. Thus, the relationship described by Bellini is given by P'(v,0) = P(v,8). [(v2+.2)1/2,0B+isinh1(p/v)] where P' and P represent the attenuationcompensated and original frequency domain sinograms, respectively. The Bellini's method requires the interpolation of the sinogram in the frequency domain to evaluate P(v,O) at the frequency value (v2 + 2)1/2 and at complex angle 0 + i.sinh1'(v/p). After interpolation in the v and 0 directions, the Fourier inversetransform of the sinogram is performed to obtain the attenuationcompensated sinogram. The advantages of this 33 method are its speed and that can be used with the FBP technique. Disadvantages are that corrections are only valid for convex regions of uniform attenuation around the center of rotation and ignores the problem of distancedependent resolution. In addition, the Bellini's method has not been tested extensively as other methods and its robustness in clinical studies is unknown (King et al., 1995). Intrinsic numerical methods incorporate a nonuniform attenuation coefficient map into an iterative statistical estimation algorithm such as the maximum likelihoodexpectation maximization (MLEM) or the weighted least squaresconjugate gradient (WLSCG) (Tsui et al., 1994b). The MLEM method assumes a Poisson distribution in the projection counts and determines the source distribution that most likely reproduces the projection data (Shepp and Vardi, 1982). Briefly, the MLEM method works as follows. An initial estimate of the transaxial source distribution is performed, e.g., using the FBP technique. Next, each transaxial pixel is updated according to the iteration scheme given by g (x,y) =g1(x,y) .B[ Pl] , P81 [r] where B is the backprojection operator and the index i enumerates the iteration. The previous transaxial source distribution g'l(x,y) is corrected by the backprojection of the ratio of the actual projection counts P,[r] over the projection counts Pg'[r] which were calculated assuming a source distribution given by g'(x,y). In this scheme, the backprojection operator B weights the projection counts by the likelihood that they actually contribute to the source distribution. The 3D 34 implementation of this iterative scheme is straightforward. In theory, the major advantage of this method is that the physical image formation process (attenuation, geometric collimator response and even scatter) can be modeled into the projector and backprojector operators. The inclusion of degrading effects into the reconstruction scheme allows more accurate corrections and source distribution quantitation (Tsui et al., 1994a; 1994b). This method have two major disadvantages. The first is the processing time which can be excessive for clinical studies, specially if 3D collimator geometric response corrections are included. The second is the slow and image dependent convergence, in which the noise increases as the number of iterations increases. Other iterative reconstruction algorithms (Tsui et al., 1995) have been proposed to speed the convergence and reduces noise, but so far the MLEM algorithm is the most widely used and the standard for accurate quantitation comparations with other reconstruction methods. When only attenuation correction is included in the MLEM reconstruction algorithm, proper compensation may be achieved after a limited number of iterations. Recently, several commercial system (Vantage system from ADAC Laboratories, Optima Lshaped Cardiac Camera from General Electric Medical Systems, and PRISM 3000 SPECT Camera from Picker Inc.) have included the MLEM reconstruction algorithm in conjunction with transaxial attenuation maps for performing accurate nonuniform attenuation corrections. Attenuation coefficient maps have been simultaneously determined using point or line sources in multidetector SPECT systems. At the present time, these systems are in process of clinical validation (Case et al., 1996; Cullom et al., 1996). 35 Chang Attenuation Correction Algorithm Postprocessing techniques are based on the original Chang attenuation correction algorithm (Chang, 1978). This is the best known attenuation correction procedure and was specifically designed to work in conjunction with the FBP method. Most commercial SPECT systems have incorporated this correction technique. The algorithm is implemented by dividing each reconstructed pixel value by the average attenuation factor for that pixel (Chang, 1978). The average attenuation factor is calculated as the average attenuation along all rays from the pixel to the boundary of the attenuating medium. Therefore, the correction factor for attenuation at a point (x,y) is given by c(x, y) =  exp (P. lk) k1 where M is the total number of projections and 1k is the distance between point (x,y) and the boundary point of the medium at projection k. When the source is extensively distributed, the above correction factor over or undercorrects some parts of the image, depending on the source distribution. Next, a second correction is needed. This correction is performed by projecting the corrected data to form a new set of projections. A set of error projections is obtained by subtracting each corrected projection from its corresponding original projection. Error transaxial slices are reconstructed by using the FBP method and corrected for attenuation again. These slices are added to the initially corrected slices to form the final attenuationcorrected image. This correction can 36 be repeated using more than one iterative step. Experimental studies have demonstrated that the Chang algorithm gives the best results after only one or two iterations (Tsui et al., 1989). For nonuniform attenuation, the attenuation coefficient distribution can be calculated from an attenuation map (Tsui et al., 1989; Manglos et al., 1987). This method provides a significant improvement in image quality and acceptable quantitation for practical purposes (Manglos et al., 1987; Tsui et al., 1989; 1994b). Photon Scattering Compton scattering is the dominant interaction in tissue for the energy range (50 keV to 300 keV) used in nuclear medicine (Evans,1985, p. 714). In addition, small scattering angles are more probable than large deflections for the emission energies of 201T1 and 99mTc. Hence, a relatively high number of forward deflections can pass through the collimator holes. In theory, scattered photons should be separated from those unscattered using the pulse height analyzer of the detector. However, because the finite energy resolution and small sensitivity of SPECT detectors, a certain number of scattered photons has to be counted in order to detect the largest possible number of unscattered events (O'Connor et al., 1995a; Welch et al., 1995). Photon scatter accounts for approximately 40 % and 60 % of the total counts for 99mTc and 201T1 cardiac SPECT studies, respectively (O'Connor et al., 1995a; 1995b). This undesirable effect degrades the reconstructed slices by blurring fine details and lowering contrast, and leads to inaccuracies in the quantitative estimate of perfusion defect size (O'Connor et al., 1995a; 1995b; Case et al., 1996; Cullom et al., 1996; McCartney et al., 1996). 37 For projections acquired with parallel hole collimators, photons emitted in a given transaxial plane can be detected as emitted in other transaxial plane due to Compton scattering and the finite detector resolution. This effect could be significant in cardiac SPECT where the outofplane scattered liver photons are detected in the inferior myocardial planes (Smith, 1994; King et al., 1995). A simple scatter model based on the KleinNishina (Evans, 1982, p. 683) differential probability for a scatter angle and incident energy can give an estimate of the contribution and importance of outofplane scattered photons on a given transaxial plane. It is known (Evans, 1982, p. 675) that the relation between the initial energy E. and final energy E, of a photon experiencing Compton interaction is given by E = E 1 +a. (1cos) ' where q is the angle between the initial and final direction of the photon and a is the ratio of the incident photon energy and the electron rest mass. Assuming that the detector has perfect energy resolution and that a rectangular energy window AE is centered at Eo, the maximum allowable first order scatter angle within this energy window is calculated by P.x = cos1'[1 A] , a. (2Eo, AE where for 99"Tc: E. = 140 keV, a = 0.274 and pMAx = 54'. In this simple model, a photon is allowed to be scattered only once and photon paths are only traced from each voxel containing activity towards the detection plane (Figure 25). The incident photon flux 38 b Figure 25. Detection geometry of unscattered and outofplane scattered photons. It is assumed a perfect geometry resolution and only perpendicular photons to the plane of detection are registered. The outofplane voxel is represented by "a". Voxel "b" belongs to the emission plane. Photons emitted by "a" and scattered at voxel "b" becomes an additional effective source of activity at b (Smith et al., 1993). density I (photons/cm2voxel) emitted by an outofplane voxel "a" which can reach a voxel "b" in the emission transaxial plane (Figure 25) can be expressed as @_ = .exp(p.d) , where A, is the activity of the outofplane voxel, d (cm) is the distance from the outofplane voxel "a" to the emission voxel "b" and p (cm"') is the narrowbeam attenuation coefficient of the uniform medium. It is supposed that the flux arriving to the emission voxel is scattered and becomes in an additional source of activity. From the simple geometry of Figure 25, oncescattered photons will be detected only when the 39 scattered angle is given by: W = sin'(At/d) : 54", where At is the distance between the outofplane and emission plane. Assuming a pixel size of 0.356 cm, At = 0.356 x n, where n is the number of transaxial planes between the voxel "a" and the emission voxel "b" (Figure 25). Then, the above equation can be represented by A,, sin2, S= A. sn2 exp(0.0534n/sin/) , 1.593n2 Notice that the scatter angle has been chosen such that scattered photons propagate only along a perpendicular projection ray accepted by the detector (Figure 25). Therefore, the total number of outofplane photons which are scattered by the voxel "b" and registered by the detector, is calculated by adding the contributions of all outofplane source voxels multiplied by the probability of scattering at the corresponding angle ]((i).( AY) .At,. where Aarv/A4i is the differential cross section for photon scattering at angle Wi per unit of scattering angle and per voxel. This is calculated from (Evans, 1982, p. 690) .av) AlO l AQ ) ( A *p Am. ( > where, p, = 3.343 x 1023 electrons/g is the electron density of water (Attix 1986, p. 531), mV = 1 g/cm3 x (0.356 cm)3 is the mass of a voxel, 40 AD/Ap, = 2xsiny, is the differential solid angle per unit of scattering angle (Evans, 1982, p. 690). Finally, a.o/AQ is the differential scattering cross section at angle per unit solid angle and per electron and is calculated from (Evans, 1982, p. 683) A0_ r (1+cos22) (1+C (1cos) ) +a2 (1cos) 2 A 2 2 (1+a (1cosO)) where ro = 2.8179 x 1013 cm is the classical electron radius (Attix, 1986, p. 525). Figure 26 shows the relative contribution of outofplane scattered photons per voxel as a function of the scattering angle for different out ofplane distances. According to this model, the relative contribution of nearest planes to the outofplane detected scattered photons at the voxel "b" are 73%, 16.6%, 6.63% and 3.53% for the first, second, third and fourth neighboring planes, respectively. Table 2.1 illustrates the contribution of outofplane scattered photons per voxel as a function of the activity ratio of outofplane and emission voxel. The number of photons actually emitted by the emission voxel (unscattered photons) and registered by the detector was calculated as: A,/4w, where Ab is the activity of the emission voxel "b" (Figure 25). When the emission voxel has the same activity as outofplane neighboring voxels, the total contribution of scattered photons can be 20.1%. If the activity of outofplane voxels increase, then their contribution can reach more than 50% (Table 5.1). Notice that this model does not include attenuation along the path between the emission voxel "b" (Figure 25) and the detector. This is because it is assumed that attenuation along this 41 120 S80 2 60 E 40 20 ~4O 75 0 0 10 20 30 40 50 60 Angle of scattered photon (degrees) Figure 26. Plot of the relative number of outofplane scattered photons per voxel as a function of scattering angle. The contributions of the first, second, third and fourth neighbouring planes are represented from the highest to lowest curves, respectively. 42 Table 2.1 Contribution of outofplane scattered photons per voxel as a function of the activity ratio of outofplane to emission voxels (ROE) for different neighboring planes. Distance between planes: 0.356 cm. SROE = 0.5 ROE = 1.0 ROE = 1.5 ROE = 2.5 ROE = 3.5 First 7.2% 14.3% 21.5% 35.8% 50.1% Plane Second 1.8% 3.6% 5.4% 9.0% 12.6% plane Third 0.75% 1.5% 2.3% 3.8% 5.3% plane Fourth 0.35% 0.71% 1.1% 1.8% 2.5% plane Total 10.1% 20.1% 30.3% 50.4% 70.5% path is corrected by the reconstruction process. However, photon attenuation along the path between voxels "a" and "b" (Figure 25) is included, because correction methods in SPECT only compensate attenuation along those projection rays which are perpendicular to the detector. In this model, it was assumed that all outofplane neighboring voxels have the same activity A,. It is obvious that this assumption increases the number of scattered photons in the calculations. However, this approximation could be valid in those situations in which high activity is taken by large organs nearby to the region or organ to be imaged (e.g., the high liver uptake of sg9Tclabeled radiotracers that is observed in cardiac SPECT studies). Data derived (Table 2.1) from this simple scatter model confirm the significant contribution of outofplane A 43 scattered photons and illustrate the need of scatter compensation methods for improving accuracy in SPECT. Scatter Correction Methods The main benefit of any scatter correction method is to improve quantitation accuracy by removing those events with inaccurate positional information. Over the last several years a number of technological advances in instrumentation have occurred that reduce scatter (O'Connor et al., 1995a; Case et al., 1996). Major advances have been the improved detector energy resolution, with some systems achieving values of 8% to 10% at 140 KeV (O'Connor et al., 1995a) and the implementation of scatter correction hardware/software on some of the newer SPECT systems (Buvat et al., 1994; Case et al., 1996; Cullom et al., 1996; McCartney et al., 1996). Scattering correction techniques can be broadly considered into two groups: energy window based methods and restoration based corrections. Energy window based methods estimate the scatter component by using the data of one or more energy windows abutted to the photopeak window. This approach to scatter correction uses the data acquired with energy windows located on the Compton portion of the spectrum to estimate the scatter contribution to the photopeak region. There is a number of procedures for acquiring and processing energy spectral information, but in general the compensation consists of subtracting the estimated scatter component from the photopeak data (Jaszczak et al., 1985; Hademenos et al., 1993). The method was originally proposed by Jaszczak et al. (1985) who reconstructed transaxial slices from a second window on the Comptom region (92 to 125 keV for egmTc). Next, a fraction of the Compton image (Is) is subtracted 44 from the photopeak image (Ipp) to obtain the corrected image (Ic) Ic = Ipp kIs where the factor k is empirically determined for each system and radionuclide. Although the subtraction process could increase the noise, these methods have demonstrated good clinical performance and have been incorporated into many commercial SPECT systems (Case et al., 1996; Cullom et al., 1996; McCartney et al., 1996). Restorationbased methods attempt to return scattered photons to their original emission sites. They can be divided into two categories: scatter operator based corrections and filtering restoration. The first category is based on modeling the scattering process by a mathematical operator and including it into an iterative reconstruction algorithm (Buvat et al., 1994). Since Compton scattering is a 3D effect that depends on the position within the medium and the source distribution, exact scatter response functions can be only derived from Monte Carlo calculations (Frey and Tsui, 1991; Buvat et al., 1994). The large memory and extensive computation required for Monte Carlo simulations on each patient make these methods impractical for clinical use. Filtering restoration methods are more easily implemented and do not require extra memory or time for acquiring and processing SPECT data. They compensate the degrading effects of both scatter and detector geometric response, while also decreasing noise (King et al., 1991; Boulfelfel et al., 1992). These methods are based on the deconvolution of SPECT data with a representative point source response function of the imaging system. 45 Because the illconditioned nature of the restoration problem (Andrew & Hunt, 1977, p. 114), a number of different regularization approaches has been developed (Madsen and Park, 1985; King et al., 1991; Boulfelfel et al., 1992). So far, the simplest and most effective methods in SPECT are those that filter the data using the Metz or Wiener filters (King et al., 1991; Boulfelfel et al., 1992). Restoration based methods are discussed in detail in Chapter 3. Artifacts in the Inferior Myocardial Wall Two main types of artifacts have been observed in the inferior myocardial wall as consequence of the photon attenuation, scatter, geometric collimator response and high background in abdominal organs. First, the high liver activity in 99mTc cardiac SPECT studies produces a hypoperfusion artifact in some patients (Johnson et al., 1991; 1994; Matsunari et al., 1996b). The origin of this artifact has been associated with the photon attenuation (King et al., 1995; Nuyts et al., 1995) and the spatial domain negative sidelobes of reconstruction filters (Figure 2 2) (Ljunjberg et al., 1994; Zeng et al., 1995). On the other hand, experimental studies (Ljunjberg et al., 1994; King et al., 1995) have demonstrated that when photon attenuation is compensated, the hypoperfusion artifact changes to an overcounting artifact in the inferior myocardial wall. This second type of artifact has also been observed in 201T1 studies with high liver uptake of the radiotracer (Case et al., 1996; McCartney et al., 1996). 46 The Hypoperfusion Nyocardial Artifact This type of artifact is typical of the FBP method when the structure to image is close to a very high activity region. A similar artifact has been observed in SPECT studies of the skeleton where high bladder activity produces very low counts in the hips (Bunker et al., 1990). This is due to the convolution of the very high pixel values with the negative sidelobes of the convolution kernel (Figure 22). In theory, a convolution kernel with very narrow sidelobes could solve the problem (Zeng et al., 1995). However, when the negative sidelobes are reduced, blurring is increased. The hypoperfusion artifact can then be eliminated, but an artificial increase of the counts will appear as a consequence of the blurring of the image. Taking a different point of view, Nuyts et al. (1995) based the origin of the hypoperfusion artifact on the fact that attenuated projection data are inconsistent. According to Chang (1978), inconsistent projection data means that there does not exist any realistic object able to produce such a set of projections in the absence of attenuation. Nuyts el at. (1995) demonstrated that the FBP method does not introduce significant negative values for consistent projections, because the backprojection operator compensates the negative convolved values. They concluded that the main cause of the artifact is the lack of attenuation correction. The Overcounting or Spillover Artifact The conclusion of Nuyts et al (1995) have been confirmed by King et al. (1995) in a 3D simulated antrophomorphic phantom. However, when the projections are corrected for attenuation, a spurious increase or 47 spillover effect in the number of counts appears in the inferior myocardial wall (Ljunjberg et al., 1994; King et al., 1995). This effect has been mostly attributed to the scattered photons from radioactivity in nearby organs, such as liver and stomach (Case et al., 1996; Maniawski et al., 1996; McCartney et al., 1996). It is known that in cardiac SPECT the most important source of scattered photons is the liver, while the most affected region is the inferior myocardial wall. Also, the position of the lungs make the scatter effect higher around the inferior myocardial wall than in superior regions (McCartney et al., 1996). These observations have indicated that optimal and accurate quantitation of cardiac SPECT requires scatter compensation in addition to attenuation correction. Furthermore, the practical availability of transmissionemission methodologies for accurate nonhomogeneous attenuation correction (Case et al., 1996; Maniawski et al., 1996; McCartney et al., 1996), has signified the overcounting or spillover artifact and renewed the interest for scatter compensation techniques in cardiac SPECT. CHAPTER 3 THE RESTORATION PROBLEM IN SPECT Image formation is the process in which an energy source distribution s(x',y',z') is mapped onto a 2D function f(x,y) which is proportional to the intensity of a physical magnitude, such as density, activity or brightness. The source distribution is called the object or input image, while the resultant function is referred as the degraded or output image. Any image formation system imposes some degrading effects upon the output image, e.g., blurring, sampling artifacts or noise addition. Image restoration attempts to recover the input image from the degraded image by using a priori knowledge of degrading effects. Therefore, restoration methods are based on modeling the image formation system by a degradation operator H which operates on the object s(x',y',z') to produce the output f(x,y) f(x,y) =H[s(x',y',z')] (3.1) where the input image should be recovered by applying an inverse operator. In deterministic models, the degradation operator H is represented by the point spread function (PSF) of the imaging system. This function describes the response of the system to a point source in the object of interest. The PSF usually is a spatially shiftvariant function. This means that the function depends on the position of the point source inside the object. When the source is embedded in a scattering medium, the PSF is 48 49 also object dependent (Andrews & Hunt, 1977, p.62). Hence, the PSF can be written as: (3.2) PSPF=h(x, z Il,,x,y,s(xt,yt,zl)) (3.2) Since the resultant image is the superposition of all the individual point source responses, Equation (3.1) becomes: f (x, y) = fffh(x, y', z',x y, s (x', y', z') ) dx'dy'dz. This integral equation represents a nonlinear model with respect to the source distribution. There is no analytical solution for this type of model and it must be approximated to a linear representation in which the inversion process can be performed by using concepts of linear systems theory. This chapter first introduces the restoration problem in SPECT in the most general and realistic way by using a nonlinear spatially shift variant degradation model. Next, the assumptions required to approximate this model to a linear spatially shiftvariant or invariant representation are described and justified. Also, different practical and experimental approaches to the restoration problem in SPECT are discussed. Finally, the design of the restoration filter and the experimental determination of the representative PSF of the SPECT system used in this work are described. Degradation Models in SPECT Nonlinear ShiftVariant Model In contrast to attenuation that is independent of the source distribution, scattering is affected by the spatial density of the source 50 (Buvat et al., 1994). No pointwise scattering correction exists because the scatter at one position (x',y',z') depends upon the source distribution at all other positions. The scatter correction has to involve an integral operator rather than a multiplicative operation as those used in attenuation correction (e.g.: Chang's method). Therefore, a realistic degradation model in SPECT has to be nonlinear with respect to the source distribution. In this model, projections p,(x,y) at an angle 0 can be represented by PO(x, y) = ffh(x', y, z',x,y, s) (x',y, z' dx/dy'dz + n (x,y) where h,(x',y',z',x,y,s) is the source dependent PSF of the system at the projection angle 0 and n,(x,y) is a Poisson noise component added to the projection data to model the statistical nature of the radioactive disintegration. This integral equation represents a continuouscontinuous model because both the projection and the source distribution are continuous functions. However, projections and source distributions in SPECT are sampled and discretized in pixels and voxels, respectively. For practical purposes the continuouscontinuous model has to be converted into a discretediscrete representation XI Y' ZI x 1 y 1 '1 where: x = 1,2,...., X and y = 1, 2,..., Y, and the number of projection pixels in axis x and y, are X and Y, respectively. The number of pixels of the source distribution in axis x', y' and z' areX', Y' and Z', respectively. 51 To simplify the notation, the voxel (x',y',z') will be represented by a voxel index r' and the projection pixel (x,y) by a projection index t. Then, the discretediscrete representation for each projection angle becomes Ps.h' e h..tst .s,, + ng, 1 1 and the description of the SPECT problem is given by a system of M equations of this type, where M is the number of projection angles. The complete system of equations can be expressed in matricial notation as P= H(S) .S+N, where P is the vector of all projection data measurements with MXY elements. H(S) is a huge (MXYX'Y'Z' elements), nonsparse transfer matrix in which each element is proportional to the probability that a photon emitted by a source voxel r' will be detected by projection pixel t. Vectors of source voxel data and noise component are S and N, respectively. The transfer matrix H(S) contains a 3D full description of the PSFs of the system at all projection angles and includes the effects of attenuation, scattering, depthdependent collimator response and irregular body contour. The solution of the above equation allows simultaneous SPECT reconstruction and corrections for attenuation, photon scatter and collimator geometric response. Because each element of the transfer matrix depends on the source distribution, the well known methods of linear algebra can not be applied to calculate the vector of source voxel data. 52 In this case, the inverse solution is a very difficult task that must be incorporated into a complex iterative reconstruction scheme. In addition, the transfer matrix consumes an excessive amount of computer memory: 64 GB is required in a typical SPECT study with 60 projections of 64 x 64 pixels, source distribution sampled by 64 x 64 x 32 voxels and 2 bytes per integer word. On the other hand, modeling the transfer matrix H(S) is a very difficult problem which does not have an analytical solution. Despite considerable effort by several investigators (Buvat et al., 1994), a complete and practical solution has not as yet been obtained. The only theoretical possible approach is to simulate by Monte Carlo the transfer matrix for each projection. This requires an extraordinary amount of computer time, that has been estimated at 190 days with the present SPECT computer technology (Frey & Tsui, 1994). Because PSFs are object dependent, the simulation must be repeated for each patient. The impracticability of the nonlinear model has lead to some approximate assumptions in which the PSF is object independent and spatially shiftvariant or invariant (King et al., 1991; Frey & Tsui, 1990; Smith, 1994). These more conventional representations allow one to solve the restoration problem using linear algebra and iterative statistical estimation techniques. Table 31 shows the main three different approaches to the restoration problem in SPECT. Linear ShiftVariant Model This model assumes that scattering is independent of the source distribution. The physical meaning of this assumption is that the system response to any arbitrary source distribution is equivalent to the linear combination of the response to individual point sources. This 53 TABLE 3.1 Restoration models in SPECT and their approximate solutions. Model Assumptions Solution Nonlinear and PSF is object and There is no practical spatially position dependent. solution. shift variant. Linear and The PSF is object The PSF is parameterized spatially independent, and incorporated into an shift iterative reconstruction variant. Source is uniformly algorithm. distributed or can be approximated by The model has only been a slowly varying used in simulated and function. physical phantoms. Linear and The imaging system Pre or post reconstruction spatially is described by a filtering restoration with shift representative PSF. the FBP technique. invariant. Source is uniformly Applicable in clinical distributed in a studies. uniform attenuating medium. PSF: Point spread function; FBP: Filtered backprojection 54 approximation is justified when the source is uniformly distributed, or at least it is a very slowly spatial variant function (Andrews & Hunt, 1977, p. 63) in an isotropic homogeneous medium. The model can be represented by a matrix formulation of the projection of the source onto detection planes (Frey et al., 1993; Frey and Tsui, 1994; Smith, 1994) or in terms of the photon transport theory (Egbert and May, 1980). Matricial Formulation. In this model the transfer matrix H does not depend on the source distribution. The nonlinear problem has been solved but not the large memory required to store the degrading operator in matricial form. Also, it is necessary to know the PSF for every point and for every projection angle. PSFs have been calculated using Monte Carlo techniques (Floyd et al., 1985) and the results stored in computer memory for use during the reconstruction process. Due to computational limitations, this method have been only applied to very simple geometries and source distributions. This method is impractical for the complex geometries and source distributions found in clinical studies. To reduce the computational burden, Frey and Tsui (1990; 1991) introduced the concept of a parameterized scatter response function. Then, the elements of the transfer matrix can be calculated from analytical functions. These functions were determined for different source locations by fitting Monte Carlo simulated data or experimental measurements to empirical functions composed of a Gaussian plus a pair of half Gaussians with different widths (Frey & Tsui, 1990; 1991; Beckman et al., 1994). They demonstrated that the fitting parameters vary in a systematic manner with the source position in a cylindrical water filled phantom. They also developed a functional description of the scatter response function by parameterizing 55 the Gaussian functions (Frey & Tsui, 1991). The addition of the parameterized scatter response function to a projector/backprojector iterative reconstruction scheme has demonstrated superior absolute quantitative accuracy compared to compensation for only uniform attenuation and geometric collimator response in simulated cylindrical phantoms (Frey & Tsui, 1993; Frey et al., 1993). The concept of a parameterized scatter response function has been recently extended to nonuniform attenuators by using water equivalent distances instead of geometric distances (Frey & Tsui, 1994). The use of parameterized scatter response functions eliminates the need for Monte Carlo simulations to determine the transfer matrix and large memory required to store it. One inconvenience of this method is that the parameterized functions have to be determined for each SPECT system, collimator and photon energy. An additional drawback is that direct solution of the model by matrix inversion is not practical due to to the large number of matrix elements. The solution requires a projector/backprojector iterative reconstruction technique with its inherent limitations such as the dependence of accuracy on the number of iterations and its lack of application in the clinical routine due to the slow convergence and relatively long processing time. Smith (1994) proposed an alternative procedure in which scattered photons were compensated using the dual window energy subtraction technique (Chapter 2, p. 44). In this approach, the 3D problem is approximated by series of coupled 2D transaxial reconstructions. The source distribution in each transaxial plane is calculated iteratively using a small transfer matrix which only contains the contribution of the 56 unscattered photons from a small number of nearby planes. Therefore, the size of the transfer matrix and memory requirements are significantly reduced. In each iterative step, the source distribution is estimated by subtracting the contribution of the outofplane unscattered photons from the projection data (Smith, 1994). The increase in noise level due to subtraction corrections makes this method unapplicable to the count limited clinical SPECT data (Smith, 1994). Photon Transport Theory. In terms of the photon transport theory the problem can be represented as follows: let us assume a stationary process given by a monoenergetic gamma ray isotropic source distributed in a homogeneous and isotropic medium of volume V (Figure 31). Then, the integral transport equation can be expressed as the one speed approximation (Duderstadt and Hamilton, 1976, p. 128) given by (r) = fS(r') K(r, r) d / + L (r .K(r, r) .d3r , where O(r) is the particle flux (photons/cm2.s) at r, s (cm') is the Compton attenuation coefficient (Attix, 1986, p. 132) of the medium and K(r,r') is the kernel of the integral which represents the particle flux at r due to a unit strenght, monoenergetic, steady and isotropic point source at r'. The above integral equation is a homogeneous Fredholm linear equation of the second kind whose solution is given by the Neumann expansion (Zwillinger, 1996, p. 430) (r)= f r(r') .K(r,r') .d3r + P k(rr') .S(r') .dS'r, 57 X z' A x D Figure 31. Coordinates for source distribution and detection plane D. Shaded area represents the volume V in which the source is distributed. where K(rrf) =fKn1(r,r/).K(r/,r").d3zr" for n = 2,3,4,... and where KI (r, ') =K(.r, ) . The physical meaning of each nth term is the contribution of the n times collided photons to the total flux. The first term ( n = 0) represents the uncollided flux, the second one (n = 1) represents the oncescattered photons and so forth. In SPECT studies the narrow energy window around the photopeak ( 10 %) and the small acceptance angle of the 58 collimator holes exclude significantly the contribution of more than once scattered photons. It has been calculated that oncescattered photons constitute more than 90 % of the scattered detected gammarays in ""Tc SPECT studies (Welch et al., 1995). Hence, in the analysis that follows, only first order collisions (n = 1) will be considered. Then, the solution of the integral transport equation becomes (r)C =S (r') .K(r, r ') d&r 1+p vK(r, z") .K(r", r') .S(r) d 3r1d3r11, where the projection p(x,y) of the source distribution onto the detection plane D (Figure 31) is obtained from the line integral of the flux along a projection ray normal to the detection plane. If photon attenuation is included, the line integral is the attenuated Radon transform of the flux (Tretiak and Metz, 1980). Expressing the above equation and the attenuated Radon transform in operator notation p(x,y)= RU[S] +RC[S] , where R is the attenuated Radon transform operator, U and C are integral operators representing the uncollided and once collided flux, respectively. An approximative solution S' can be proposed as S'= GRU[S] , where G is a pseudoinverse operator which represents the reconstruction process and includes attenuation correction, the Radon inverse transform and collimator geometric response correction. Then, by applying G to the 59 projection p(x,y) S'=G[p (x, y) ] GRU[ S where the first term is the source distribution estimated without scatter correction, while the second term contains the exact source distribution that is unknown. The solution can be approximated by an iterative scheme such that S=0 , S=G [p(x,y)] S2= G[p(x,y)] GRC[S1] S = G[p(x,y)] GRC[SI1] , where S1 is the reconstruction obtained with the FBP method and corrected for attenuation (e.g., Chang's method), S2 includes a first order correction for oncescattered photons and so forth. Now, the problem is how to estimate the combined operator GRC and the large computer memory required to store it as an array of elements for each projection angle 0. Any practical solution implies an oversimplification of the model, such as the operator representation by a biexponential empirical function (Egbert and May, 1980). This simple model has given good results but only for a small source in the center of a simulated cylindrical water tank. 60 Linear ShiftInvariant Model In order to avoid iterative methods and adapt restoration to the clinical setting by using the FBP technique, the PSF has been considered as a 2D shiftinvariant function (Boulfelfel et al., 1992; King et al., 1991). Then, if the axis z' is normal to the detection plane D (Figure 3 1), the discretediscrete representation for each projection angle becomes: X' Y' Z' Notice that in this model, he,x',y',x,y is depth independent. Therefore, the above equation can be expressed as x' T' Pxy = .xy.x.y* Qx' + . x1 y 1 where Q,.y. is the projection of the source distribution onto the plane x',y'. This is expressed by: Z/ X1 Now the PSF has a different meaning from those in the previous models. Rather than an operator representing a full 3D physical process, the PSF is a mathematical linear operator that maps a 2D object projection onto a degraded projection. In this model restoration must recover the undegraded or object projection Q,,,, from the degraded projection px, for each projection angle 0. After that, the source 61 distribution can be calculated by applying the FBP method (Chapter 2, p. 20) to restored projections Q,,x.,.. Because the attenuation effect is not included in the PSF, attenuation corrections must be performed before or after reconstruction (Chapter 2, p. 30). The condition of spatially shiftinvariance of the PSF along both directions x' and y', results in the following property (Andrews and Hunt, 1977, p.70) 1)6.xl.y.x.y = h.xxr.yy ' where the physical meaning of this relation is that the PSF response at any point in the object depends only on the value at this point and not on the position of the point (Gonzalez and Wintz, 1977, p. 185). Then, the distorted projection becomes: x IY A'.,y = ]P :P xX/,yy/ QS ,J ,y1 +J"Jo.,y . It can be assumed that projections are discrete periodic functions with periods X' and Y' in the axis x' and y', respectively. Then, the double summation in the above equation is the expression of the 2D discrete convolution of the PSF with the object projection (Andrews and Hunt, 1977, p. 66). The advantage of this representation is that in the Fourier frequency domain, the convolution operation becomes a simple multiplicative operation (Gonzalez and Wintz, 1977, p. 65). By applying the discrete Fourier transform (DFT) and the convolution theorem (Gonzalez and Wintz, 1977, p. 65), the above equation can be expressed in the 62 frequency domain as P(vx,Vy) =MTF(VX,Vy) .Q(V, Vy) +N(Vx, Vy) where capital letters indicate functions in the frequency domain, vx and vy are frequency variables. The function MTF is the modulation transfer function of the imaging system and represents the normalized amplitude of the frequency components of the PSF. This is the representation of the PSF in the spatial frequency domain and is experimentally calculated as the normalized discrete Fourier transform of the PSF (King et al., 1991). In order to simplify the notation, onedimensional continuous variables will be used in the discussion that follows. The conversion to twodimensional discrete indices is straightforward. The main limitation of this model is that no single PSF or MTF(v) can describe the image formation process. The PSF depends on the source detector distance and the depth of the source in the attenuating medium. However, the assumption of shiftinvariance has been experimentally supported by the fact that PSFs, and their corresponding MTFs(v), resulting from either the arithmetic or the geometric mean of conjugate opposing projections are approximately invariant in an uniformly attenuating medium, except near the boundaries of the medium (Larssen, 1980; King et al., 1983; 1987; 1991). King et al. (1991) have shown an experimental procedure for calculating the representative MTF(u) of an imaging system from combined conjugate opposing projections in a uniformly attenuating medium and demonstrated that the function provides a reasonable approximation to a shiftinvariant system (King et al., 1991). Figure 32 depicts the variations of the PSF with the sourcedetector 63 LSF FOR 3mm SOURCE OF "9TC Air Water Arithmetic Mean 36 (A, *A2)/2 ~S /_e, (Wo/,e) \ A 00 1" 28 FWMt DeIh .1 S10 0.07 Sc$m 1 cm SCm*Ieem DISTANCE Figure 32. Line spread functions (LSF) as a function of the distance from the face of the collimator (depth) in air and in a water tank of 20 cm in diameter. Line spread functions at right correspond to the arithmetic means of conjugate opposed projections with the source in water. Notice that the arithmetic mean of conjugate line spread functions is approximately depth independent (from Cho et al., 1993, p. 175). 64 distance, the effect of an attenuating medium and the spatial invariance when the PSF is calculated as the arithmetic mean of conjugate opposing projections. Restoration Filters In a linear shiftinvariant model the restoration problem can formally be expressed in the frequency domain as: given a degraded projection P(v), a linear spatially shiftinvariant degradation operator MTF(v) and the additive noise N(v), calculate the original undegraded projection Q(v). This problem can only be solved if there exists an unique nonsingular inverse operator MTF'(v) such that: MTF(v) .P(v) Q*(v) . In SPECT systems the PSF can be approximated by the linear combination of two 2D Gaussians (Nuyts et al., 1993). Consequently, by the properties of the Fourier transform, the MTF(v) has a 2D Gaussian shape and is equivalent to a lowpass frequency filter (Figure 33). Since the inverse MTF(v) is a singular function (for high frequencies MTF(V) > 0 and MTF1'() > w), there is no solution for the above equation. The physical meaning of this singularity and the absence of an accurate solution is that the power spectrum of the degraded image is restricted by the MTF(v) of the image formation system. The higher frequency components of the object were eliminated during the image formation process and then they can not be recovered by restoration of the degraded image. In addition, because of even small amounts of noise, restoration is an ill conditioned problem. This means that trivial perturbations in P(v) (i.e., 65 0.8 t0.6 0.4 0.2 0 0 10 20 30 40 50 60 Pixel (frequency domain) Figure 33. Plots of modulation transfer function (MTF) at three different distances from the face of the collimator. Distances are 18 cm, 24 cm and 30 cm from the highest to the lowest plot, respectively. 66 noise) can produce nontrivial perturbations in the solution (Andrews and Hunt, 1977, p. 114). Also, the presence of noise does not allow an unique solution. If the noise component is nonzero, there is not an unique association between P(v) and Q(v). Therefore, the object Q(v) can only be approximated by a solution Q'(v) which satisfies some criterion defined a priori. A number of restoration approaches has been developed to obtain practical satisfactory solutions (Andrew and Hunt, 1977, p. 114 117). So far, the simplest and most effective methods in SPECT are those that use restoration filters instead of the inverse MTF (King et al., 1991; Boulfelfel et al., 1992). The filtering is carried out by multiplying the Fourier transform of the degraded image by a filter function and then Fourier inversetransforming the result. Restoration filters are bandpass frequency filters which combines the inverse of the representative MTF(V) for restoring low frequencies, with a lowpass frequency filter for removing noisy high frequency components. They avoid the singularity of the inverse MTF(Y) at high frequencies (regularization). In other words, the inverse MTF(v) is regularized by a lowpass frequency filter. Resolution Recovery Filters In general, any bandpass frequency filter applied to the proper frequency range, produces the amplification of some frequency components which were depressed by the MTF(v). Therefore, resolution is partially recovered and the quality of the image improved. Some bandpass frequency filters have been defined by empirical functions such as Gaussians (Madsen and Park, 1985) or their combinations instead of using the inverse MTF(v) regularized by lowpass frequency filters. The parameters of Gaussian 67 functions have to be optimized for each image by means of interactive trials. Because these filters do not require knowledge of the system MTF(V), they are considered as resolution recovery rather than restoration filters. A well known filter of this type is the Canterbury filter (Corefield et al., 1975) defined by C(v) =Ga(v)K.Gm(v) , where G(v) is a Gaussian function and n, m and K filter parameters that must be varied interactively to optimize the resolution recovery. A proper choice of these parameters can improve significantly image resolution. This method was used by Franquiz et al. (1982a; 1982b; 1983) to improve resolution and calculate right ventricular ejection fraction (RVEF) in equilibrium gated radionuclide ventriculography studies (Franquiz et al., 1983). In these studies, the finite detector resolution blurs cardiac structures producing the overlapping of the pulmonary outflow tract with the right ventricle. This is the main difficulty in calculating the RVEF using the equilibrium gated technique (Winzelberg, 1981; Franquiz et al., 1982a). Canterbury filter significantly improved resolution allowing a more accurate calculation of RVEF (Franquiz et al., 1982b). Major disadvantages of this method are the optimization process of filter parameters, that has to be performed by trialanderror, and the increase in the noise level due to the subtraction operation. Wiener Filter Wiener filter is a restoration filter derived from the constrained leastsquares regularization approach to the restoration problem (Gonzalez and Wintz, 1977, p. 197). This is a well known and effective method for 68 obtaining approximate solutions to restoration problems. This method estimates a solution Q'(v) which satisfies a criterion of goodness given by the minimization of a quadratic functional defined by J(v) = IP(v) MTF(v) .Q'(v) lIN(v) 12+A. IC(v) (v) 2 , where III. is the norm, X is a regularization parameter and C(v) is a regularization operator (Galatsanos and Katsaggelos, 1992). By differentiating J(u) with respect to Q'(u) and setting the result equal to zero, the approximate solution can be calculated from: S(v) A MTF(v) .P(v) MarF(v) +.C2(v) The operator C(v) is an increasing function that avoids the singularity of the filter. For high frequencies, X.C2(v) > MTF2(u), the filter goes to zero and the noise is suppressed. For low frequencies, when v > 0, then C(v) > 0 and the filter becomes the inverse MTF(v). The compromise between inverse restoration and noise suppression is determined by the values of C(v) and X. Both values depend on the noise, MTF(v) and the object and degraded images (Galatsanos and Katsaggelos, 1992). The classical Wiener filter (Andrew and Hunt, 1977, p. 151) has been defined by setting X = 1 and C2(v) = N2(v)/Q2(v), where N2(v) and Q2(v) are the power spectrum of the noise and the object, respectively. Then, the Wiener filter W(v) can be expressed as: w(v) = MTF(v) MTF2 (v) +N2 (v )/ (V) 69 The main inconvenience of this filter is that N2(v) and Q2(v) are unknown functions and must be estimated from the degraded image (Penney et al., 1990). In spite of this drawback, the Wiener filter has demonstrated good performance in deblurring and improving image contrast in physical phantoms and clinical studies (King et al., 1984; Yanch et al., 1988; Penney et al., 1990; Boulfelfel et al., 1992). However, the requirement of estimating the object and noise power spectra for each image, has imposed a real limitation in its practical use. Metz Filter The performance of the Metz filter M(v) is equivalent to that of Wiener (King et al., 1984; Penney et al., 1990), but with the advantage that is simpler and easier to implement. It is given by the expression M(v) (1MTF2(N))x MTF(v) where the exponent X is a power factor that could depend on the total number of counts in the image (King et al., 1983; 1984) and determines the frequency at which the filter goes to zero amplitude (Figure 34). Experimental studies (King et al., 1983; 1987; 1991) have shown that a more convenient filtering, in terms of noise and image quality, has resulted when instead of the representative MTF(v) the following generalized exponential function is used B(v) =exp(i ) where p and S are parameters to be optimized by the minimization of the 70 2.5 2 S1.5 n 1 0.5 0 0 10 20 30 40 50 60 Pixel (frequency domain) Figure 34. Plots of Metz filter for power factors of 10, 4 and 2, from the highest to the lowest curves, respectively. 71 normalized meansquare error (King et al., 1991) between the restored and object image. Another method (King et al., 1987) for determining p and S uses interactive visual feedback to allow the operator to select the optimal function from among a family of possible optimal filters. The Metz filter, either with the MTF(v) or with the exponential function, has been successfully used for compensating scattering and geometric collimator response in liverspleen studies (King et al., 1991), cardiac and brain SPECT (Tsui et al.,1994a) and simulated projections of cardiac studies with high liver activity (Franquiz and Shukla, 1996). Tsui et al. (1994b) have reported improved accuracy and better separation between myocardium and the top of the liver by using Metz filter in simulated and clinical 20'T1 cardiac SPECT. Although the advantages of Metz filter in terms of easy implementation and favorable results, it has seldom been used in clinical studies most likely due to the lack of a consistent methodology for determining the optimal filter parameters. This is especially important because in a restoration operation results are very sensitive to the parameters of the filter function. As with the Wiener filter, the definition of the Metz filter is image dependent and somewhat arbitrary. The main theoretical inconsistency of the Metz filter is that its rolloff, which is given by the exponential factor X (Figure 34), only depends on total image counts (King et al., 1983; 1987; 1991). Since images with very different power spectra can have the same number of total counts, this criterion is not sufficient for determining where the filter amplitude should go to zero. For example, 99mTcMDP bone scan and 9smTcHMPAO brain projections can have the same number of counts but very different behavior in both the spatial and frequency domain. 72 Computer Implementation of the Metz Filter For comparative purposes, conventional restoration of SPECT images was conducted with a Metz filter in this research. It was implemented with the representative MTF(v) of the SPECT system. Since restoration filters are 2D isotropic functions, v represents the radial spatial frequency given by: v = (vx2 + y 2)1"2, where Vx and Vy are the horizontal and vertical spatial frequencies, respectively. All filter algorithms and operations were programmed in MATLAB code (Appendix A). The Filter performance was validated using simulated 2D myocardiumliver projections blurred by convolution with the experimental PSF for s9Tc and 201Tl line sources. Determination of the Modulation Transfer Function The first step was the experimental determination of analytical expressions for the representative MTF(v) of the SPECT system. Modulation transfer functions were determined for the energies of 201T1 and sm"Tc. Also, MTFs(v) were determined with the radioactive source in water and in free air. The determination in water includes the influence of scattering and detector geometric response, while the determination in free air only includes the detector influence. Theory. Modulation transfer functions were calculated using a double Gaussian model for the PSF. Nuyts et al. (1993) have demonstrated that the PSF of a SPECT system can be approximated by: PSF(x,y) =A.exp( x2+y )+B.exp(2( 2 2S2 2522 73 In the above equation A, B, S, and S2 are parameters of the Gaussian functions that fit the PSF counts. However, line sources were employed in experimental measurements because they are easier to handle and prepare than point sources. A line spread function (LSF) expression can be easily obtained from a PSF by integrating in one axis LSF(x) = fPSF(x, y) dy LSF(x) = C1.exp( ) +C.exp( ) , 2S1 2S2 where: C, =v4x.AS1, C2 = i. BS 2 Since PSF and MTF(v) are 2D radially symmetric functions, the Fourier transform of the LSF is equivalent to the cross section of the 2D Fourier transform of the PSF. Then, the MTF(v) at any radial direction can be obtained from the amplitude normalized to zero frequency of the Fourier transform of LSF(x): FT[LSF(x)] = v (CqS1.exp(2 (ivS) 2 + '.exp (2 (v)2)) LSF Measurements. Conjugate opposed SPECT projections (128 x 128 pixels) were acquired with capillary sources of 25 cm long filled with approximately 80 MBq of 99"Tc or 201T1. Sources were positioned at the 74 center and in the long axis of the phantom used in experiments. This was a circular crosssection phantom of 25 cm in diameter. Acquisitions were performed with a threeheaded SPECT system (Trionix Research Laboratories, Inc., Twinsburg, OH) with the high resolution parallel hole and the low energy general purpose collimator for ""Tc and 20T1 sources, respectively. Line source images were acquired with the same energy window used in clinical and experimental studies. Five million counts were acquired in each image. The distance sourcedetector was equal to the radius of rotation (24 cm) used in experimental SPECT acquisitions. LSF numerical data were determined from the arithmetic mean of conjugate projections as the counts in a profile line in the x axis crossing the image of the capillary source. Care was taken to align the source along the y axis in the center of the detector field of view. Line spread functions were determined for each detector. The final LSFs were calculated by adding the contribution of each detector. Curve Fitting. Fitting of LSF parameters to experimental count data was carried out using the nonlinear least squares algorithm of Levenberg Marquardt (Appendix A.1). Line sources in water were fitted to a double Gaussian model, while sources in free air were fitted to a simple Gaussian function. The X2 value between experimental LSF count data and their corresponding fitted functions was calculated as a mesure of the goodness of fit. Results. Table 3.2 presents the PSF and LSF fitting results for 9s"Tc and 20T1 line sources in free air and in water. Figure 35 shows the experimental LSF counts overlayed with the fitted functions. Table 3.3 gives the calculated analytical expressions for the MTFs(v) used in this 75 TABLE 3.2 Line spread function (LSF) and point spread function (PSF) parameters fitted to experimental data from line sources of 9"'Tc and T1. LSF LFTLAI LFTLWA LFTCAI LFTCWA C, 5695 1309 35742 10817 S1 2.316 2.354 1.519 1.522 A 981 222 9387 2835 C2 119 426 S 12.126 8.704 B 3.9 19.5 FWHM 5.60 5.91 3.62 3.71 FWTM 9.63 12.44 6.37 7.27 X2 3.2 3.9 3.9 4.7 LFTLAI: Line source of 201Tl in free air. LFTLWA: Line source of 201Tl in water. LFTCAI: Line source of 99mTc in free air. LFTCWA: Line source of "9mTc in water. A, B and Ci are expressed in relative counts and Si in number of pixels. FWHM: Full width at half maximum in pixels. FWTM: Full width at tenth maximum in pixels. 76 201TI in Water 201TI in Air 1600 01600 I  1200 z 1200 O 0 W 800 W 800 400 400 w wJ S0 0 30 20 10 0 10 20 30 30 20 10 0 10 20 30 POSITION IN PIXELS POSITION IN PIXELS 99mTc in Water 99mTc in Air 1600 1600 z z 400 400 0 0 30 20 10 0 10 20 30 30 20 10 0 10 20 30 POSITION IN PIXELS POSITION IN PIXELS Figure 35. Fitted line spread functions (continuous line) and experimental counts (dots). Data obtained with the line source in water were fitted to a double Gaussian model, while data obtained with the line source in air were fitted to a simple Gaussian function. 77 TABLE 3.3 Analytical expressions of fitted modulation transfer functions. Modulation Transfer Function TLAI exp(0.00646n2) TLWA 0.68 exp(0.00668n2) + 0.32 exp(0.17721n2) TCAI exp(0.00278n2) TCWA 0.82 exp(0.00279n2) + 0.18 exp(0.09127n2) TLAI: "'T1 in free air. TLWA: 201T1 in water. TCAI: ""Tc in free air. TCWA: 9"Tc in water. n is the radial pixel number in the frequency domain. 78 research. Figure 36 depicts MTFs(v) functions calculated from the fitted parameters. Results demonstrated that the experimental data were fit very well by the empirical analytical functions. Broader PSFs were observed for 20'T1 sources (Figure 35). This is due to the lower emission energy, that produces a higher scattering contribution for the source in water, and to the collimator which has lower resolution than that used with the 99mTc source. The effect of broader PSFs is to produce MTFs(v) with lower cutoff frequencies (Figure 36), and consequently more blurred and degraded images. Performance of the Netz Filter The performance of the Metz filter was assessed in simulated cardiacliver projections. The optimal power factor of the filter was determined by trialanderror and based on the visual quality of images and minimization of the normalized mean square error (NMSE) between the filtered image and the original nondegraded object image. Simulation of CardiacLiver Projections. Simulated images were equivalent to the left anterior oblique projection (LAO) at 60 degrees. Images were simulated as frames of 128 x 128 pixels in zoom mode (amplification x 2) and pixel size of 0.356 cm. Simulation included four steps. First, it was simulated the myocardium, second the liver, third the blurring due to scattering and detector geometric response, and fourth the addition of statistical noise. The myocardium was simulated by the projection in a plane of two concentric paraboloids. The myocardial wall was described by the volume between the two paraboloids with dimensions of: Myocardial wall : 1.068 cm (3 pixels), 79 1 1 0.8 201TI in Water 0.8 201TI in Air 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 10 20 30 40 50 60 0 10 20 30 40 50 60 PIXEL (frequency domain) PIXEL (frequency domain) 1 1 0.8 99mTc in Water 0.8 99mTc in Air 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 10 20 30 40 50 60 0 10 20 30 40 50 60 PIXEL (frequency domain) PIXEL (frequency domain) Figure 36. Plots of modulation transfer functions calculated from fitted line spread functions. Modulation transfer functions in water include the effects of scattering and detector geometric response, while those calculated in free air only include the effect of the detector geometric response. 80 Outer paraboloid, long axis (L): 8.2 cm (23 pixels), Inner paraboloid, short axis (S): 6.05 cm (17 pixels). The relationship between Y and X dimensions was given by: Y = L 4.L. X/S2. The simulation was carried out with the MATLAB function heart.m (Appendix A.2). The background was 15% of the myocardium activity. The liver projection was simulated by a rectangle of 10.3 cm length (29 pixels) separated from the myocardium apex by 1.068 cm (3 pixels). The livertoheart activity ratio was 1.5. This was simulated using the MATLAB function liver.m (Appendix A.2). The simulated heart and liver images were added to obtain the ideal projection free of scattering, spatial blurring and noise. This was considered as the object image (Figure 37). The blurring due to scattering and the detector geometric response was obtained from the convolution of the object image with the 2D PSF determined from the 99"Tc or 201T1 line source. Point spread functions were calculated from the parameters of Table 3.2 and using MATLAB functions psftc.m and psftl.m (Appendix A.2). Statistical noise was added by the MATLAB function noise.m (Appendix A.2). This function adds the value RANDN(a(x,y)112) to each pixel (x,y), where RANDN is a normally distributed random number and a(x,y) the pixel content before adding noise. Images were scaled to myocardial counts of 100 and 50 per pixel. This represents myocardial statistical errors of 10% and 14%, respectively. These noise levels are similar to those encountered in clinical studies. 81 Figure 37. Simulated cardiacliver projection. Top left: Object image without any degrading effect. Top right: Degraded image by convolution with the 201Tl point spread function and addition of noise. Bottom left: Degraded image by convolution with the 99mTc point spread function and addition of noise. Bottom right: Restored image (99"Tc) after Metz filtering with the optimal power factor. 82 Metz Filtering. Twodimensional Metz filters were implemented by MATLAB functions tcmetz.m and tlmetz.m (Appendix A.3). The filtering process was performed in the Fourier frequency domain. The 2D discrete fast Fourier transform of the degraded image was multiplied by the 2D Metz filter. The restored image was obtained by Fourier inverse transforming the result (Appendix A.3). The power factor was initially selected by an interactive visual feedback procedure in which the criterion was the visual quality of restored images. A second fine adjustment was performed based on the minimization of the NMSE calculated from NMSE [R(x,y) O(x,y)]2/1 O'2(xy) x1 y1 x1 y1 where R(x,y) and O(x,y) are the restored and object images, respectively. The MATLAB function nmse.m (Appendix A.3) was used for calculating NMSE values. Results and Discussion Table 3.4 and Figure 38 show NMSE values as a function of the power factor X for 9mTc and 201T1 projections with different noise levels. Notice that the optimal power factor value depends on the degradation operator (PSF) rather than on the number of counts (Table 3.4 and Figure 38). The best restoration was obtained with X = 4 for images degraded by the s""TcPSF. This value was independent on the noise level (Table 3.4 and Figure 38). Restored images were very blurred when the power factor was below the optimal value, while above this value there was a significant noise increase without further resolution improvement. In contrast, when 83 Table 3.4 Normalized Mean Square Error (NMSE) as a function of the power factor X in image restoration using the Metz filter. X 99mTc 99mTc 201T1 201T1 (r.e. = 10%) (r.e. = 14%) (r.e. = 10%) (r.e. =14%) NR 7.198 10.535 8.762 12.357 1 7.583 9.395 10.637 12.920 2 6.212 7.834 8.266 10.281 3 5.928 7.538 7.539 9.430 4 5.878 7.536 7.235 9.058 5 5.901 7.635 7.087 8.876 6 5.954 7.776 7.008 8.787 7 6.020 7.936 6.963 8.748 8 6.093 8.107 6.938 8.740 9 6.168 8.282 6.925 8.752 10 6.245 8.460 6.921 8.777 11 6.323 8.639 6.923 8.811 12 6.401 8.819 6.929 8.853 13 6.480 8.999 6.939 8.900 14 6.557 9.178 6.952 8.951 15 6.634 9.357 6.966 9.006 16 6.711 9.535 6.983 9.063 17 6.786 9.712 7.000 9.123 18 6.861 9.889 7.019 9.184 19 6.935 10.064 7.039 9.246 20 7.009 10.239 7.059 9.310 NR: No restoration; r.e.: Relative statistical error per pixel. 84 13 13 11.5 11.5 j 10 I 10 z 8.5 z 8.5 5.5 5.5 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 POWER FACTOR (X) POWER FACTOR (X) 13 13 11.5 11.5 U 10 \ 10 z8.5 1 8.5 i..A .. < 7 7 5.5 5.5 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 POWER FACTOR (X) POWER FACTOR (X) Figure 38. Plots of normalized mean square error (NMSE) as a function of Metz power factor X for 9SmTC (top graphs) and 2 T1 (bottom graphs) projections with different noise levels. Left graphs correspond to relative statistical error of 10%, while right graphs to relative statistical error of 14%. The horizontal line represents the NMSE value obtained without Metz filtering. 85 projections were degraded by 20'TlPSF, optimal restoration was obtained for power factors of X = 10 (for 10% statistical error) and X = 8 (for 14% statistical error). In these images, power factors higher than the optimal value did not significantly increase the NMSE (Figure 38). This can be explained by the flat response of the power spectrum at high spatial frequencies (King et al., 1983) and because in these images most significant image features are contained in the lowest frequency region of the spectrum. It is known that in this frequency region the amplitude of the Metz filter is relatively independent on the power factor X (Figure 3 4). This behavior is significantly different from that of 99"Tc restored projections in which NMSE values critically depend on the power factor (Figure 38). The explanation of this different behavior is given by the fact that 99mTc images contain significant features in the intermediate frequency region, where the amplitude of the Metz filter is critically dependent on the power factor (Figure 38). Another important fact is that as the number of count increases more resolution recovery occurs (lower NMSE) (Figure 38). By increasing image counts the object can be restored from the blurred noise at higher spatial frequencies without significant amplification of statistical noise. This means that as the number of counts increases, the cutoff frequency can be moved to higher frequencies. However, this effect is not only dependent on the noise level, but also on the image structure and degradation operator (PSF). Conclusions Results illustrate the well known fact that Metz filter restoration is image dependent and their parameters (MTF and power factor X) establish 86 a critical compromise between deblurring and noise suppression. In addition, optimal restoration using Metz filtering has to be determined by trialanderror or interactive selection of filter parameters by an experienced operator (King et al., 1983; 1987; 1991). These problems are due to the lack of a consistent and general methodology for regularization and denoising in single scale representations where it is difficult to separate noise from image features (Laine et al., 1994). An alternative restoration approach based on multiresolution analysis could overcome these drawbacks. This approach and its application to image restoration in SPECT is presented and discussed in the next Chapter. CHAPTER 4 NULTIRESOLUTION RESTORATION ALGORITHM The basic idea of multiresolution analysis is to decompose an image into a coarse approximation or low resolution image and the image details for successive higher spatial resolutions (Vetterli and Kovacevic, 1995, p. 75; Wornell, 1995, pp. 16 21). Notice that in this context (multiresolution analysis) spatial resolution means the number of samples or pixels which represent the image. Multiresolution decomposition in the frequency domain is equivalent to dividing the frequency spectrum of the image into a lowpass frequency subband image (the coarse approximation) and a set of bandpass subband images. A particular class of multiresolution analysis uses orthonormal wavelet basis or analyzing functions defined in both the space and frequency domain to decompose images into a set of subband images of constant frequency bandwidth in a logarithmic scale (Daubechies, 1982, p. 7 11; Vetterli and Kovacevic, 1995, p. 69 71). Then, spatial resolution varies with frequency allowing higher sampling rates at higher frequencies (Mallat, 1989). The main advantage of this representation is that it provides orthonormal bases whose components have good localization properties in both the spatial and frequency domain. Variations in resolution for each frequency band allow one to "zoom" image details and localize them through a spacefrequency representation (Daubechies, 1992, p. 27). In other words, this representation provides the spatial location 87 88 of the frequency components for each frequency band (Mallat, 1989). This dual representation has been demonstrated to be a powerful methodology for image analysis and processing (Laine et al., 1994; 1995; Healy et al., 1995; Qian et al., 1995; Sahiner and Yagle, 1995). One particular application has been the denoising of images while keeping or amplifying their significant details at high frequencies (Healy et al., 1995; Laine et al., 1995; Sahiner and Yagle, 1995). It is known that the high frequency content of images is a combination of local features, such as edges, and noise. In a conventional singlescale representation, noise is reduced by using a linear, spatially invariant lowpass filter. This operation also eliminates those high frequency components which correspond to significant image features. Because in some situations highfrequency details are important, their elimination by lowpass filtering may be unacceptable. The dual spacefrequency multiresolution representation provides an improved scheme for noise removal by using spatiallyadaptive operations in each frequency band (Laine et al., 1995). These spatially located operators eliminate noise while preserving image details at high frequencies (Sahiner and Yagle, 1995). Requirements for noise removal in SPECT data are much more simple than those in high resolution images, such as mammography, CT or MRI scans (Laine et al., 1995; Sahiner and Yagle, 1995). Nuclear medicine images are low resolution images in which the most significant features are in low and intermediate scale levels, while the noise is located in high frequency channels. This suggests a simple denoising and restoration procedure by simply discarding channels of high frequency and performing the deconvolution restoration in low and intermediate channels where the 
Full Text 
PAGE 1 A MULTIRESOLUTION RESTORATION METHOD FOR CARDIAC SPECT By JUAN M. FRANQUIZ A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1997 PAGE 2 ACKNOWLEDGMENTS The author expresses special thanks to Dr. Mario Ariet for his support and for providing access to the computer facilities at the Division of Computer Sciences of the College of Medicine, University of Florida, for this research work. The author also thanks Mr. Stefan Lupkiewicz for providing technical help in installing and using part of the software employed in this research. The author is grateful to his academic advisor, Dr. Shailendra Shukla, for the support, criticisms, technical advice and assistance offered throughout this research. The author also would like to express his gratitude to the staff of the Department of Nuclear Medicine at the Veterans Administration Medical Center in Gainesville, for supporting the experimental part and providing access to their nuclear instruments and radioactive sources. The author is specially grateful to Dr. Michael Fagein for his comments and Mr. Morris Thompson for preparing the radioactive sources used in this research and his general help in the experimental work. The author thanks Dr. Nils Diaz and Dr. David Hintenlang for their support in the beginning of the graduate course work at the University of Florida. The author also thanks Dr. Andrew Laine of the Department of Computer and Information Science and Engineering for his initial comments and revision of the proposal of this research and for providing part of the software used in this work. Finally, the author wishes to thank PAGE 3 supervisory committee members Dr. Wesley Bolch, Dr. Walter Drane, Dr. Janice Honeymann and Dr. William Properzio for their assistance and the time dedicated to the revision, comments, criticisms and recommendations throughout this research. PAGE 4 TABLE OF CONTENTS Page ACKNOWLEDGMENTS 11 LIST OF TABLES vi LIST OF FIGURES vi i i ABSTRACT xi CHAPTERS 1 INTRODUCTION 1 The Problem 3 Hypothesis and Objectives 8 Content of Chapters 9 2 CARDIAC SPECT 11 Radiopharmaceuticals 11 SPECT Data Acquisition 14 SPECT Reconstruction 19 Myocardial Tomographic Slices 25 Degrading Factors in SPECT 28 Artifacts in the Inferior Myocardial Wall 45 3 THE RESTORATION PROBLEM IN SPECT 48 Degradation Models in SPECT 49 Restoration Filters 64 Computer Implementation of the Metz Filter 72 4 MULTIRESOLUTION RESTORATION ALGORITHM 87 Wavelet Multiresolution Representation 89 Multiresolution Restoration 99 5 RESTORATION IN A REALISTIC CARDIAC CHEST PHANTOM Ill Material and Methods 112 Results and Discussion 119 Conclusions 156 iv PAGE 5 6 CONCLUSIONS AND FUTURE WORK 159 Other Potential Applications 161 Future Developments 163 APPENDIXES A COMPUTER IMPLEMENTATION OF A METZ FILTER 170 Al Calculation of MTF Parameters 170 A2 Simulation of a CardiacLiver Projection 175 A3 Metz Filtering 177 B MULTIRESOLUTION DECOMPOSITION AND RESTORATION ALGORITHM 179 Bl Multiresolution Decomposition 179 B2 Multiresolution Restoration 181 C MULTIRESOLUTION RESTORATION OF SPECT PROJECTIONS 184 CI Restoration of SPECT Projections 184 C2 Analysis Programs 187 REFERENCES 192 BIOGRAPHICAL SKETCH 202 v PAGE 6 LIST OF TABLES Table Page 2.1 Contribution of outofplane scattered photons per voxel as a function of the activity ratio of outofplane to emission voxels (ROE) for different neighboring planes 42 3.1 Restoration Models in SPECT and their approximate solutions 53 3.2 Line spread function (LSF) and point spread function (PSF) parameters fitted to experimental data from line sources 75 3.3 Analytical expressions of fitted modulation transfer functions .. 77 3.4 Normalized mean square error (NMSE) as a function of the power factor X in image restoration using the Metz filter 83 5.1 Quantitative analysis of myocardial bull's eye polar maps 123 5.2 Quantitative analysis of a midventricular shortaxis slice 125 5.3 Normalized chi square of bull's eye polar maps as a function of the attenuation coefficient and livertoheart activity ratio (LHAR) 129 5.4 Uniformity of bull's eye polar maps as a function of the attenuation coefficient and livertoheart activity ratio (LHAR) for the cardiacliver phantom in water 130 5.5 Normalized chisquare calculated with different reconstruction protocols as a function of the livertoheart activity ratio 135 5.6 Uniformity of bull's eye polar maps calculated using different reconstruction protocols 137 5.7 Quantitative analysis of a midventricular shortaxis slice of the cardiac insert in water without liver activity 141 5.8 Quantitative analysis of bull's eye polar maps of the cardiac liver insert in water 142 vi PAGE 7 5.9 Quantitative analysis of the effect of multiresolution restoration on the normalized chi square of bull's eye polar maps for different livertoheart activity ratio (LHAR) and reconstruction protocols 143 5.10 Quantitative analysis of the effect of multiresolution restoration on uniformity of bull's eye polar maps as a function of the livertoheart activity ratio for different reconstruction protocols 146 5.11 Contrast of myocardial defect measured in a shortaxis slice for different reconstruction protocols 148 5.12 Results of geometric calculation of true defect size on the bull's eye polar map 150 5.13 Defect size in percent of the polar map area calculated from the histogram of counts of bull's eye polar maps as a function of the threshold level 150 5.14 Threshold values for calculating defect size in cardiac SPECT 154 6.1 Data of the three patients included in a preliminary evaluation of the multiresolution restoration algorithm in attenuationcorrected cardiac SPECT scans 165 vi i PAGE 8 LIST OF FIGURES Figure Page 21 Density function g(x,y) in the transaxial plane XY and its projection P 9 [r] at an angle 8 20 22 Filtered backprojection technique 21 23 SPECT reconstruction convolution kernels and their frequency response 26 24 Myocardial tomographic slices 27 25 Detection geometry of unscaterred and outofplane scattered photons 38 26 Plot of the relative number of outofplane scattered photons per voxel as a function of scattering angle 41 31 Coordinates for source distribution and detection planeD 57 32 Line spread function (LSF) as a function of the distance from the face of the collimator (depth) in air and in a water tank of 20 cm in diameter 63 33 Plots of modulation transfer function (MTF) at three different distances from the face of the collimator 65 34 Plots of Metz filter for power factors of 10, 4 and 2 70 35 Fitted line spread functions and experimental counts 76 36 Plots of modulation transfer functions calculated from fitted line spread functions 79 37 Simulated cardiacliver projection 81 38 Plots of normalized mean square error (NMSE) as a function of Metz power factor X for projections with different noise levels 84 vi i i PAGE 9 41 The "mexican hat" function and its translations and dilatations (wavelets) 93 42 Frequency response of analyzing functions and MTF 101 43 Multiresolution decomposition of a high resolution MRI brain image into 5 resolution levels 102 44 Multiresolution decomposition of a low clinical brain SPECT projection into 5 resolution levels 103 45 Square 12 norm expressed in percent as a function of the subband image m 107 46 Count profiles through the simulated projection 110 51 Decomposition of a projection into 5 images or resolution levels 115 52 Bull's eye polar maps of the cardiac phantom without defects and LHAR = 0 120 53 Histograms of polar map counts for the cardiac insert without liver activity 122 54 Horizontal profile lines through the center of a midventricul ar shortaxis si ice 126 55 Graphs depicting the effect of variation in the attenuation coefficient value upon the uniformity of reconstructed bul 1 s eye pol ar maps 131 56 Graphs depicting the effect of variation in the attenuation coefficient value upon the normalized chi square of polar maps ... 132 57 Normalized chi square of bull's eye polar maps as a function of the livertoheart activity ratio (LHAR) and cutoff frequency for four different reconstruction protocols 136 58 Bull's eye polar maps of the cardiac liver insert (LHAR = 3.5) in water 138 59 Graphs depicting the effect of restoration on uniformity of bull's eye polar maps as a function of the livertoheart activity ratio (LHAR) 144 510 Effect of restoration on bull's eye polar maps 145 511 Myocardial tomographic slices for the cardiac insert in water with apical and inferoapical defects 149 ix PAGE 10 512 Graphs depicting the calculated defect size as a function of the threshold level 151 61 Bull's eye polar maps of patient number 3 166 x PAGE 11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A MULTIRESOLUTION RESTORATION METHOD FOR CARDIAC SPECT By Juan M Franquiz December, 1997 Chairman: Shailendra Shukla, Ph.D. Major Department: Nuclear and Radiological Engineering Singlephoton emission computed tomography (SPECT) is affected by photon attenuation and image blurring due to Compton scatter and geometric detector response. Attenuation correction is important to increase diagnostic accuracy of cardiac SPECT. However, in attenuationcorrected scans, scattered photons from radioactivity in the liver could produce a spillover of counts into the inferior myocardial wall. In the clinical setting, blurring effects could be compensated by restoration with Wiener and Metz filters. Inconveniences of these procedures are that the Wiener filter depends upon the power spectra of the object image and noise, which are unknown, while Metz parameters have to be optimized by trial and error. This research develops an alternative restoration procedure based on a multiresolution denoising and regularization algorithm. It was hypothesized that this representation leads to a more straightforward and xi PAGE 12 automatic restoration than conventional filters. The main objective of the research was the development and assessment of the multiresolution algorithm for compensating the liver spillover artifact. The multiresolution algorithm decomposes original SPECT projections into a set of subband frequency images. This allows a simple denoising and regularization procedure by discarding high frequency channels and performing inversion only in low and intermediate frequencies. The method was assessed in bull's eye polar maps and shortaxis attenuationcorrected reconstructions of a realistic cardiacchest phantom with a custommade liver insert and different 99m Tc livertoheart activity ratios. Inferior myocardial defects were simulated in some experiments. The cardiac phantom in free air was considered as the gold standard reference. Quantitative analysis was performed by calculating contrast of shortaxis slices and the normalized chi square measure, defect size and mean and standard deviation of polar map counts. The performance of the multiresolution method was also assessed in 201 T1 clinical SPECT studies. Phantom results demonstrated that the multiresolution algorithm compensated the liver spillover artifact, yielded uniform polar maps and improved significantly the accuracy in calculating myocardial defect size. The procedure does not require operator intervention and can be easily implemented in the clinical setting by using Fourier transform techniques. Finally, the extension of the multiresolution method to other SPECT procedures is discussed and recommended. xi i PAGE 13 CHAPTER 1 INTRODUCTION Singlephoton emission computed tomography (SPECT) is the clinical procedure most commonly used in nuclear medicine. The objective of SPECT is the quantitation of the threedimensional (3D) distribution of a radiopharmaceutical in an organ or region of interest. This is accomplished by calculating, from a set of images (projections) about the human body, the radioactivity concentration in twodimensional (2D) transaxial slices. Radiopharmaceuticals are usually labeled with gamma emitters with emissions between 80 keV and 300 keV. The projection images are acquired with one or more rotating scintillation or gamma cameras. Quantitation in SPECT studies includes two broad categories of experimental procedures. The first is absolute quantitation, which means the accurate and precise measurement of certain quantities in absolute units. Examples are measurements of organ volumes and concentration of the radiotracer expressed in units of activity per unit of volume or mass. Absolute quantitations are needed in absorbed dose calculations for target and normal tissues when radionuclides are administered with therapeutic purposes (Zanzonico et al., 1989; Qian and Clarke, 1996). The other category is relative quantitation, in which relative comparisons of the radiotracer uptake in two regions or at different times are needed. This is the type of quantitation used in diagnostic studies. Examples are the assessment of brain SPECT where the radiotracer uptake in regions of the l PAGE 14 2 left and right hemispheres are compared, or myocardial SPECT where the radiotracer uptake in the cardiac muscle with the patient at rest is compared with the uptake during a physical or pharmacological stress (Berman et al., 1991; Leppo, 1996). Absolute quantitation could be an ideal goal in SPECT studies and has demonstrated great utility in physiological research (Maddahi and Czermin, 1996) and dosimetry calculations (Zanzonico et al 1989). However, experience has demonstrated that absolute quantitation in SPECT is not required for clinical studies. The major use of SPECT in nuclear medicine is for diagnostic purposes. Thus, relative quantitation is a much more relevant problem in SPECT than absolute quantitation. In addition, the particular problem studied in this research is related to myocardial SPECT in which absolute measurements are not needed (Maddahi and Czermin, 1996). Therefore, absolute quantitation problems in SPECT will not be treated in this research. In all the material that follows, SPECT quantitation means relative quantitation. The 3D reconstruction problem in SPECT has been approached using iterative algorithms based on parameter estimation techniques (Shepp and Vardi, 1982; Tsui et al 1994a; 1994b) and noniterative Fourier transform based methods (Larssen, 1980). Iterative algorithms have the potential advantage that accurate quantitative corrections for the three major degrading factors in SPECT (photon attenuation, Compton scatter and finite detector geometric resolution) can be included into the reconstruction process (Buvat et al 1994; Tsui et al 1994a). However, these corrections are done at the expense of such large hardware and computational requirements, that iterative algorithms have remained PAGE 15 3 unsuitable for routine clinical use. Recently, some commercial systems (Case et al., 1996; Cullom et al., 1996; McCartney et al 1996) have included iterative algorithms but only for attenuation correction with results which are simpler than that obtained with scatter and geometric detector response compensations. Fourier transform based methods are computationally efficient but they have the limitation that accurate corrections can not be included into the reconstruction process (Larssen, 1980). Any correction for SPECT degrading factors has to be performed before or after reconstruction (Chang, 1978; Gull berg and Budinger, 1981; King et al 1987; Buvat et al 1994). Among these methods, the filtered backprojection algorithm (Larssen, 1980) is the most easily implemented and so far it is the only SPECT reconstruction technique used routinely in clinical studies. Therefore, because this research is focused on the analysis and solution of a particular clinical problem, only reconstruction and correction methods based on the filtered backprojection technique are considered. The Problem SPECT quantitation is significantly affected by photon attenuation and Compton scatter. Both effects can induce inaccuracies in clinical diagnosis, mostly in myocardial SPECT studies (Nuyts et al 1995; O'Connor et al., 1995a; Ficaro et al 1996; Maniawski et al 1996). Because tissue attenuation is not uniform, the attenuation of photons in some regions produces a larger reduction in the number of counts than in other regions (Tsui et al., 1994a; 1994b). It is known that in 99m Tc SPECT studies, the attenuation effect in some regions can reduce the detected PAGE 16 counts to 25% of the unattenuated counts (Tsui et al 1994a). This artifactual reduction of counts can be interpreted as heterogeneous or pathological uptake distributions. A low clinical specificity (high number of false positive studies) could be obtained if apparent uptake defects due to attenuation are attributed to a true reduced radiotracer uptake. Similarly, a reduced clinical sensitivity (low number of true positive studies) could be obtained if true radiotracer uptake defects are attributed to attenuation. Another important degrading effect is Compton scatter (O'Connor et al 1995a; Welch et al 1995) which produces blurred images and degradation of lesion contrast (Buvat et al., 1994). Consequently, radiotracer defect severity and extent can be obscured or can appear within normal limits because of the scattered radiation that reduces clinical sensitivity. An additional and very important blurring effect is due to the limited spatial resolution of the detector (Tsui et al., 1994a). Photon attenuation as well as Compton scattering depends on the photon energy, the body size of the patient and its composition. This is why they are considered patientdependent effects (Tsui et al 1994a). In addition, the physiological distribution of the radiotracer also influences the magnitude of the scatter effect. This is because the scatter contribution depends on the source distribution (Buvat et al., 1994). Finally, Compton scattering is also dependent on the width and offset of the energy window, the energy resolution of the imaging system and the acceptance solid angle of the detector collimator (O'Connor et al., 1995a; 1995b). PAGE 17 5 It is a widely recognized fact that photon attenuation is the major factor that affects the accuracy of quantitations in SPECT (Tsui et al., 1994a; Ficaro et al 1996). This effect has been an important source of false positive studies of coronary artery disease in myocardial SPECT studies (Ficaro et al., 1996). Although myocardial SPECT without attenuation and scatter corrections has demonstrated high sensitivity in detecting coronary heart disease (Leppo, 1996), its specificity in assessing the absence of disease has been reported as suboptimal in some groups of patients (Ficaro et al., 1996). This has been attributed to the photon attenuation in the thorax which produces false uptake defects in patients with normal myocardial perfusion (Tsui et al 1994a; Ficaro et al., 1996; McCartney et al., 1996). As a result of intensive investigative efforts during the last decade, the new commercial SPECT systems includes hardware and software for correcting attenuation in myocardial studies (Case et al., 1996; Cullom et al 1996; Ficaro et al 1996; McCartney et al., 1996). Although these systems are in the process of clinical validation, preliminary data have demonstrated a significant improvement in specificity when myocardial SPECT is attenuationcorrected (Ficaro et al., 1996; McCartney et al 1996). However, the spillover of activity from nearby organs into the inferior myocardial wall resulting from photon scatter has become a significant problem in attenuationcorrected cardiac SPECT (Ljunberg et al 1994; King et al 1995; Nuyts et al 1995; Cullom et al., 1996; Ficaro et al 1996; McCartney et al 1996). This effect could obscure myocardial defects in the inferior wall yielding false negative results (lower sensitivity). PAGE 18 6 The main source of scattered photons in myocardial SPECT is the radiotracer uptake in liver and bowel (Case et al 1996; King et al., 1995; Nuyts et al 1995). The effect of photon scatter from the liver has been a recognized problem in the assessment of inferior myocardial wall perfusion in some patients, even without attenuation correction (Beller and Watson, 1991; Heo et al 1992; Johnson, 1994). However, this undesirable effect is more severe in attenuationcorrected studies (Ljunberg et al 1994; King et al 1995; Nuyts et al 1995). The reason is that attenuation correction methods can not distinguish between unscattered and scattered photons. Therefore, scattered photons are artificially amplified with attenuation correction. Although the spillover effect has been reported in 201 T1 myocardial SPECT (Case et al 1996; McCartney et al 1996), the effect is particularly important in studies with 99m Tclabeled compounds which have a higher liver and gall bladder uptake (Beller and Watson, 1991; Heo et al 1992; Johnson et al 1994). The general consensus is that some scatter compensation has to be used in conjunction with attenuation correction in order to solve the liver spillover artifact. Different methods have been proposed to compensate for scatter in SPECT (Buvat et al 1994), but only the simplest approaches can be used in the clinical setting because of constraints in computer capacity and processing time. Although there is not a completely satisfactory solution to the scatter problem, even the simplest methods could provide adequate compensation, at least for the accuracy required in most clinical studies (King et al 1987; Frey and Tsui, 1995). PAGE 19 7 One common approach for compensating blurring due to both scatter and detector geometric response is restoration by deconvolution of SPECT data with a representative point spread function of the imaging system (King et al 1983; 1987; Frey and Tsui, 1995). Assuming a linear and shiftinvariant imaging system, the restoration operation is usually performed in the spatial frequency domain, where the deconvolution operator is the inverse modulation transfer function of the imaging system (King et al 1987). However, due to the singularity of the inverse modulation transfer function and the illconditioned nature of the restoration problem, the presence of even small amounts of noise can lead to unacceptable restorations (Andrew and Hunt, 1977, p. 114). Regularization approaches which overcome this drawback have led to the Wiener and Metz restoration filters in nuclear medicine (King et al., 1983; 1987; 1991a). The Wiener filter depends on the power spectra of the object image and noise, which are unknown (King et al 1987; Penney et al., 1990). The Metz filter is easier to implement but there is not a general and consistent methodology for determining its optimal parameters (King et al 1987; 1991). In addition, both types of filters are image dependent. Therefore, the problem of SPECT restoration can be summarized as the optimal design of a lowpass frequency filter for denoising and regularization operations. Optimal filter design is important because the tradeoff between blurring and noise suppression depends critically on filter parameters (King et al., 1987; 1991). This is a common problem in single scale representations where it is difficult to completely separate noise from significant image features (Laine et al 1994). PAGE 20 8 An alternative and relatively new restoration procedure based on mul ti resolution analysis by wavelet expansion (Daubechies, 1992, p. 7 11) is proposed in this research. The advantage of multiresolution representation is that it provides a general spacefrequency framework in which significant image features are more easily identified and segmented from noise (Laine et al 1994). Projections can be decomposed into a set of subband images corresponding to frequency channels of different bandwidth which allow a combined spacefrequency representation of images. In the case of SPECT projections, high frequency or noisy subbands can be eliminated while restoration is performed only in those subband images where significant image features are present. This approach has been previously proposed by Wang et al (1995) for general image restoration and Qian and Clarke (1996) for restoration of bremsstrahlung images of beta emitters. In this research, a multiresolution restoration algorithm is developed for deblurring cardiac SPECT images. The practical contribution of this new approach is that denoising and regularization operations can be performed in a more straightforward and automatic way, without trial anderror optimization (Metz restoration) or power spectra estimation (Wiener filtering). An automatic and general restoration methodology, which does not require operator intervention, would greatly help its practical implementation and application in the clinical setting. Hypothesis and Objectives The ultimate goal of this research is the development and critical assessment of a multiresolution regularization and denoising algorithm for compensation of the liver spillover artifact in attenuationcorrected PAGE 21 cardiac SPECT. The investigation is based on the hypothesis that a multi resolution restoration algorithm can compensate scatter and geometric detector response in noniterative cardiac SPECT reconstructions and with sufficient accuracy to be used in clinical studies. The demonstration of this hypothesis and the assessment of the multiresolution restoration algorithm were performed using a realistic cardiacchest phantom in which a custommade liver insert was included. The methods, algorithms and experimental results are presented in the chapters that follows. All the source code of the software written specifically for this research is included in appendixes. Content of Chapters Chapter 2 describes in detail the methodology of cardiac SPECT reconstructions, correction methods usually employed in the clinical setting and the physical degrading effects that could produce clinical artifacts. Special emphasis is dedicated to those artifacts that are due to the scattering of photons emitted within the liver in attenuationcorrected SPECT scans. Chapter 3 discusses the restoration problem in SPECT and the mathematical models in which restoration is based. The discussion starts with the most realistic and complete approach, which is the nonlinear shiftvariant model. Then, it is demonstrated how practical limitations in the clinical implementation of complex models have led to approximate solutions in which the imaging system is represented by a linear and shiftinvariant model. This chapter also includes a discussion on restoration filters, the implementation and validation of the Metz filter, and the methodology and results of the experimental determination of the modulation transfer functions used in this research. Chapter 4 PAGE 22 10 briefly presents the theory behind multi resolution wavelet representations and discusses its application in restoration of SPECT projections. This Chapter also describes the algorithms used in this research for multiresolution decomposition and restoration, the criteria for denoising and restoration and their validation in a computer simulated cardiac projection. Chapter 5 presents the experimental results obtained with the application of multiresolution restoration to attenuationcorrected SPECT reconstructions of a realistic cardiacchest phantom with different 99m Tc livertoheart activity ratios. Experiments simulated the liver spillover artifact with myocardial homogeneous distribution of the radiotracer (normal subject) and with the inclusion of inferior myocardial defects in the cardiac phantom (pathologic lesions). The performance of the multiresolution algorithm was compared with that of the conventional Metz restoration. Finally, Chapter 6 summarizes the conclusions of this research, discusses a preliminary application of the multiresolution restoration method to clinical studies and makes recommendations for future clinical and physical work. Special emphasis is dedicated to the extension of the multiresolution restoration method to PET and SPECT using a linear shiftvariant restoration model. PAGE 23 CHAPTER 2 CARDIAC SPECT The objective of cardiac SPECT is the relative quantitation and 3D display of the distribution of a radiopharmaceutical in the cardiac muscle. Myocardial radiopharmaceuticals reach the viable myocardium by a blood flow dependent process and tomographic slices provide an accurate indication of the amount of tissue perfusion from the coronary arteries. This information plays an important role in the diagnosis, followup and risk stratification of patients with suspected or known coronary artery disease (CAD), which constitutes the single most common cause of death in the United States (Cullom, 1995). In the last decade, major improvements have been made in radiopharmaceuticals, instrumentation and processing methods in order to improve accuracy in cardiac SPECT (Gait, 1994). Nevertheless, some inaccuracies persist as a consequence of the photon absorption and scatter in tissue. This chapter reviews the basic principles of cardiac SPECT, the newest systems, the physical degrading factors and the correction methods mostly used in routine clinical settings. Specific artifacts encountered in the inferior myocardial wall are discussed in more detail. Rad i opharmaceut i cal s Cardiac SPECT has been performed using 201 T1 for almost 15 years. This is a monovalent cation which is taken up by the viable myocardium by 11 PAGE 24 12 a blood flow and passive transport process (Beller, 1994). The myocardial uptake of 201 T1 starts immediately after its intravenous administration and also begins to redistribute into other tissues and the myocardium itself. Significant change in distribution is observed after more than four hours (Garcia et al 1985). An initial cardiac SPECT depicts the myocardial perfusion at the time of radiotracer injection. A later cardiac SPECT reflects the initial perfusion and the rate of regional radiotracer loss. This biokinetics behavior has been used in the study of transient changes in myocardial perfusion. Clinical protocols usually include two cardiac SPECT scans for each patient (Mahmarian et al 1993; Matsunari et al., 1996a). In the first study, the radiotracer is administered at the peak of a pharmacologic or exercise stress. The second one is performed at rest and after redistribution has occurred. If both stress and rest cardiac SPECT scans show homogeneous distribution of the radiotracer, they are interpreted as normal studies. Regions with significant decreased activity (regional defects) are considered abnormal. If the regional defect is present at rest, this is interpreted as suggestive of myocardial infarction. If the regional defect is present after stress but is either not present or less apparent at rest, then it is considered as indicative of regional ischemia. Unfortunately, the physical properties of 201 T1 are less than ideal for cardiac SPECT imaging. Thallium201 decays by electron capture to 201 Hg emitting characteristic xrays in the 60 keV to 83 keV range. This low energy produces low quality images. Therefore, accurate quantitations are very difficult to make because of photon absorption and scatter in tissue. In addition, the relatively long halflife of 201 T1 (73 hours) does not PAGE 25 13 allow the administration to the patient of high activities, resulting in noisy SPECT images. On the contrary, 99m Tc, the most widely used radionuclide in nuclear medicine, has a monoenergetic gamma emission of 140 keV which is ideal for gamma camera imaging. The shorter halflife (6 hours) of 99m Tc allows the administration of activities that are up to 10 times higher than those of 201 Tl This yields higher quality images in a shorter time period. Another advantage is that 99m Tc is available from a 99 Mo99m Tc generator 24 hours a day, while 201 T1 is cyclotrongenerated and requires offsite delivery. These superior properties of 99m Tc motivated in past decades the search for different compounds to be used as 99m Tclabeled cardiac perfusion agents. In December 1990, the Food and Drug Administration approved the clinical use of two of those compounds: 99m Tcsestamibi (DuPont Nemours, Inc.) and 99m Tcteboroxime (Squibb Diagnostic, now Bracco Diagnostics). Both have shown successful results in clinical trials (Beller and Watson, 1991; Berman et al 1991; Johnson, 1991; Van Train et al 1994). More recently, a new compound, 99m Tctetrofosmin (Amersham International, pic), has also been approved and is used commercially (Braat et al., 1994; Matsunari et al 1996a). The compound 99m Tcsestamibi is a monovalent cation that is taken by the viable myocardium by a passive transport process, similar to the uptake of 201 T1 (Beller and Watson, 1991). This compound has a long myocardial residence time (approximately 5 hours) and little redistribution (Berman et al 1991). In contrast, 99m Tcteboroxime is a neutral lipophilic compound that diffuses rapidly across phospholipid membranes in a manner similar to freely diffusible radiotracers such as PAGE 26 14 133 Xe (Johnson, 1991). Its myocardial washout is flow dependent and shows a rapid bi exponential clearance with halftimes of 3 to 6 minutes and 60 minutes for the early and late components, respectively (Johnson, 1994). The myocardial agent 99m Tctetrofosmin is a lipophilic cationic compound with a rapid heart uptake and relatively slow clearance (Braat et al., 1994). These different redistribution properties have required different clinical protocols from those used for 201 T1 SPECT (Berman et al 1991; Braat et al 1994; Johnson, 1994; Matsunari et al 1996a). One adverse aspect of 99m Tclabeled myocardial agents is the hepatobiliary excretion of these compounds. Photon scatter from hepatic activity up into the heart can interfere with interpretation and quantitation of inferoapical segments, specially in obese patients or in individuals with high diaphragms (Johnson, 1991). Artifacts in the inferior myocardial wall due to liver and abdominal background have been reported in clinical studies with 99m Tcsestamibi (Middleton and Williams, 1993; DePuey, 1994), 99m Tcteboroxime (Johnson, 1991; 1994; Chua et al 1993) and 99m Tctetrofosmin (Braat et al 1994; Matsunari et al 1996b). In spite of their physical advantages, 99m Tc compounds have not replaced 201 T1 in clinical studies. Both radionuclides are currently used in cardiac SPECT. SPECT Data Acquisition The first step in SPECT studies is the acquisition of a set of planar images or projections about the human body using one or more rotating gamma cameras. Typical SPECT studies acquire 60 projections at intervals of 6 degrees. The gamma camera orbits the patient in a stepand PAGE 27 15 shoot motion mode. This means that there is no data acquisition during the time in which the detector is moving between two positions. Projections are converted and stored as computer files with formats of 64 x 64, 128 x 128 or 256 x 256 pixels. The most common format in cardiac SPECT is 128 x 128 pixels. Gamma Camera The gamma camera is an image formation device that uses a scintillator detector with a large Nal(Tl) crystal coupled to an array of photomultipl ier tubes (Cho et al 1993, p. 165). The Nal(Tl) crystal can be either circular, with diameter of 40 cm, or rectangular, with dimensions of 40 cm x 20 cm. The thickness of the crystal is approximately 1.25 cm. The number of photomultipl ier tubes depends on the manufacturer, but presentday gamma cameras use from 37 to 101 tubes arranged in an hexagonal pattern. A multihole collimator, made of a high atomic number substance, is attached to the external surface of the crystal and confines the direction of the incident photons to an extremely small acceptance solid angle. Thus, the crystal receives gammarays only from sources directly in front of it, making a projection from the 3D source distribution onto the 2D face of the scintillation crystal. When a gamma photon is absorbed within the crystal, the photomultipl ier tubes closest to the interaction position receive an amount of light proportional to their distance to the event. The output signals of the photomultipl ier tubes are processed by an analog or digital circuitry to calculate a pair of coordinates (x,y) representing the 2D position of the interaction on the crystal. Finally, the output from each tube is summed to obtain a signal proportional to the total energy absorbed within the crystal. A PAGE 28 16 pulseheight analyzer retains only those events whose signal amplitude lies within a selected small window (+ 10%) around the amplitude of the photopeak. This energy discrimination excludes a high percent of Compton scattering events that occur in the patient, the collimator or in the crystal itself. Spatial Resolution and Sensitivity Two important variables that characterize the performance of a gamma camera are spatial resolution and sensitivity. Spatial resolution indicates the ability of the system to identify the exact location at which a photon has been emitted. This variable is commonly measured in terms of the full width at halfmaximum (FWHM) of the line spread function (LSF) obtained from the count profile recorded across the image of a line source. The FWHM value critically depends on the energy resolution of the detector, the energy window, the sourcetodetector distance, the geometric response of the collimator and the presence of a scattering medium. The best spatial resolution is nearest the collimator. By increasing the distance from the detector to the source, the spatial resolution is significantly degraded. The FWHM in air at most typical organ depths (~ 10 cm) for lowenergy high resolution (LEHR) and lowenergy general purpose (LEGP) collimators are approximately 7 mm and 9 mm, respectively (Tsui et al 1994a). In SPECT studies, at the typical distance from the collimator to the center of rotation (~ 18 cm), the FWHM is approximately 14.5 cm. The presence of a scattering medium (e.g., the patient) should increase the FWHM up to 16.5 mm. Sensitivity is defined for a point source and expressed in terms of counts per minute per unit of activity per meter from the center of the PAGE 29 17 camera detector (Nichols and Gait, 1995). Unlike spatial resolution, sensitivity decreases only slightly with distance. When distance increases, the number of photons per unit of area decreases by the inverse square of the distance, but then the area of the crystal exposed to the source increases by the square of distance. This is valid only in the absence of an attenuating medium. Conventional parallel hole collimators used in SPECT permit only about 0.015% of the emitted photons to interact with the scintillation crystal. For most SPECT studies, the average count rate is on the order of 10 4 counts per second. Because of the limited number of counts per pixel, projections can be severely corrupted by random noise. Sensitivity could be improved by increasing the size of collimator holes or by using more holes of smaller diameter but thinner septal thickness. Unfortunately, there is an inverse relationship between spatial resolution and sensitivity. Both solutions to increase sensitivity allow undesirable photons to penetrate the collimator with the subsequent degradation in spatial resolution. The possibility of increasing the number of counts by increasing the acquisition time is not technically viable. Long acquisition times increase the probability of patient involuntary motion which can give rise to image artifacts. On the other hand, the biokinetics of some radiotracers commonly used in cardiac SPECT do not allow long acquisition times. Increasing the amount of radiotracer also is limited by dosimetry considerations and the dead time of SPECT systems (< 1 /js). Practical solutions to increase the number of counts in cardiac SPECT while preserving spatial resolution have been proposed by modifying PAGE 30 18 acquisition modes (Gait, 1994; Gait and Germano, 1995). The most common modalities have been the following: a) The introduction of elliptical orbits to minimize the distance from the camera to the body throughout all the acquisition is one modality, b) Acquisition of projections in an orbit of 180 degrees instead of the full 360 degrees rotation is another (this method excludes those projections which are severely distorted by photon attenuation and scatter. The acquisition time is then confined to count acquisitions for the most reliable projections), c) The third modality is rotation of the gamma camera and data acquisition in continuous mode rather than the conventional stepandshoot. In this mode, counts are not wasted during the time the gamma camera is moving. In the stepandshoot mode the waste of counts can be significant because detectors take 2 to 4 seconds to move between two positions (Gait and Germano, 1995). Multi Detector SPECT Systems The most important advance for improving sensitivity has been the development of multi detector SPECT systems: the higher the number of detectors, the larger the number of counts that are acquired in the same period of time. The higher sensitivity of these systems has also improved spatial resolution by allowing the use of high resolution collimators (Gait and Germano, 1995). Two opposed detector systems have been used for many years. They were developed for acquiring simultaneous anterior and posterior projections in body skeleton imaging. Threedetector rotating SPECT systems were introduced by Trionix Research Laboratory (Twinsburg, OH) and Picker International (Cleveland, OH) in the late 1980s (Gait, 1994). Today, almost every manufacturer produces threeheaded SPECT systems which are mostly used for cardiac and brain studies. Other systems PAGE 31 19 designed specifically for cardiac SPECT are those with two rectangular gamma cameras mounted at 90 degrees, for example, the Optima gamma camera from General Electric Medical Systems (Milwaukee, WI), the Cardial, from Elscint (Hackensack, NJ), and others (Gait, 1994; Case et al 1996; Cullom et al ., 1996). SPECT Reconstruction The Inverse Problem in SPECT The 3D reconstruction or inverse problem in SPECT can be stated as follows: Given a set of collimated projections {P e [r]} around a transaxial plane t of constant thickness and along an axis r at right angles to the axis of rotation, calculate the density function of the radioactivity concentration g(x,y) in the transaxial plane t that produces {P e [r]> (Figure 21). This problem has been extensively investigated using iterative algorithms and Fourier transform based methods. However, currently only the filtered backprojection (FBP) technique has been routinely used in clinical studies and universally incorporated in commercial SPECT systems. The Filtered Backprojection Technique The FBP is a Fourier transform method based on the application of two mathematical operators to projection data. The first one filters projections using a reconstruction bandpass frequency filter. The second operator, which is the backprojector or backprojection operator, assigns each projection filtered pixel values to all pixels along a perpendicular corresponding line through the reconstructed transaxial image plane. This is done for all the pixels in the projection and over all projection PAGE 32 20 Pe M g (x.y) r = x cose + y sin e Figure 21. Density function g(x,y) in the transaxial plane XY and its projection P[r] at an angle 6. angles. The reconstruction is achieved by the superposition of the filtered projection values for all angles (Figure 22). Briefly, the method can be explained as follows: For a narrow ray, assuming ideal detector response, no photon attenuation and no photon scatter, the relationship between g(x,y) and the set {P e [r]} (Figure 21) is given by the Radon transform (Tretiak and Metz, 1980): where R is the Radon operator. The xy coordinate system is stationary and defines the transaxial image plane. The rs system is rotating and gives g(x,y) 6 (xcos8 +y sin6 r) dxdy PAGE 33 21 0.3 0.2 0.1 0 0.1 0.2 i Original Point Profile Reconstruction (Ramp) Filtered Profile Single Backprojected Profile 10 5 0 5 10 pixel position Superposition of Backprojected Profiles from Multiple Angles (a) (b) Figure 22. Filtered backprojection technique, (a) Convolution kernel or reconstruction filter in the spatial domain; (b) Filtering, backprojection operation and reconstructed transaxial image plane of a central point source. PAGE 34 22 the reference for the projection P[r] (Figure 21). Equations for coordinate transformation between the two systems are given by r = xcos8 + y sin8 s = xsin8 +ycos8. From the Radon transform and equations for coordinate transformation, it can be demonstrated (Rogers and Clinthorne, 1987) that the onedimensional (1D) Fourier transform (FT) of the projection P[r] is equal to the FT of g(x,y) along a radial line in the frequency domain at the same angle as the projection. This is the "Central Slice" theorem, which constitutes the basis of Fourier reconstruction methods. Mathematically, it can be expressed as FT[P 9 [z]] = jfg(x,y) exp [2%j (w x x + w y y) ] dxdy where w x = w.cosfl and w Y = w.sinfl. Then, the above equation can be expressed as FT[P B [r]] = FT[g(x,y) ] e = G(w.B) where this equation indicates that g(x,y) in the frequency domain G(w x ,w Y ) can be synthesized from the radial samples G(w,0), so that: g(x.y) = FT' 1 [G(w x w y ) ] PAGE 35 23 Expressing w x and w Y in polar coordinates (w,0), the above equation becomes g(x,y) = fjG(w, 6) .\w\.exp [2n jw(x. cosd + y. sinQ)] .dwdd where w is the Jacobian of the coordinate transformation. The above integral can be arranged as g(x,y) = fPftlx. cosd +y.sin8] d0 (2.1) where Pftlr] = FT' 1 [FT[P e [r]] \w\] (2.2) Equations (2.1) and (2.2) indicate that g(x,y) is calculated by applying two operators to the projection data. Equation (2.2) corresponds to a filter operator, where P Fe [r] is the filtered projection and w is the ramp filter (Larssen, 1980). Equation (2.1) represents the backprojection operator and indicates the backprojection of P F9 [r] onto the plane xy. Because of the divergent nature of w, the integral (2.2) can not be evaluated. Ramachandrian and Lakshminaryanan (1971) solved this problem by considering that FT[P e [r]] is band limited and then w can be truncated at some value w MAX . Therefore, the ramp filter was modified according to L(w) = \w\ if \w\ < 2nf N L{w) =0 Otherwise, PAGE 36 24 where f N is the Nyquist frequency (Gonzalez and Wintz, 1977, p. 70) defined by and Ax is the onedimensional pixel size. The filtering operation in the frequency domain is equivalent to the convolution of P 9 [r] with the ramp filter in the spatial domain (Gonzalez and Wintz, 1977, p. 61). The convolution kernel L(r) = FT 1 [L(w)] is a symmetrical function with negative sidelobes (Figure 22a) that introduces negative or null values in the filtered projection. These negative values compensate the blurring and star effects produced by the backprojection operation (Figure 22b). Smoothing Window Functions There are two practical limitations in the use of the ramp filter. First, L(w) is a highpass band filter. This is an ideal function to reconstruct noisefree data, but not real SPECT projections which are degraded by statistical noise. Second, the truncated ramp filter introduces rim and overshoot artifacts (Gibbs phenomenon) into the image. This is because the hard cutoff of the ramp filter produces ripples in the filter response to image edges (Algazi et al 1995). Both drawbacks are overcome by multiplying L(w) by a lowpass frequency filter or smoothing window function. Thus, the highpass frequency ramp filter becomes a bandpass frequency filter which gradually rolls off the high frequencies rather than cutting off sharply at f N The most commonly used window functions in clinical SPECT are the Butterworth, Hann, Hamming and Parzen PAGE 37 25 lowpass frequency filters (Larssen, 1980). Figure 23 shows examples of some convolution kernels used in SPECT and their frequency response. The exact form of the convolution kernel depends upon how much lowpass filtering is acceptable in order to reduce the noise in the images while preserving the resolution. The width of the central positive peak determines the amount of smoothing and the negative sidelobes determines the amount that is subtracted from the neighboring backprojected rays. In practical reconstruction algorithms, the window function is a 2D radially symmetric lowpass frequency filter that is applied to projection images before reconstruction. Then, the ramp filter is applied onedimensional ly to the transaxial projection data, row by row. The objective of the first operation is to avoid streaking artifacts in sagittal and coronal slices, which are calculated from transaxial data (Larssen, 1980). These artifacts appear when projection data are smoothed in just one direction. The 2D prereconstruction smoothing avoids streaking artifacts but increases blurring in both transaxial and axial directions. Myocardial Tomographic Slices Coronal and sagittal slices are calculated from transaxial slices for further evaluation of the 3D distribution of the radiotracer (Larssen, 1980). They are displayed in three perpendicular axes parallel to the natural vertical long, horizontal long and short axes of the heart (Figure 24). Axes reorientation is performed by two successive image rotations through angles selected by the operator (Borello et al 1981; Garcia et al 1985). The most important slices from the clinical point of view, are those in the short axis. These slices are normal to the left ventricular long axis and provide a set of crosssectional slices of the PAGE 38 26 Figure 23. SPECT reconstruction convolution kernels (left) and their frequency response (right). From top to bottom: Hamming, Hann and Parzen filters. PAGE 39 27 Figure 24. Myocardial tomographic slices, (a) Standardized display according to the natural axes of the heart (Berman et al 1991); (b) Computergenerated bull's eye polar map display of circumferential maximal count profiles of short axis slices. Apex is derived from both short and vertical long axis slices (Matsunari et al 1996b). PAGE 40 28 cardiac muscle. This reorientation standardizes the cardiac image display and facilitates qualitative and quantitative analysis by exploiting the symmetry of normal left ventricles in the short axis (Berman et al., 1991). Quantitative Analysis Quantitative estimates of the radiotracer distribution are commonly performed using circumferential maximal count profiles of shortaxis slices (Berman et al 1991). An alternative semiquantitative method is the bull'seye plot or display of SPECT data in the form of a polar map (Garcia et al 1985). This procedure converts the circumferential maximal count profiles of short axis slices into polar coordinates profiles. Slices are displayed as a series of concentric circles, with the apex at the center, the inferior wall at the bottom, the anterior wall at the top, the septum at the left and the lateral wall at the right (Figure 24). Pixel counts inside all rings are normalized to the maximum pixel counts and expressed in percent of this value. This representation results in a simple 2D display of the 3D data and allows the comparison of patients's bull'seye with reference normal standards (Ficaro et al 1996; Matsunari et al., 1996b). Degrading Factors in SPECT The FBP method provides accurate 3D reconstruction for idealized projections free of the degrading effects due to the geometric response of the detector, photon attenuation and Compton scatter. When the method is applied to real projections the reconstruction is significantly limited in terms of accurate quantitation, spatial resolution and contrast (Tsui et al., 1994a; 1994b). Therefore, a number of preand postreconstruction PAGE 41 29 semiquantitative techniques have been proposed for compensating distortions when the FBP technique is used in SPECT reconstructions (Chang, 1978; Larsson, 1980; King et al 1991; Buvat et al 1994; Tsui et al., 1994a; Case et al 1996). The discussion that follows, describes the three main degrading factors and the correction methods mostly used in experimental and clinical cardiac SPECT studies. Geometric Response of the Detector The geometric response of the detector is the main cause of blurring and loss of spatial resolution in SPECT studies (Tsui et al., 1994a). This degrading factor critically depends on the geometric response of the collimator and the intrinsic resolution of the scintillator crystal. However, the major contributor to the detector response is the collimator. The number of collimator holes and their acceptance solid angle are the determining factors in the spatial resolution of a SPECT system. The geometric response of the detector is experimentally characterized by the point spread function (PSF), which fully represents the blurring characteristics of the SPECT system. This is a spatially variant function, which rapidly broadens as the distance from the collimator increases. Photon Attenuation Photon attenuation is the most important source of inaccuracy in SPECT quantitations (Tsui et al 1994a; 1994b) and is considered the major cause of falsepositive cardiac SPECT studies (King et al 1996). The attenuation effect can reduce the detected counts to 25 % and 40 % of the unattenuated counts for 99m Tc and 201 T1 respectively (Tsui et al 1994a). Because photon attenuation is a function of the thickness and composition of the absorbing medium, the attenuation effect is not uniform PAGE 42 30 for an imaged organ. This is specially important in cardiac SPECT due to the geometric position of the heart and because surrounding tissues (lung, muscle, adipose tissue, bone) have very different attenuation coefficients. When cardiac SPECT is performed in supine position, there is a significant decrease in myocardial inferior wall counts as a result of attenuation by the left hemidiaphragm (Segall and Davis, 1989; Wallis et al., 1995) and the overlaying right ventricle and right ventricular blood pool (DePuey and Garcia, 1989). However, the problem is not limited to the reduction of counts in the inferior wall. Significant artifacts are induced in the anterior wall due to attenuation from the breast tissue in females (DePuey and Garcia, 1989). Also, in obese patients the accumulation of adipose tissue in the lateral chest wall may result in an apparent perfusion defect in the lateral wall (Wallis et al 1995). Attenuation Correction. The attenuation problem in SPECT can be mathematically described by the inclusion of a negative exponential term in the Radon transform (Tretiak and Metz, 1980). The resultant projection is known as the attenuated Radon transform (Gullberg and Budinger, 1981). So far, an analytical solution to the attenuated Radon transform has not been found. Different attenuation correction numerical algorithms have been proposed, but a complete solution to the problem has not as yet been proposed (Gullberg and Budinger, 1981; Tanaka et al 1984; Tsui et al., 1994a; King et al 1996). However, considerable progress has been achieved in the last decade due to improvements in correction algorithms, hardware of SPECT systems and increased computing power available in current SPECT computers. PAGE 43 31 Depending on their application before, during or after reconstruction, there are three broad approaches for attenuation correction: a) preprocessing correction; b) intrinsic correction; and c) postprocessing correction. Preprocessing methods are based on the assumptions of a uniform concentration of activity and a constant attenuation coefficient in the absorbing medium. One method corrects by hyperbolic functions the arithmetic or geometric mean of conjugate opposing projections (P e [r] and P e+(7 [r]) before filtering and backprojection (Keyes, 1976; Larssen, 1980; Tsui et al 1994a). A practical inconvenience is that the body contour is needed. The simplest and easiest prereconstruction correction makes the reconstruction for 180 with the arithmetic or geometric mean of the conjugate opposing projections which were taken over 360. This correction is employed in commercially available SPECT systems. Intrinsic compensation methods include attenuation correction directly into the reconstruction process. They can be divided into two categories: a) analytical and b) numerical. The first group is based on an approximate solution of the attenuated Radon transform with uniform attenuation coefficient (Gullberg and Budinger, 1981; Tanaka et al 1984). Analytical intrinsic methods use a first correction for the body outline. This is accomplished by multiplying each projection pixel by the inverse of the attenuation along a ray from the projection pixel to the center of rotation. After that, a second correction compensates for the actual attenuation distance. This can be performed in various ways, but the best known method is the Gullberg and Budinger (1981) attenuationweighted backprojection technique. This method uses window functions which depend PAGE 44 32 on a uniform attenuation coefficient (Gullberg and Budinger, 1981). Filtering is performed with a modified ramp filter in which its value is zero in the frequency range below /i/2w (Gullberg and Budinger, 1981). Then, filtered projections are backprojected and multiplied by an attenuation weighting factor which depends on the position of the pixel in the reconstruction matrix. One important limitation is that these algorithms amplify noise and introduce artifacts in lowcount images. Bellini et al (1979) have proposed an analytical solution for the attenuated Radon transform, which relates the nonattenuated and attenuated projection sinograms in the frequency domain. The sinogram is an image which represents the projection data as a function of projection angles. The sinogram is formed by taken the same row of data from all projection images and creating a new image in which the xcoordinate values are those of the projection and the yvalues are determined by the projection angle. The image is called "the sinogram" because a single point should follow a sine curve in this type of representation. Thus, the relationship described by Bellini is given by P'(v,8) =P(v / 8) [(v^J^e+isinhMh/v)] where P' and P represent the attenuationcompensated and original frequency domain sinograms, respectively. The Bellini's method requires the interpolation of the sinogram in the frequency domain to evaluate P{v,6) at the frequency value {v 2 + n 2 ) V2 and at complex angle 6 + i .sinh" 1 (j//jt) After interpolation in the v and 0 directions, the Fourier inversetransform of the sinogram is performed to obtain the attenuationcompensated sinogram. The advantages of this PAGE 45 33 method are its speed and that can be used with the FBP technique. Disadvantages are that corrections are only valid for convex regions of uniform attenuation around the center of rotation and ignores the problem of distancedependent resolution. In addition, the Bellini's method has not been tested extensively as other methods and its robustness in clinical studies is unknown (King et al 1995). Intrinsic numerical methods incorporate a nonuniform attenuation coefficient map into an iterative statistical estimation algorithm such as the maximum likelihoodexpectation maximization (MLEM) or the weighted least squaresconjugate gradient (WLSCG) (Tsui et al 1994b). The MLEM method assumes a Poisson distribution in the projection counts and determines the source distribution that most likely reproduces the projection data (Shepp and Vardi, 1982). Briefly, the MLEM method works as follows. An initial estimate of the transaxial source distribution is performed, e.g., using the FBP technique. Next, each transaxial pixel is updated according to the iteration scheme given by P Trl gHx,y) g H (x,y).B[ l 3 where B is the backprojection operator and the index i enumerates the iteration. The previous transaxial source distribution g'' 1 (x,y) is corrected by the backprojection of the ratio of the actual projection counts P[r] over the projection counts P e '" 1 [r] which were calculated assuming a source distribution given by g ,1 (x,y). In this scheme, the backprojection operator B weights the projection counts by the likelihood that they actually contribute to the source distribution. The 3D PAGE 46 34 implementation of this iterative scheme is straightforward. In theory, the major advantage of this method is that the physical image formation process (attenuation, geometric collimator response and even scatter) can be modeled into the projector and backprojector operators. The inclusion of degrading effects into the reconstruction scheme allows more accurate corrections and source distribution quantitation (Tsui et al 1994a; 1994b). This method have two major disadvantages. The first is the processing time which can be excessive for clinical studies, specially if 3D collimator geometric response corrections are included. The second is the slow and image dependent convergence, in which the noise increases as the number of iterations increases. Other iterative reconstruction algorithms (Tsui et al 1995) have been proposed to speed the convergence and reduces noise, but so far the MLEM algorithm is the most widely used and the standard for accurate quantitation comparations with other reconstruction methods. When only attenuation correction is included in the MLEM reconstruction algorithm, proper compensation may be achieved after a limited number of iterations. Recently, several commercial system (Vantage system from ADAC Laboratories, Optima Lshaped Cardiac Camera from General Electric Medical Systems, and PRISM 3000 SPECT Camera from Picker Inc.) have included the MLEM reconstruction algorithm in conjunction with transaxial attenuation maps for performing accurate nonuniform attenuation corrections. Attenuation coefficient maps have been simultaneously determined using point or line sources in multi detector SPECT systems. At the present time, these systems are in process of clinical validation (Case et al., 1996; Cullom et al., 1996). PAGE 47 35 Chang Attenuation Correction Algorithm Postprocessing techniques are based on the original Chang attenuation correction algorithm (Chang, 1978). This is the best known attenuation correction procedure and was specifically designed to work in conjunction with the FBP method. Most commercial SPECT systems have incorporated this correction technique. The algorithm is implemented by dividing each reconstructed pixel value by the average attenuation factor for that pixel (Chang, 1978). The average attenuation factor is calculated as the average attenuation along all rays from the pixel to the boundary of the attenuating medium. Therefore, the correction factor for attenuation at a point (x,y) is given by cix.y) = £ £ expfu..!*) where M is the total number of projections and l k is the distance between point (x,y) and the boundary point of the medium at projection k. When the source is extensively distributed, the above correction factor overor undercorrects some parts of the image, depending on the source distribution. Next, a second correction is needed. This correction is performed by projecting the corrected data to form a new set of projections. A set of error projections is obtained by subtracting each corrected projection from its corresponding original projection. Error transaxial slices are reconstructed by using the FBP method and corrected for attenuation again. These slices are added to the initially corrected slices to form the final attenuationcorrected image. This correction can PAGE 48 36 be repeated using more than one iterative step. Experimental studies have demonstrated that the Chang algorithm gives the best results after only one or two iterations (Tsui et al., 1989). For nonuniform attenuation, the attenuation coefficient distribution can be calculated from an attenuation map (Tsui et al., 1989; Manglos et al 1987). This method provides a significant improvement in image quality and acceptable quantitation for practical purposes (Manglos et al., 1987; Tsui et al 1989; 1994b). Photon Scattering Compton scattering is the dominant interaction in tissue for the energy range (50 keV to 300 keV) used in nuclear medicine (Evans, 1985, p. 714). In addition, small scattering angles are more probable than large deflections for the emission energies of 201 T1 and 99m Tc. Hence, a relatively high number of forward deflections can pass through the collimator holes. In theory, scattered photons should be separated from those unscattered using the pulse height analyzer of the detector. However, because the finite energy resolution and small sensitivity of SPECT detectors, a certain number of scattered photons has to be counted in order to detect the largest possible number of unscattered events (O'Connor et al 1995a; Welch et al 1995). Photon scatter accounts for approximately 40 % and 60 % of the total counts for 99m Tc and 201 T1 cardiac SPECT studies, respectively (O'Connor et al 1995a; 1995b). This undesirable effect degrades the reconstructed slices by blurring fine details and lowering contrast, and leads to inaccuracies in the quantitative estimate of perfusion defect size (O'Connor et al 1995a; 1995b; Case et al 1996; Cullom et al 1996; McCartney et al 1996). PAGE 49 37 For projections acquired with parallel hole collimators, photons emitted in a given transaxial plane can be detected as emitted in other transaxial plane due to Compton scattering and the finite detector resolution. This effect could be significant in cardiac SPECT where the outofplane scattered liver photons are detected in the inferior myocardial planes (Smith, 1994; King et al 1995). A simple scatter model based on the KleinNishina (Evans, 1982, p. 683) differential probability for a scatter angle and incident energy can give an estimate of the contribution and importance of outofplane scattered photons on a given transaxial plane. It is known (Evans, 1982, p. 675) that the relation between the initial energy E 0 and final energy E f of a photon experiencing Compton interaction is given by B = E £ 1 + a (1 cosq>) where
exp(0 0534.n/sin.9) / 1.593.n 2 ) = p..m v ( ) ( Aft > Aft where, p e = 3.343 x 10 electrons/g is the electron density of water (Attix 1986, p. 531), m v = 1 g/cm 3 x (0.356 cm) 3 is the mass of a voxel,
r' + \i a f$(r') .K(r.r') .d*r' where 0(r) is the particle flux (photons/cm 2 s) at r, (cm 1 ) is the Compton attenuation coefficient (Attix, 1986, p. 132) of the medium and K(r,r') is the kernel of the integral which represents the particle flux at r due to a unit strenght, monoenergetic, steady and isotropic point source at r'. The above integral equation is a homogeneous Fredholm linear equation of the second kind whose solution is given by the Neumann expansion (Zwillinger, 1996, p. 430)
