I M P L E M E N T A T I O N OF A N A N A L Y T I C A L L Y B A S E D S C A T T E R C O R R E C T I O N IN S P E C T R E C O N S T R U C T I O N S By Eric Vandervoort B.Sc, University of Guelph, 2000 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F P H Y S I C S A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A January 2004 Â© Eric Vandervoort, 2004 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Eric Vandervoort 23/02/2004 Name of Author (please print) Date (dd/mm/yyyy) Title of Thesis: Implementation of an Analytically Based Scatter Correction in SPECT Reconstructions Degree: Master of Science Year: 2004 Department of Department of Physics and Astronomy The University of British Columbia Vancouver, BC Canada A B S T R A C T Photon scatter is one of the main effects that degrade the quality and quantitative accu-racy of SPECT images. The objective of the work performed in this thesis was to develop and validate a scatter correction technique that uses an accurate physics-based model of photon scatter as part of an iterative image reconstruction algorithm. This work was performed in two distinct stages. The first stage was to develop a scatter calculation based on a simplified version of an existing technique, the analytic photon distribution (APD) method. The scatter distributions generated using the approximate technique were compared to those obtained using the original APD method. On average the differences between the two methods were small relative to the statistical uncertainty in SPECT projection data (between 6 to 30 times smaller depending on the activity distribution) for a typical perfusion study. The second stage was to implement this scatter correction in an image reconstruction algorithm, and test its performance using computer simulated projection data, experi-mental data obtained from physical phantoms, and patient data. Images corrected for scatter, attenuation and collimator blurring were compared to images corrected only for attenuation and collimator blurring. In the simulation studies results were compared to an ideal situation in which only the primary (unscattered) photon data is reconstructed. In all cases, the scatter corrected images demonstrate substantially improved image con-trast relative to no scatter correction. In simulations the scatter corrected and primary-only images had similar contrast (to within 3%) and noise properties. Images with no scatter correction had contrast values 29% poorer on average than the ideal correction. Currently this scatter correction takes 3 to 4 hours for a typical clinical myocardial study. In the future additional code optimization and some improvements in computer processor speed could allow this scatter correction to be implemented clinically. ii T A B L E OF CONTENTS Abstract ii Table of Contents iii List of Figures vii List of Tables xiv Acknowledgments xvii 1 Introduction 1 1.1 Outline of Thesis 3 1.2 Photon Interactions in Matter 5 1.2.1 Photoelectric Absorption 6 1.2.2 Rayleigh Scatter 6 1.2.3 Compton Scatter 7. 1.3 Photon Attenuation 8 1.4 Nuclear Medicine . . . 10 1.5 SPECT Imaging 10 1.6 Clinical Applications of SPECT 13 1.6.1 Myocardial Imaging 13 1.7 Image Reconstruction Methods 14 1.7.1 Filtered Back Projection Reconstruction 14 1.7.2 Iterative Reconstruction Algorithms 16 iii 2 Literature Review 22 2.1 Holistic Compensation for the Effects of Scatter 23 2.2 Energy-Based Techniques 25 2.2.1 Single Photopeak Methods 25 2.2.2 Dual and Triple Energy Window Methods 26 2.2.3 Multiple Energy Window Methods 28 2.3 Spatial Modeling of the Scatter Distribution 30 2.3.1 Spatially Invariant Scatter Point Spread Functions 31 2.3.2 Slab Derived Scatter Estimation 32 2.4 Theoretically-Based Scatter Corrections 33 2.4.1 Effective Scatter Source Estimation 34 2.5 Monte Carlo-Based Scatter Corrections 40 2.6 Complete Analytical Scatter Corrections 44 2.7 Conclusion 45 3 The Method 47 3.1 APD Method 48 3.1.1 Primary Photon Distribution 48 3.1.2 First Order Compton Scatter 51 3.1.3 Second Order Compton Scatter 54 3.1.4 Higher Order Compton Scatter 56 3.1.5 Rayleigh Scatter 57 3.2 The APD interpolation (APDI) method 58 3.2.1 Defining the Activity Map Subset 60 3.2.2 Estimating the SPSFs using Tri-Linear Interpolation 61 3.2.3 The Interpolated Total Scatter Distribution Function 66 iv 3.2.4 Iterative Calculation of the Scatter Distribution 66 3.3 Implementation of Scatter Corrections in SPECT Reconstructions . . . . 67 3.3.1 The OSEM-BB Method 67 3.3.2 The OSEM-ISR Method 68 3.3.3 The OSEM-APDI Method 68 3.3.4 The OSEM-APDt Method 72 4 Validation of A P D I and O S E M - A P D I Using Computer Simulations 74 4.1 Digital Phantoms and Simulation Acquisition Parameters 76 4.2 Comparison of the full APD method to APDI 81 4.2.1 Results 81 4.2.2 Discussion 86 4.2.3 Conclusions 89 4.3 Optimizing the OSEM-APDI Reconstruction Algorithm 90 4.3.1 Data Analysis 90 4.3.2 Analysis of the OSEM-APDI(0%)P Method 92 4.3.3 Analysis of the OSEM-APDI(15%) Method 93 4.3.4 Analysis of the OSEM-ISR, OSEM-APDt, and OSEM-BB recon-structions 94 4.3.5 Results 94 4.3.6 Discussion 100 4.3.7 Conclusions 108 4.4 Coronary Perfusion Simulations 109 4.4.1 Results 110 4.4.2 Discussion I l l 4.4.3 Conclusions 112 v 5 Experiments 114 5.1 Generation of the Experimental APD Look-up Tables 115 5.1.1 Results 119 5.1.2 Discussion 120 5.1.3 Conclusions 124 5.2 Phantom Experiments 124 5.2.1 The Experiment 126 5.2.2 OSEM-APDI reconstruction specifications 127 5.2.3 Results 129 5.2.4 Discussion 133 5.2.5 Conclusions 141 5.3 Patient Experiments 141 5.3.1 Results 142 5.3.2 Discussion 145 5.3.3 Conclusions 146 6 Conclusions and Future Work 147 6.1 Overview 147 6.2 Current Applicability 149 6.3 Future Work 151 A Comparison of APD with SimSET 155 A. l Results 156 A.2 Discussion 160 A.3 Conclusions 162 Bibliography 163 vi LIST OF FIGURES 1.1 Schematic representation of a typical SPECT imaging system . . . . . . . 11 1.2 A schematic diagram of the line-integral projection data acquisition model is shown in figure (a). Figure (b) shows the back projection of this data. Intersection of back-projected lines through the object correspond to pos-sible source locations 16 1.3 Two dimensional representation of a system matrix that models collimator blurring and photon attenuation. The system matrix element Cy repre-sents the probability of detecting a photon emitted from voxel j in the detector bin i. This probability is the product of the collimator response function Pc(x,y) and the attenuation probability Pa(l) 18 2.1 Simulated Energy Spectra for a 9 9 m T c line source in a water cylinder, based on data from [1]. The solid line is for all photons, both scattered and unscattered. The dotted line is for only the scattered photons 26 2.2 Simulated spatial projection data for a line source in a water cylinder, based on data from [1]. The solid line is for all photons, both scattered and unscattered. The dotted line is for only the scattered photons 32 3.3 Schematic diagram of three different photon paths. The primary or unscattered photon is shown as the dash-dotted line. The primary photon travels directly from the source position: s, to the detector at n. The first order scattered photon, shown as the solid line, originates at source position s, and scatters once at t, and is then detected at n. The second order scattered photon (dashed line) scatters twice, once at position u and again at t, and is then detected at n. 49 vii 3.4 Diagram of the parameterization of the APD look-up tables. The parameters are denned with reference to a plane parallel to the detector surface, shown as the x-y plane in the diagram. The vectors s, t, and n are the positions of the source voxel, the last scattering voxel and the detector pixel, respectively. The detector normal is parallel to the z-axis. The APD parameters are four lengths (a-d) and one angle (e). They are shown as dashed lines in the diagram. The parameters (a) and (b) are the distances from the center of the detector pixel to the center of the final scattering voxel measured perpendicular (a) and parallel (b) to the detector normal. The parameters (c) and (d) are the distance from the center of the source voxel to the center of the final scattering voxel, measured parallel (c) and perpendicular (d) to the detector normal. The parameter (e) is the angle between the vector from s to t and the vector from t to n projected into the x-y plane 53 3.5 Transaxial slice through a digital activity distribution for the heart. Activity intensity values are shown on a linear scale from white to gray. Overlaid upon the activity distribution in black is the interpolation grid with 4x4x4 spacing. 61 3.6 Schematic diagram of to illustrate tri-linear interpolation as it was used in the APDI method to estimate the SPSF, represented here as f(s), for an activity voxel centered at s with spatial coordinates sx,sy, and sz. Points for which the value of the two dimensional SPSF(s},n*j) array has been determined are labeled as open circles. Points for which the function is estimated using linear interpolation are shown as solid circles. The detector position is fixed in the y-z plane at x=0 62 3.7 Schematic diagram of a reconstruction procedure that includes an APDI based scatter correction 69 viii 4.8 Schematic diagram to illustrate the steps involved in the APD (or APDI) pro-jection calculation for the data used in this chapter. The figure shows all the intermediate steps which are required to accurately model the simulated data acquisition 80 4.9 Comparison of the (a) and bias value (b) as a function of projection angle, using the activity distributions for the heart, liver and body 82 4.10 Comparison of full APD and approximate APDI scatter distribution functions for the MCAT phantom (sum of heart, liver and body activity distributions) with the detector at the anterior view (360Â°). Horizontal profiles through the maximum count value are shown using linear and logarithmic scales in (a) and (b) respectively 84 4.11 Same as figure 4.10, except that the projections shown are for detector at the posterior view (180Â°) 85 4.12 Schematic diagram to illustrate the steps involved in OSEM-APDI(0%) P and its analysis 92 4.13 Schematic diagram to illustrate the steps involved in the analysis of the OSEM-APDI(15%) reconstruction technique 93 4.14 Schematic diagram to illustrate the steps involved in the analysis of the OSEM-ISR, OSEM-APDt, and OSEM-BB reconstruction techniques 95 4.15 Comparison of the and bias values as a function of projection angle for each scatter iteration of the OSEM-APDI reconstructions. The APDI data for both the OSEM-APDI(15%) and iterative OSEM-APDI(0%) P reconstruc-tions are shown. Figures (a) and (b) are for the noiseless projection data, while (c) and (d) are for the SimSET projection data 97 ix 4.16 Comparison of the contrast and NSD image quality measurements as a func-tion of OSEM iteration, for the OSEM-ISR, OSEM-APDt (where applicable), OSEM-APDI(15%), OSEM-APDI(0%) P and OSEM-BB reconstructions. The contrast between the healthy myocardium and the defect are shown for the noiseless projection data in figure (a) and for the SimSET projection data in figure (b). Figure (b) and (d) show the liver NSD values for the noiseless and SimSET data, respectively 98 4.17 Comparison of short axis slices for the true MCAT activity distribution with each of the reconstructions considered for the noiseless simulated data: OSEM-ISR, OSEM-APDt, OSEM-APDI(0%) 3, OSEM-APDI(15%) and OSEM-BB. . 101 4.18 Comparison of vertical profiles through inferior wall defect region of the true MCAT activity map and each of the reconstructions considered for the noiseless simulated data 102 4.19 Comparison of short axis slices for the true MCAT activity distribution with each of the reconstructions considered for the SimSET simulated data: OSEM-ISR, OSEM-APDI(0%) 3, OSEM-APDI(15%) and OSEM-BB 103 4.20 Comparison of vertical profiles through inferior wall defect region of the true MCAT activity map and each of the reconstructions considered for the SimSET simulated data 104 4.21 Comparison of contrast between the healthy myocardium and defect VOIs for each of the five coronary defect positions using the (a) the noiseless and (b) the SimSET projection data I l l 5.22 Schematic diagram of a hexagonal collimator hole 117 x 5.23 Comparison of the experimentally acquired projection data (a) and the ex-panded 128x128 APDI photon distribution (b) for the line source experiment with sources located at 0, 10.15 and 20.38 cm from the collimator surface. The experimental projections have been reduced to a 128x128 projection matrix size to allow them to be compared to the APDI projections 121 5.24 Comparison of the experimentally acquired projection data (a) and the ex-panded 128x128 APDI photon distribution (b) for the line source experiment with sources located at 39.22, 49.45 and 59.60 cm from the collimator surface. The experimental projections have been reduced to a 128 x 128 projection matrix size to allow them to be compared to the APDI projections 121 5.25 Schematic diagram to illustrate an error introduced when bi-linear interpolation is used to expand from a 64x64 to a 128x128 matrix size. Figure (a) shows a profile extracted from the original 64x64 matrix with pixel length Ax. This profile represents an impulse function of height h. Figure (b) shows the resulting function when linear interpolation is used to increase the grid sampling by a factor of two (i.e., pixel length = ^p). The resulting triangular function has a F W H M equal to Ax 122 5.26 Schematic diagram to illustrate the steps involved in the APDI projection cal-culation for the data used in this chapter. The figure shows all the intermediate steps which are required to accurately model clinical data acquisition 130 xi 5.27 Plot of the ^j- and bias values as a function of projection angle for the phan-2 torn experiments. The fj values are shown for the high and normal statistics projections in figures (a) and (c), respectively. The bias values are shown for the high and normal statistics projection data in figures (b) and (d), respec-tively. The experimentally measured projections were compared to the sum of the primary and scattered APDI(1%) projection data. 131 5.28 Comparison of reoriented vertical long axis slices through the OSEM-APDI(l%), and OSEM-BB reconstructions which used the high statistics projection data. The two images are displayed using different gray scales normalized to the max-imum count value in each reconstruction. The position of the horizontal profiles shown in the lower figure is shown as the dashed line on the reconstructed slices. 134 5.29 Comparison of reoriented vertical long axis slices through the OSEM-APDI(l%), and OSEM-BB reconstructions which used the normal statistics projection data. The two images are displayed using different gray scales normalized to the max-imum count value in each reconstruction. The position of the horizontal profiles shown in the lower figure is shown as the dashed line on the reconstructed slices. 135 5.30 Comparison of the APDI(1%) estimate of the combined primary and scattered photon distribution with the experimental projections acquired for the high statistics phantom study. The image displayed in the top left is the exper-imental projection and the top right is the APDI(1%) estimate. Horizontal profiles through the heart region are shown using linear and logarithmic scales in (a) and (b), respectively. The detector was positioned at projection angle 5.31 Plot of the (a) and bias values (b) as a function of projection angle for the 50.625Â° 139 patient data, 143 xii 5.32 Comparison of reoriented short axis slices through the OSEM-APDI(0.1%) and OSEM-BB reconstructions for the patient data. The plots in the lower portion of the figure show horizontal and vertical profiles through the slices, indicating regions of the plot corresponding to the left ventricle (L.V.), right ventricle (R.V.) and gall bladder (G.B.) 144 6.1 Schematic diagram of the subset of projection angles for which APD look-up table elements would need to be calculated for two different incident photon directions. The problem here is restricted to 2D, the detector energy response and collimator spatial resolution are assumed to be perfect, and data is collected in a symmetric 20% energy window around 140 keV 153 A . l Comparison of the full APD and SimSET projection data. The projection images are shown for the sum of scattered and unscattered photons, for the MCAT phantom with an inferior wall defect and the detector at the anterior view (360Â°). Horizontal profiles through the maximum count value are shown using linear and logarithmic scales in (a) and (b), respectively. Data in the plots are shown separately for all photons, unscattered photons, first order and second order Compton scattered photons 157 A.2 Same as figure 4.10, except that the projections shown are for detector at the posterior view (180Â°) 158 A.3 Plot of (a) and bias (b) as a function of projection angle to compare the SimSET and APD scattered photon data. The solid and dashed lines are for the APD data using scaling factors determined using the weighted average and total sum of the counts, respectively 160 L I S T OF T A B L E S 4.1 Comparison of mean and the mean bias, averaged over all 64 projec-tions and their standard deviations for the heart, liver and body activity distributions 83 4.2 Timing results for one scatter projection calculation, for the APD and APDI methods, using the heart, liver and body activity distributions. The calculations were performed on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM 86 4.3 Comparison of APD and APDI projections using the SimSET /J, values and p values from the MCAT program 88 4.4 Timing results for each step of the OSEM-APDI(15%) and OSEM-APDI(0%)P reconstructions, for a 64x64x24 reconstruction matrix running on a com-puter with a 1.7 GHz processor and 1024 Mb of RDRAM 99 4.5 Comparison of NSD image quality measures for the liver and body VOIs for each reconstruction technique using the noiseless and SimSET projection data. 112 5.6 Comparison of mean FWHMs fit to profiles through both experimental and APDI projection data. The APDI FWHMs were fit to projection data which was generated using 64x64 APDI projection matrix and expanded to 128x128 using linear interpolation. The experimental FWHMs were fit to the 128x128 projection matrix size, which had been reduced from 1024x1024, to allow for its comparison with the APDI projections 120 5.7 Comparison of mean FWHM fit to profiles through APDI line source pro-jection data for the original 64x64 and the expanded 128x128 matrix size 123 xiv 5.8 Summary of volumes, injected activities, and activity concentrations cor-rected for decay for both high and normal statistics scans for each of the phantoms simulated organs 127 2 5.9 Comparison of mean and the mean bias, averaged over all 64 pro-jections, and their standard deviations for the high and normal statistics phantom experiments. The experimentally acquired projection data was compared to the sum of the primary and scattered APDI(1%) photon distributions 130 5.10 Summary of image quality measures calculated from the known activity con-centrations in the phantom with those measured for the OSEM-APDI(l%) and OSEM-BB reconstructions using both high statistics and normal statistics pro-jection data 133 5.11 Timing results for APDI(1%) and OSEM-APDI(l%) running on a com-puter with a 1.7 GHz processor and 1024 Mb of RDRAM. OSEM recon-structions were performed using 128x128x38 matrix size, while the APDI calculation was performed on a 64x64x19 sized matrix and expanded us-ing linear interpolation. There were negligible differences in the CPU time required for high and normal statistics projection data 136 2 5.12 Comparison of mean and the mean bias averaged over all 64 projections and their standard deviations for the patient data 143 xv 5.13 Timing results for the patient data for each of the OSEM-BB, APDI(0.1%) and OSEM-APDI(0.1%) calculations running on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM. OSEM reconstructions were per-formed using 128x128x38 matrix size, while the APDI(0.1%) calculation was performed on a 64x64x19 sized matrix and expanded using linear interpolation 145 6.1 Timing results (1.7 GHz Processor, 1024 Mb RDRAM) for each step of the scatter correction and reconstruction for the simulation, phantom and patient data sets 151 A . l Comparison of the total sum of the counts for the SimSET and APD projection data. Also included is the contribution from primary photons, first and second order Compton Scattered photons. The APD data using both scaling factors, the weighted average (WA) and the total sum of the counts (TSC), are included 159 A.2 Comparison of the SimSET and APD scattered photon data (with both scaling factors) using the mean *j and the mean bias over all 64 projec-tions and their standard deviations 159 xvi A C K N O W L E D G M E N T S First and foremost, I would like to thank my supervisor, Anna Celler. She has been extremely helpful and enthusiastic throughout my entire thesis. She has done an excellent job preparing me to be an independent researcher and she has truly treated me as a colleague. She knows exactly when to encourage her students, when to motivate them to work harder, and when to give them freedom to pursue their own ideas. I know her research will continue and that recent short-sighted bureaucratic decisions will not discourage her. I would also like to thank Glenn Wells who gave me free reign to use all of the work from his Ph.D. and who answered all my questions thoughtfully and promptly. Without his research this thesis would not have been possible. I also appreciate the help of Javed Iqbal as examining committee member. My thanks also go out to members of the MIRG research group including Kat Dixon, Stephan Blinder and Yanxin Pang; Kat for her advice and for all her assistance in the lab, Stephan for teaching me everything I need to know about programming and computer networks, and Yanxin for greatly improving the speed of my programs and making my method as "user-friendly" as possible. I also appreciate the financial support provided by NSERC and Siemens Medical Solutions USA. I would like to thank John Engdhal and Hans Vija for making my stay in Chicago as comfortable and interesting as possible, and for allowing me to experience what industrial research is really like. My thanks also go to my parents, whose unconditional support has really kept me on track and provided me with a real sense of confidence. Dad, I haven't forgotten our deal: When I become a millionaire physicist - funny I haven't met any of those yet -your retirement is on me. My love and thanks also go out to Ally. She is the most kind and understanding person I know, and no matter how crabby research and thesis writing made me, coming home to her always made me feel happy and loved. xvii C H A P T E R 1 I N T R O D U C T I O N Single photon emission computed tomography (SPECT) allows for visualization of three-dimensional distributions of radio-labeled molecules in the body. These distributions permit physicians to investigate an organ's function and changes in metabolic activity due to disease. SPECT currently provides clinicians with mainly qualitative informa-tion obtained from visually inspecting images. In some situations relative quantitative analysis is being used, however in their present form reconstructed SPECT images are not adequate for accurate absolute quantitative evaluation. Clinically viable quantita-tive analysis of SPECT images will create new opportunities in diagnostic medicine for assessment of in-vivo organ function, physiology, and biochemistry. In SPECT imaging a patient is administered a radio-labeled compound. The com-pound may be injected ingested or inhaled and are then selectively metabolized or other-wise taken up by organs or organ systems in the body. A typical SPECT system consists of one or more scintillation detectors, which rotate around the patient to acquire two di-mensional distributions of detected photons (projections) from a number of angles. The acquired projection images are used to reconstruct the three dimensional distribution of radio-tracer in the patient's body. 1 Chapter 1. Introduction 2 Two important effects that degrade the accuracy and quality of S P E C T images are photon attenuation and scatter. Both effects are a result of photon interactions with matter. Attenuation refers to the fraction of photons that are completely removed or diverted from the path to the detector. Attenuated photons are removed through pho-toelectric absorption or are diverted through scattering in which the photons energy and direction of propagation change. A subset of those scattered photons may be recorded in an incorrect part of the camera, providing false information about the source loca-tion and intensity. Both effects depend strongly on the distribution of the attenuating medium in the body, which can be represented by an attenuation map. These maps are best obtained through a transmission scan. In a transmission scan an external source of radioactivity, called a transmission source, is used to measure the attenuation of the object being imaged. Another factor that can negatively influence the quantitative ac-curacy of S P E C T images is collimator blurring. Since the detectors cannot determine the direction of the incoming photons, they are equipped with collimators. Collimators severely limit the sensitivity of the imaging system and cause distance dependent image blurring. The main objective of the work performed in this thesis was to develop an iterative image reconstruction algorithm for S P E C T that employs exact, physics-based, patient-specific corrections for attenuation, collimator blurring, and photon scatter. The goal was to develop an scatter correction for S P E C T image reconstruction that employs the computed analytic scattered photon distribution (APD) method [2, 3]. The A P D method analytically evaluates the probability that a scattered photon will be detected in a given detector bin in S P E C T projections, using the exact Klein-Nishina formula and an accu-rate model to describe the S P E C T data acquisition process. The A P D method is quantitatively accurate and faithfully reproduces the spatial Chapter 1. Introduction 3 distribution of the unscattered and scattered photons. However, in its original form this method was too slow for practical clinical use. In order to implement a clinically viable S P E C T reconstruction algorithm that accurately models the effects of photon scatter, it was therefore necessary to investigate simplifications to the A P D method and determine an optimal way to include this scatter correction in a reconstruction algorithm. This research focused on cardiac applications. There are a number of advantages to using quantitatively accurate S P E C T images, as opposed to the typically used qualitative images, for the diagnosis of coronary artery disease. In S P E C T images in which scatter effects are not taken into account, there is the potential for misleading or ambiguous reconstructed activity distributions. For example, in images of patients with an elevated liver, it is possible that photons originating from the liver are scattered and may be falsely detected as having come from the heart. The apparent elevated uptake in the heart wall could obscure the presence of a perfusion defect in the heart. This ambiguity could be avoided using quantitatively accurate S P E C T images that have been corrected for scatter. 1.1 OUTLINE OF THESIS This thesis consists of six chapters. The following summary briefly describes the purpose and contents of each chapter. C H A P T E R 1: I N T R O D U C T I O N This chapter begins with a discussion of the photon interactions with matter for energies relevant to S P E C T imaging. Also included is an overview of S P E C T imaging, in which the equipment used and a number of common clinical applications are reviewed. The chapter concludes with an introduction to S P E C T image reconstruction. Chapter 1. Introduction 4 C H A P T E R 2: L I T E R A T U R E R E V I E W Chapter 2 consists of a review of current scatter correction techniques which have been investigated in the clinical and research environment. C H A P T E R 3: T H E M E T H O D This chapter begins wi th a review of the A P D method for calculating the scatter compo-nent in S P E C T projections as it was developed by Wells et al . This chapter then describes the A P D I interpolation ( A P D I ) approximation method. This method is a modification of the full A P D method that substantially reduces its calculation time. Final ly, the O S E M - A P D I method is described. This method incorporates the A P D I scatter estimate into the ordered subset expectation maximizat ion ( O S E M ) [11] iterative reconstruction algorithm. C H A P T E R 4: V A L I D A T I O N O F A P D I A N D O S E M - A P D I U S I N G C O M P U T E R S I M U L A -T I O N S This chapter investigates and validates the O S E M - A P D I scatter correction technique for S P E C T reconstructions. The validation was performed in two steps. The first validation was for the A P D I method and consisted of a comparison between the distr ibution of scat-tered photons in S P E C T projections that were calculated using the original A P D and A P D I methods. The second stage of this validation was for the O S E M - A P D I method and was performed by comparing reconstructed images that used computer simulated projection data. The O S E M - A P D I reconstructed images were compared to other scatter compensation techniques. Several measures of image quality and quantitative accuracy were compared for each reconstruction technique. A l l calculations in this chapter used Chapter 1. Introduction 5 realistic computer simulated phantoms that modeled the activity and attenuation distri-butions for perfusion defect studies. C H A P T E R 5: EXPERIMENTS This chapter describes a number of experiments performed to further validate the OSEM-APDI method. These studies used a physical phantom that modeled the heart with a perfusion defect in a realistic thorax. The OSEM-APDI reconstructed images were compared with respect to image quality and their quantitative accuracy relative to the known activity concentrations that were injected into the phantom. Finally, OSEM-APDI was applied using clinical patient data. In this case, the true distribution of activity was not known, so this reconstruction could only be evaluated qualitatively. C H A P T E R 6 : CONCLUSIONS AND F U T U R E W O R K This thesis will conclude with a summary of the work done thus far and the applicability of the OSEM-APDI method to other nuclear medicine applications. A discussion of future work that could further improve the utility and clinical viability of the OSEM-APDI method is also provided. 1 . 2 PHOTON INTERACTIONS IN MATTER Photon scatter and attenuation are a result of gamma ray interactions with atoms be-tween the point of emission in the patient's body and absorption in the SPECT gamma ray detector. There are three main types of interactions that are important for the range of energies used in nuclear medicine. These are the photoelectric effect, Rayleigh scat-tering, and Compton scatter. For low Z organic materials and photons with energies in the range of interest in nuclear medicine (~ 100-500 keV), Compton scattering is the Chapter 1. Introduction 6 dominant process. 1.2.1 P H O T O E L E C T R I C ABSORPTION The photoelectric effect (PE) occurs when a photon is completely absorbed by an atom, which is ionized in the process. The PE effect can occur if the photon energy is sufficiently high to eject an orbital electron from the atom leaving a vacancy in the atomic shell. A transition then occurs in which an electron with a lower binding energy fills this vacancy. This process is accompanied by the emission of characteristic x-rays or Auger electrons. The energy of the characteristic x-ray corresponds to the difference between the binding energies of the two atomic shells involved in the transition. This energy can also be directly emitted in the form of an electron, called an Auger electron, ejected from another of the atomic shells. The probability r, that photoelectric absorption occurs (or cross-section) per atom is approximately proportional to Z 3 , where Z is atomic number, for higher Z materials and to Z 4 , 8 for low Z materials. The PE cross-section is also approximately proportional to E ~ 3 , where E 7 is the incident photon energy. For biological materials and photon energies used in SPECT, the contribution from PE to the total cross-section (probability for all interactions) is small (1.9% for water at 140 keV). 1.2.2 R A Y L E I G H S C A T T E R In Rayleigh scattering, also referred to as coherent scattering, the oscillating electric field of an incident photon interacts with the electrons surrounding the nucleus. This causes the electrons to momentarily vibrate. These oscillating electrons emit radiation with the same wavelength as the incident photon. This radiation is emitted as a scattered photon with the same energy as the incident photon but moving in a different direction. The Chapter 1. Introduction 7 differential Rayleigh scattering cross-section represents the fraction of photons scattered per electron per unit solid angle, and can be expressed as: d ^ { 6 ) = rj(l + cos>e)T(x,Z), (1.1) where 0 is the Rayleigh scattering angle (measured with respect to the direction of the incident photon), and ro is the classical electron radius. The total probability of Rayleigh interaction (total cross-section) per electron ocoh can be obtained by integrating equa-tion 1.1 over the solid angle. The term J-(x, 2), represents the coherent scattering atomic form factor. This factor depends on 2, the atomic number (or effective atomic number) of the scattering medium and on x, the momentum transfer of the scattered photon. The momentum transfer can be written: sin f , where A is wavelength of the scattered photon. The contribution of Rayleigh scatter to the total cross-section for biological materials and photon energies used in SPECT is also small (0.7% for water at 140 keV). 1 . 2 . 3 C O M P T O N S C A T T E R Compton or Incoherent scattering is the most significant photon interaction for energies of interest in nuclear medicine. For example, for 140 keV photons in water Compton Scatter accounts for 97.4% of the total cross-section. In Compton scattering a photon interacts with a free or weakly bound orbital electron. Part of the photon energy is transferred to the electron and both the energy and direction of the scattered photon changes. The Klein-Nishina differential scattering cross-section for Compton scatter can be written as: Chapter 1. Introduction 8 / a 2 ( l - c o s 0 ) 2 \ â€ž , - x where: â€˘ 0 is the Compton scattering angle, â€˘ a is energy of the incident photon divided by the rest mass of the electron, â€˘ S(x, Z) is the incoherent scattering atomic form factor and â€˘ r0 is the classical electron radius. The original Klein-Nishina derivation of the scattering cross-section was for free elec-trons and did not include the atomic form factor S(x, Z). This correction factor depends on Z, the effective atomic number of the scattering medium, and the momentum transfer of the scattered photon x. Integrating equation 1.3 over the solid angle yields oinc, the total probability of Compton interaction per electron. The energy of the Compton scattered photon can be obtained from the principle of momentum and energy conservation. This energy can be written as: E E' = f l 4) where E' and E are the energy of the scattered and incident photon respectively. 1.3 PHOTON ATTENUATION Attenuation can introduce significant artifacts in reconstructed S P E C T images. Specif-ically the reconstructed intensity of regions lying deeper in the body will be underes-timated if no correction is made for attenuation. Photons from regions deeper in the Chapter 1. Introduction 9 body will be attenuated more strongly than those closer to the surface. Attenuation also causes inconsistencies in the measured projections, that is, the number counts recorded are different when viewed from different angles. These inconsistencies cause geomet-rical distortion of shapes in reconstructed images. For these reasons, it is necessary to accurately model the process of photon attenuation in S P E C T reconstruction algorithms. If a mono-energetic beam of photon flux with intensity $ 0 passes through a thickness dx of attenuating material, the fraction of the flux d$ that is removed from the beam or is attenuated is given by where /i is a constant of proportionality referred to as the linear attenuation coefficient. Since attenuation occurs as a result of photon interactions /J, represents the fraction of photons that interact per unit thickness of attenuating material. Each of the three photon interactions mentioned in the previous section can occur independently. The total / J for energies of interest in SPECT, for a given attenuating material, is the sum of the linear coefficients of each interaction. This coefficient can be written as: where ÂŁ i p . e . , P-coh, and p,inc, are the linear attenuation coefficients for the photo- electric effect, Compton scatter and Rayleigh scatter respectively. These values are related to the total cross sections discussed previously. For example, the linear attenuation coefficient of Compton scattering can be obtained from <7;nc, the total cross-section per electron, as follows: d$ = -fx^0dx (1.5) M f^coh ~t~ l^inc ~t~ A^ p.i (1.6) t^inc â€” Â®inc Peleci (1.7) Chapter 1. Introduction 10 where peiec is the electron density per unit volume of the attenuating medium. If equation 1.5 is integrated over thickness d of an attenuating medium, the following expression is obtained for the fraction of photon flux that is transmitted through the attenuating medium: 1.4 NUCLEAR MEDICINE Nuclear medicine provides information about the function and metabolism of organs by detecting the radiation emitted by radio-labeled compounds that a patient has ingested or inhaled or that have been administered to the patient intravenously. The compounds used in nuclear medicine are selectively metabolized or otherwise taken up by organs or organ systems in the body. A s a result the distribution of radio-labeled compounds in the body provides information about the functioning of the targeted organ or organ systems. A typical S P E C T imaging system consists of one or more scinti l lation detectors mounted on a rotating gantry to allow the detection of photons emitted by a distribution of radioactive compounds inside the patient's body from a number of angular views. A schematic diagram of a typical gamma camera is shown in figure 1.1. The components of a gamma camera are the collimator, detector crystal, photomultiplier tube array, position logic and energy discrimination circuits and a computer for data acquisition and analysis. The purpose of the collimator is to determine the direction of incident photons by (1.8) 1.5 S P E C T IMAGING Chapter 1. Introduction 11 Patient Radioactive Source Path of Photon Collimator Detector Crystal Photomultiplier Tube Array Position and Energy Electronics Computer Figure 1 . 1 : Schematic representation of a typical S P E C T imaging system limiting the angular distribution photons that reach the detector crystal surface. The collimator consists of a sheet of dense, highly-attenuating material (often lead or tung-sten) that contain many small circular or hexagonal holes. These holes only transmit photons that are traveling at small angles relative to a vector perpendicular to the camera surface. The maximum angle for which photon detection is possible that is determined from the collimator hole geometry is often referred to as the acceptance angle. There is a small probability that photons with incident angle greater than the acceptance angle will be detected due a septal penetration. In this process photons pass through the walls between collimator holes (septa). A typical S P E C T collimator has an efficiency on the order of approximately 1 detected photon for every 1 0 4 photons emitted. Collimators with larger holes provide higher photon detection efficiency, but increase the acceptance angle and thereby lower the detector spatial resolution. Collimators containing smaller holes have lower detection efficiency, but improved spatial resolution. Photons that pass through the collimator are detected in a scintillation crystal. The most commonly used crystal in S P E C T is thallium-activated sodium iodide (Nal(Tl)). These crystals are highly efficient at converting photon energy into visible light photons. Chapter 1. Introduction 12 The amount of light produced by a Nal(Tl) crystal is proportional to the energy of incident photon, allowing for discrimination of detected events by energy. However, these crystals have relatively poor energy resolution making it difficult to distinguish scattered from unscattered photons from energy information alone. Semiconductor crystals, such as Ge(Li) or high purity germanium, can provide substantially better energy resolution but are not used commercially because of their high cost, low detection efficiency and complicated maintenance (i.e., liquid nitrogen cooling). Although Nal(Tl) crystals are highly efficient at converting radionuclide emission pho-tons into light photons, the amount of light produced is still quite small. The light signal is amplified and converted to an electrical signal by means of an array of photomultiplier tubes (PMT). These tubes form a tightly packed hexagonal array on the back surface of the crystal. The PMT array also provides the spatial localization of the detected photon on a discretized two dimensional grid of camera bins. The amount of electrical signal pro-duced by a given PMT is proportional to its distance from the photon absorption event in the crystal. Therefore, the spatial location of the photon detection can be determined by combining the electrical signal from each photomultiplier tube in a logic circuit that weights the signals appropriately. A computer is used to specify acquisition parameters such as the number and size of detector bins, the number of projection angles, the temporal duration of each view, and other various features related to the motion of the gamma camera. The computer also digitizes and stores the collected data. It is often also equipped with software for image reconstruction and visualization of the projection or reconstructed data. Chapter 1. Introduction 1 3 1.6 CLINICAL APPLICATIONS OF S P E C T It is possible to investigate the function and physiology of a wide variety of organs using SPECT imaging techniques. Some common SPECT applications include imaging of pul-monary emboli, abnormal bone growth, renal function, cancerous lesions and coronary perfusion defects. Inhaled radio-pharmaceuticals such as 1 3 3 Xe or 1 2 7 Xe gas can be used to detect airway obstructions known as pulmonary emboli. These radio-pharmaceuticals bind to the collagen present in an embolus. The resulting SPECT scan provides infor-mation about the size and location of the embolus. In bone imaging small fractures and fissures in bone tissue can be visualized using phosphate complexes (e.g.; 9 9 m T c -phosphate). Phosphate is more highly metabolized in the formation of new bone tissue. The distribution of these radio-pharmaceuticals allow physicians to detect abnormal bone growth. Renal plasma flow and glomular filtration rates can also be investigated with the use SPECT imaging or dynamic planar nuclear medicine scans. Low renal blood flow and filtration rates can be an early indicator of kidney damage. The most common SPECT application and the focus of this research is coronary perfusion studies. These studies will be described in further detail in the following section. There are numerous other applications of SPECT imaging and a more complete description of these and other applications can be found in the following references [4, 5, 6]. 1 . 6 . 1 M Y O C A R D I A L IMAGING The focus of this research has been to evaluate a scatter correction technique with an emphasis on the improvements in separation of viable from nonviable myocardium tis-sue in reconstructed images that have been accurately corrected for scatter. Radio-pharmaceuticals used in myocardial imaging are taken up by muscular heart tissue in proportion to blood flow. In healthy individuals blood flow increases during physical Chapter 1. Introduction 14 exertion. In patients with coronary artery disease blood perfusion at rest is normal but the ability to provide a larger blood supply during exercise may be compromised. By comparing SPECT images for a patient in states of rest and physical stress, a physi-cian can distinguish between normal, under-perfused and necrotic tissue. For example, necrotic tissue demonstrates reduced perfusion in both stress and rest studies but re-versible ischaemia demonstrates decreased uptake only in the stress scans. One of the most common radio-nuclei used for myocardial perfusion studies is 9 9 m T c . This radionu-clide is a low energy (140 keV) gamma emitter with relatively short physical half-life (6 .022 hours). In all the work done in this thesis, 9 9 m T c was the radionuclide under investigation. 1.7 IMAGE RECONSTRUCTION METHODS In SPECT, image reconstruction refers to the techniques used to obtain a three dimen-sional distribution of radio-tracer in the patient's body from a set of projection data, which represents the two dimensional distribution of recorded photons from various an-gular views around the patient. The accuracy of the resulting reconstruction strongly depends on how well the reconstruction technique models the SPECT imaging process. An accurate SPECT reconstruction must include realistic modeling of the data acquisi-tion process, the statistical nature of radioactive decay and photon interactions. 1 .7 .1 F I L T E R E D B A C K P R O J E C T I O N R E C O N S T R U C T I O N Filtered Back Projection (FBP) [7], is commonly used in the clinical environment and employs a very simple model of the SPECT imaging process. FBP is an analytical tech-nique that assumes that the number photons recorded in a detector bin represent a one dimensional line integral through the activity distribution along a vector perpendicular Chapter 1. Introduction 15 to the detector surface. The projection process is shown schematically in figure 1.2 (a). If the activity object distribution being imaged is \(x,y), then the number of counts Y(x',9) in projection bin position x' at angle 9, can be written as: Y(x'6) = j_\(x,y)dx\ (1.9) where x' = x cos6 + y sin6 and / is a vector perpendicular to the detector surface from the detector pixel position x' through the activity object. Back projection refers to the redistribution of these counts back along a line through the object. The points of intersection for the back-projected lines correspond to potential source locations. This process is shown schematically in figure 1.2 (b). The figure shows that a back projected object contains low density streaks radiating out from the points of intersection. The true activity distribution can be recovered only if projections are acquired for an infinite number of angular views between 9 = 0 and 6 = n, and the projection data is convolved with a ^ filter prior to back projection. Here r refers to the polar coordinate r = \/x2 + y2. The expression for the FBP reconstructed activity can be written concisely in terms of the inverse Fourier transform of the original projection data, Yift(i>,9), as: \(x,y)= / \i/\Yift{ut d) exp (2niu(xcose+ ysin9)) dudjd, (1.10) JO J-oo where v is the polar spatial frequency and the term \v\ is the inverse Fourier transform of the i filter. r The main advantage of FBP is the speed with which a reconstruction can be per-formed. However, there are a number of problems associated with FBP reconstructions. One difficulty is that statistical noise in data acquisition is not accounted for. The high frequency noise in the projection data is enhanced by the ramp filter. This noise can be reduced by modifying the ramp filter to suppress high frequency components in the Chapter 1. Introduction 16 Y(x'e) (a) (b) Figure 1.2: A schematic diagram of the line-integral projection data acquisition model is shown in figure (a). Figure (b) shows the back projection of this data. Intersection of back-projected lines through the object correspond to possible source locations. projections. However, these techniques tend to introduce artifacts in the reconstructed images and inevitably suppress high frequency components that are actually present in the imaged object. The most serious problem in F B P is the inherent lack of quantita-tive accuracy in the reconstructed images due to the algorithm's unrealistic model of the SPECT imaging process. For example, there is no way to accurately incorporate photon attenuation or the distance dependent spatial response of the gamma camera into the line integral model of projection data acquisition. 1.7.2 I T E R A T I V E R E C O N S T R U C T I O N A L G O R I T H M S Iterative reconstruction algorithms provide an alternative to analytically formulated re-construction procedures such as FBP. These methods iteratively search for the best esti-mate of the activity distribution taking into account the measured projection data and a particular system model of the S P E C T data acquisition process. The main advantage Chapter 1. Introduction 17 of iterative methods are that effects such as collimator blurring, photon attenuation and scatter can be included directly in the system model. The SPECT data acquisition process can be thought of as a transformation that maps the three dimensional distribution of activity in the patient's body into a two dimensional distribution of counts in the projection data. We can write this transformation as: Yi = J2Cii*j> (i-ii) where Aj is the amount of activity present in object voxel j, Yi are the number of counts detected in detector bin i and the summation over J is over all voxels in the object containing activity. The term is the system model matrix element which represents the probability that a photon emitted from object voxel j is detected in detector bin i. The system matrix can contain geometrical considerations such as the collimator response and other data acquisition parameters. It can also account for effects such as attenuation and scatter which are patient dependent. An example of a simple two dimensional representation of a system matrix, that includes the effects of photon attenuation and collimator blurring, is shown schematically in figure 1.3. The probability of detecting a photon emitted from voxel j in detector bin i is the product of the collimator response function Pc(x,y) and the probability of photon attenuation Pa(l). We can approximate the response of the collimator as a Gaussian function with a distance dependent full width at half maximum (FWHM) (e.g., [8]). The FWHM depends the variable y, the distance from the source voxel to the detector surface measured along a vector perpendicular to the detector face. The exact form of the FWHM depends on the accuracy of the model employed, camera acquisition parameters and on the geometry of collimator used. The probability of detecting a photon emitted from voxel j also depends on the amount of attenuating material between the source Chapter 1. Introduction 18 voxel and the detector pixel. We can write this probability of attenuation as: PQ(f) =exp (- Jp(l')dl^ , (1.12) where / is a vector from the center of activity voxel j to the center of detector pixel i. The term p,(l') is the attenuation as a function of distance along the vector /. Probability of Pa( /) Attenuation Collimator Response Pc(x,y) Yi Attenuating Medium Figure 1.3: Two dimensional representation of a system matrix that models collimator blurring and photon attenuation. The system matrix element represents the proba-bility of detecting a photon emitted from voxel j in the detector bin i. This probability is the product of the collimator response function P c (x, y) and the attenuation probability Pa(l) One of the difficulties involved in the inclusion of an accurate model of photon scatter in an iterative reconstruction algorithm is that with scatter the system matrix becomes non-sparse (i.e., a large number system matrix elements are non-zero). In general, there is a non-zero probability that photons emitted from a particular activity voxel will interact Chapter 1. Introduction 19 in any other voxel that contains scattering material. The probability that these scattered photons are detected in a given detector bin depends on a large number of parameters including the photon interaction cross sections, the collimator acceptance function, the energy response of the detector, the amount of energy lost in the interaction and the density of electrons in the patient. Scatter correction techniques are rarely employed in iterative SPECT reconstructions due to the complicated nature of scatter distribution, and difficulties, both mathematical and technical that result from non-sparse system matrices. The system matrix mathematically describes a model of the data acquisition process. This model is then used in an iterative reconstruction algorithm to incrementally update the activity distribution. The iterations generally consist of two steps. The first is a forward projection step in which the data collection process is simulated using the system matrix model of data acquisition. A set of estimated projections is obtained using the current iteration of the activity distribution. The next step is often referred to as back projection, where the differences between the estimated projection and measured projection data are used to obtain a new estimate of the activity distribution [9]. MAXIMUM L IKELIHOOD-EXPECTATION MAXIMIZATION Many of the iterative reconstruction algorithms that are in use in the research and clinical environment are based on the Maximum Likelihood-Expectation Maximization (MLEM) algorithm [10]. This algorithm seeks the most likely activity distribution that is consis-tent with the measured projections [9]. The likelihood function is maximized when the difference between the estimated and measured projections is minimized. This method implicitly assumes Poisson distributed statistical noise in the projection measurements. The entire process can be written concisely as: Chapter 1. Introduction 20 i=l k=l where Yi is the measured number of counts in projection bin i, Aâ„˘ is the nth iterative estimate of the activity in object voxel j, and CV, is the system matrix element for activity voxel j and detector pixel i. The sum over M is over all activity voxels in the reconstructed object. The sum over N is over all projection bin measurements. The M term ^ C^A^ is the forward projection step and represents the estimated projection k=l pixel value for detector pixel i. Equation 1.13 represents the back projection step in which a weighted ratio of Yi over the estimated projection pixel value is create a new estimate of the activity distribution. The MLEM algorithm has a number of advantages over analytical methods such as FBP, including the ability to better model the emission and detection processes, improve-ment of image noise properties due to its accurate modeling of the statistical noise in the projection data, and ability to apply constraints on possible solutions based on a-priori knowledge about the activity distribution [9]. One of the disadvantages of MLEM algo-rithm is that it does not converge to a unique solution in the presence of statistical noise in the projection data. The MLEM algorithm will continue to iterate attempting to fit the best distribution to the noise in the projection data, thereby degrading the noise characteristics of the reconstructed image. As a result it is not obvious at what number of iterations the iterative process should be stopped. This optimal number of iterations cannot be known prior to reconstruction because it depends on the spatial frequencies and intensities present in the unknown activity distribution as well as the exact form of the system model used. Chapter 1. Introduction 21 O R D E R E D SUBSETS - E X P E C T A T I O N MAXIMIZATION Another serious limitation of the MLEM algorithm is the time necessary to perform a re-construction. The problem is exacerbated as the system matrix used becomes more com-plicated (less sparse) to better model the physics of the SPECT data acquisition process. A commonly used technique to accelerate MLEM is the Ordered Subsets-Expectation Maximization (OSEM) algorithm [11]. The MLEM algorithm uses all the projection angles in each iterative step of the reconstruction process. The iterative process for OSEM is the same as MLEM with the exception that each iteration is accelerated by only employing a subset of the projection angles progressively using different subsets in each iteration. The OSEM algorithm converges to an image very similar to the MLEM result in a much shorter amount of time [11, 9], provided that the subsets are ordered in such a way that differences between the information contained in subsequent subsets are maximized and the number of projections included in each subset is not too small. Unfortunately, OSEM has the same problems associated with its convergence as MLEM. As described earlier the non-sparse nature of a system matrix that accurately models scatter can increase reconstruction time substantially. One solution to this problem that has been proposed [12, 13], is the use a simplified system model in the back projection step and the complete system model in the forward projection step. This method has been referred to as dual matrix ordered subsets (DMOS) expectation maximization [14]. The differences between images obtained using the same system matrix in the forward and back projection and the DMOS algorithm have been shown to be very small [12, 13, 14, 15], hence all APD based scatter corrections performed in this thesis employ the DMOS reconstruction algorithm. C H A P T E R 2 L I T E R A T U R E R E V I E W Photon scatter can cause significant degradation of image quality and quantitative ac-curacy in reconstructed SPECT images. In general, the detection of scattered photons can lead to reduced image contrast and object-dependent distortion of the image. A few specific examples of scatter-related problems observed in SPECT imaging are the follow-ing: In images which have not been accurately corrected for scatter, the reconstructed count density in a cold sphere near the center of a water cylinder with uniform activ-ity concentration can be as high as 30% of the background activity [16]. In myocardial perfusion imaging photon scatter from nearby organs with high radiotracer uptake can decrease the uniformity of the left ventricle (LV) and thereby confound the detection of perfusion defects [17, 18]. In [18] a simulation study was performed using a phantom in which the liver activity concentration was twice that in the uniformly perfused LV. The reconstructed images that were not accurately corrected for scatter had approximately two times higher count density in the inferior basal region of LV than the anterior basal, due to scatter from the liver. The reduced contrast as a result of photon scatter can also lead to larger absolute errors in the measurement of infarct size in cardiac perfusion studies [19], relative to scatter corrected images. 22 Chapter 2. Literature Review 23 Although the benefits of accurate scatter corrections have long been known, the opti-mal way to implement these corrections is still the topic of a great deal of current research. Many approaches have been proposed, but few of them have been both accurate enough for quantitative imaging and fast enough for routine clinical use. This chapter will begin with a brief overview of some simple scatter compensation techniques, scatter corrections based on the energy spectrum of the acquired data, and some early corrections based on a empirical observations of the spatial distribution of scatter. This discussion is not meant to be a comprehensive review of these topics, since some excellent overviews and comparisons of these early scatter corrections can be found in [3, 20, 21, 22]. Instead, this review will focus on some more recent scatter corrections based on sophisticated theoretical models. 2.1 HOLISTIC COMPENSATION FOR THE EFFECTS OF SCATTER A number of techniques have been proposed that compensate for the effects of scatter without using a direct estimate of the contribution of scattered photons to the measured data. These methods treat the problem of scatter holistically, that is as a modification of another correction (frequently the attenuation correction) or as part of the larger image formation process. Techniques that would be included in the latter category include the use of a modified image-dependent filter in FBP reconstruction (e.g., [23, 24]). These filters are related to the response function of the imaging system and include parameters that are iteratively fit using data acquired in the presence of attenuation, scatter and imperfect detector resolution. These methods compensate for all the image-degrading effects associated with the data acquisition process without distinguishing their origin [20]. Although these filters can improve reconstructed image quality, the effects on quantitative accuracy are Chapter 2. Literature Review 24 difficult to predict since they will depend on type of filter and the specific parameters used. One commonly used holistic scatter compensation technique involves modification of the attenuation coefficients prior to their use in an attenuation correction. As discussed in the introduction, the linear or "narrow-beam" attenuation coefficient can be used to describe the fraction of photons that undergo interactions in an attenuating medium. However, some of the scattered photons may be detected by a SPECT gamma camera. If linear attenuation coefficients are used and no compensation for scatter is made, this leads to an over-correction of the SPECT projection data. It has therefore been proposed that one can use "effective" or "broad-beam" attenuation coefficients that are slightly lower than the narrow-beam values (in the range of 0.12 to 0.14 cm - 1 for 140 keV photons in water [25, 26, 27]), to account for the contribution of these detected scattered photons. A significant problem with this approach is the implicit assumption that the spatial distribution of these scattered photons (i.e., their point of origin in the object) are the same as that of unscattered photons, which is clearly false. Another important difficulty is that the optimal amount by which the narrow-beam value should be changed is case dependent [25, 18]. Although this compensation can qualitatively improve the uniformity of reconstructed images [27], it does not reliably produce quantitatively accurate images [16, 28]. Another approach, based on similar principles, is to replace the entire attenuation correction factor (i.e., exp(â€”ixd)) with a modified function that includes a "build-up" factor to account for scatter. These factors may be determined experimentally [29, 30] or through simulations [31, 32]. The build-up factors are highly dependent on the source position and attenuation distribution within the object being imaged and would therefore be unrealistic to obtain for patient studies. Chapter 2. Literature Review 2 5 2.2 E N E R G Y - B A S E D TECHNIQUES A second class of scatter correction techniques use information about the energy of the detected photons to either reduce the number of scattered photons detected or to estimate the contribution of scattered photons to the acquired data. These techniques use the fact that the majority of scattered photons have undergone Compton interactions and therefore have lower energy than unscattered photons. However, there are two general problems with these approaches which make it impossible to completely discriminate between scattered and unscattered photons using energy information alone. The first is that Rayleigh scattered photons have the same energy as the incident photons. The second is that a large proportion of Compton scattered photons have energies very close to the unscattered energy. The poor energy resolution of the detectors used in SPECT make it impossible to distinguish these scattered photons from the primary photons. Figure 2 . 1 shows an example of a simulated energy spectra for a 9 9 m T c line source located in the center of a water cylinder. This figure is based on data given in [1] . The dotted line represents the contribution from scattered photons and the solid line represents the contribution from both scattered and unscattered photons. 2 . 2 . 1 S INGLE P H O T O P E A K M E T H O D S One simple method to restrict the detection of scattered photons is the use of energy windows, that is to only accept those photons that deposit an amount of energy that falls within a predefined range of values. For example, the typical photopeak window used in 9 9 m T c imaging is a symmetric 2 0 % window centered around 1 4 0 keV (i.e., all photons with energies between 1 2 6 and 1 5 4 keV are accepted). This method is obviously preferable to using no energy window, since it restricts those photons that have a very high probability of having undergone Compton scattering. However, the number of scattered Chapter 2. Literature Review 26 40 60 80 100 120 140 160 Energy (keV) Figure 2.1: Simulated Energy Spectra for a 9 9 m T c line source in a water cylinder, based on data from [1]. The solid line is for all photons, both scattered and unscattered. The dotted line is for only the scattered photons. photons recorded in such a window is still quite high. A simple variation of this method is the use of an asymmetric energy window that is shifted toward higher energies to further reduce the number of scattered photons detected (e.g., [33]). Although this method can reduce the contribution from scattered photons, it also reduces detection efficiency for unscattered photons which degrades image uniformity. 2.2.2 D U A L A N D T R I P L E E N E R G Y W I N D O W M E T H O D S A larger group of energy-based corrections use one or more energy widows to estimate the component of scatter present in the acquired data. One commonly used method is the Dual-Energy Window (DEW) technique (e.g., [16]), in which data is acquired in both the conventional photopeak window and a secondary energy window at lower energies (e.g., from 92 to 125 keV [16]). The distribution in the secondary window, which consists mostly of scattered photons, is then scaled by some constant factor, k, and then subtracted from the data acquired in the photopeak window [34]. Alternatively this data can be used in post-reconstruction [16], in which a "scatter image", reconstructed from the lower energy window data, is scaled and subtracted from an image reconstructed from the photopeak Chapter 2. Literature Review 27 data. The DEW method is based on the assumption that the events recorded in the sec-ondary window have approximately the same spatial distribution as the scattered pho-tons in the photopeak window. However, in reality the photons recorded in the secondary window are more likely to have scattered through larger angles and undergone multiple-scatter interactions than those recorded in the photopeak window. A number of studies (e.g., [35, 36]), have demonstrated substantial differences between the distribution of photons in the photopeak and secondary window. The selection of the scaling factor k also poses some difficulties. The optimal value of k is highly dependent on the geometry of the object being imaged, the energy resolution of the camera, and the criteria used to select its optimal value. For example, k was assigned a constant value of 0.5 in [16], but using a more complex selection criteria in [37] the best value for k was proposed to be between 1.2 and 1.3. Three other methods that also use two energy windows are the dual photopeak window method (DPW) [35] the channel ratio method (CR) [38] and the asymmetric photopeak window (APW) method [39]. All these techniques split the conventional photopeak window into two sub-windows. In the case of DPW and CR, both sub-windows have the same width. While in APW the sub-window widths are asymmetric. All these techniques assume that a relationship exists between the number of events recorded in each sub-window and the total number of events in both windows that are due to scatter. This scatter contribution can then be removed from the combined data in both windows. Since scatter is a complicated function that depends on the object geometry, all of these techniques are approximate and frequently are very sensitive to non-uniformities in the detector response, electronic drift in the camera components, and noise in the data (see [39, 40, 38]). The Triple-Energy Window (TEW) [41] method is a similar Chapter 2. Literature Review 28 technique that employs two narrow windows (e.g., 2 keV in width) on either side of the conventional photopeak window. The shape of the contribution from scattered photons in the photopeak window is estimated by linear interpolation using the data acquired in the two adjacent windows. This technique is also only an approximation (i.e., the true scatter distribution does not vary linearly). The TEW method is very sensitive to statistical fluctuations in the data, due to the extremely narrow widows which are used [41]. 2.2.3 M U L T I P L E E N E R G Y W INDOW M E T H O D S Several more complex methods have been proposed that use more than three energy windows. One such method is the use of an artificial neural network to estimate the scatter spectra from the data collected in five photopeak sub-windows [42]. The network is trained using simulated energy spectra where the truth is known. Another method energy weighted acquisition was first proposed in [43] and later implemented in [44]. The basic principle behind this technique is to apply a different weighting factor to each detected photon that depends on the amount of energy it deposited in the detector. Since lower energy photons are more likely to have undergone Compton scatter the detector position in which they are recorded contains less information about the photons point of origin than that for unscattered photons. The events from lower energy photons are therefore statistically relocated in a 21 pixel matrix centered around the detection location. Although it has long been recognized that the values of these weighting factors should be patient-dependent [43], the values used in practical implementation are not [44]. The local energy spectrum analysis method [45, 46] also collects the data in multiple energy widows and estimates the contribution from scatter based on a particular model of the shape of the scatter spectrum. However, these models may be overly simplified Chapter 2. Literature Review 29 and the optimal choice for the energy range can depend on the specific scatter model used and may vary from pixel to pixel [47]. Another group of methods, that also employ multiple energy windows, involve multi-dimensional analysis of the energy spectra for every pixel. The energy data can be thought of as a matrix, where each row represents the energy spectrum for a given detector pixel and each column can be thought of as a vectorized image of the photons with energy in a particular spectral range. These methods exploit correlations between the spectra acquired for different image pixels to estimate the contribution of scatter. First proposed as a possible method for scatter correction in [48], the method has since been implemented using two different methodologies. The first method, Holospectral imaging [49], does not assume any particular model for the scatter and unscattered spectra. Instead, it determines the variance of the spectral matrix and uses the fact that the component of the spectra due to unscattered photons should exhibit higher variance, since it contains better information about the structure of the object being imaged than the scatter and noise components. The contamination due to scatter is then reduced through filtering. Some specific problems with this method are the large number of energy windows required (e.g., 16 in [49]) and a loss of spatial reso-lution that occurs when certain corrections are applied. These corrections are necessary to account for energy-dependent efficiencies in experimental detectors [50]. The second group of methods use factor analysis, a procedure similar to singular value decomposition (e.g., [51]), in which the energy spectra matrix is represented by a factor space whose bases are found from decomposition of the covariance matrix. The matrix of spectra are first decomposed into a set of orthogonal eigenvectors that are purely math-ematical with no physical interpretation. Secondly, a set of oblique basis vectors, called factors, are determined. These factors represent physical quantities. They may consist Chapter 2. Literature Review 3 0 of a factor for the unscattered spectrum, and one or more for scatter [52].However, this method does not permit proper estimation of the unscattered spectra. Therefore, con-strained factor analysis [54], a variation of this technique, was proposed. In this method conventional factor analysis is used to determine just two factors; the scattered and un-scattered spectra. The factor that more closely resembles the expected spectral shape for unscattered photons is replaced with a theoretical unscattered spectra. Coefficients representing the contribution of the theoretical unscattered factor and original scattered factor are then determined for each pixel. The use of only one factor representing the scatter spectra means the scatter estimate is stationary. That is, it incorrectly assumes that the shape of the scatter spectrum does not change, except by a scaling factor, from one pixel to the next. All of these multiple energy window techniques are sensitive to variations in detector efficiency at different energies [55, 56] and to statistical noise which can be very signif-icant when many narrow energy windows are used. Although these methods are more sophisticated than the dual and triple energy window techniques they are still based on approximations that reduce the accuracy of their scatter correction. They often require modified acquisition equipment or the storage of large amounts of data (corresponding to different energy windows) and therefore may be difficult to implement clinically. 2.3 SPATIAL MODELING OF THE SCATTER DISTRIBUTION A third group of scatter correction techniques attempt to directly model the spatial distribution of both scattered and unscattered photons using an estimate of the source distribution and a projection operation, which maps the activity intensity for each source voxel to a two-dimensional distribution of values for each detector pixel. The scatter correction can then be implemented by directly subtracting the contribution of scattered Chapter 2. Literature Review 31 photons from the measured data. Alternatively one can use this projection operation, that includes the scatter contribution, in the forward-projection step and possibly in the back-projection step of an iterative reconstruction algorithm (see discussion in section 1.7.2). 2.3.1 SPATIALLY INVARIANT S C A T T E R P O I N T S P R E A D FUNCTIONS A simple way to model the projection operation is in terms of a spatial convolution. The intensity value for each voxel in the activity distribution is convolved with a system point spread function (PSF). This PSF can either be decomposed into a scattered and unscattered component or only consists of the scatter contribution in which case it can be though of as scatter point spread function (SPSF). The exact form of the PSF used in this type of correction is usually based on empirical observations from simulations or experiments. Figure 2.2 shows simulated data, based on data given in [1], for an off-axis line source in a water cylinder with the detector in three different positions. Both the contribu-tion from scattered photons (dashed line) and the combined contribution from scat-tered and unscattered photons (solid line) are shown. These data represent line-spread functions (LSF) from which one can determine the system PSF and SPSF. The LSF "wings" (i.e., the portion of the LSF data further away from the peak region) corre-sponds only to scattered photons. The figures are displayed on a semi-logarithmic scale to demonstrate that the scatter wings can be approximated using an exponential func-tion. Hence, two common choices for SPSFs are one-dimensional exponentials and two-dimensional exponentials which are radially symmetric in the detector pixel coordinates (e.g., [57, 58, 59, 36]). Some more complicated spatially invariant functions have also been proposed (e.g., [60, 61]), but the underlying assumption that scatter is spatially Chapter 2. Literature Review 32 invariant is false [62, 63] and the use of an exponential function underestimates the con-tribution due to scatter in the peak region of the SPSF [3]. It is obvious from figure 2.2, that even for the simple object and acquisition geometry modeled in [1], the parameters of the SPSF should depend on the distance between the source and collimator surface. Many studies have demonstrated (e.g., [57, 59, 64]) that the shape of the SPF changes depending on both the object being imaged and the acquisition geometry. Figure 2.2: Simulated spatial projection data for a line source in a water cylinder, based on data from [1]. The solid line is for all photons, both scattered and unscattered. The dotted line is for only the scattered photons. 2.3.2 S LAB D E R I V E D S C A T T E R ESTIMATION A more sophisticated approach is the slab derived scatter estimation (SDSE) method developed in [65, 66, 63, 67, 68] and concurrently in [69, 70]. The SDSE technique also uses a PSF based model for scatter but is more flexible and accurate because it does not assume a spatially invariant PSF. This technique employs a database of scatter distribu-tion functions, experimentally acquired or simulated, for a line source at various distances from the detector surface and at various depths within a homogeneous rectangular water phantom (or "slab" phantom). These data represent camera-specific, object-dependent Chapter 2. Literature Review 33 LSFs. The LSFs are then parameterized by fitting the data with a combination of Gaus-sian and/or exponential functions to account for both the peaked central region and exponential wings of a typical scatter LSF (e.g., see figure 2.2). These measurements can be performed quickly for a wide variety of depth and distance combinations using a tri-angular phantom [70]. The scatter LSF for a more complex object is then approximated by replacing the object's attenuation distribution with an equivalent incremental series of water slabs and appropriately combining the parameterized LSFs from the database. Al-though SDSE was originally derived for homogeneous objects, some work has been done to modify the method for inhomogeneous scattering media [71]. However, the errors in the scatter-to-primary ratio (SPR) were as high as 15% [71] for line sources in a simple inhomogenous water slabs (with air gaps). The SDSE method can also introduce SPR errors greater than 100% [72] for lower energy 2 0 1T1 photons when the source is located close to the edge of a homogeneous scattering media. 2.4 T H E O R E T I C A L L Y - B A S E D SCATTER CORRECTIONS Another approach to model the spatial distribution of scattered photons is to use a theoretically-based description of the interactions of photons with matter (both in the patient and in the detector). Since the formulae that describe these phenomena are well known, one can, in principle, very accurately determine the scatter distribution for a given SPECT detector that corresponds to a particular activity and attenuation distribution. However, the equations that describe the propagation of photons through the scattering medium from the point of emission to the detection location do not in general have analytical solutions. As a result, many of the proposed theoretically-based methods reduce the complexity of the problem, by making some physically-justifiable approximations and then use statistical simulation or numerical evaluation techniques to Chapter 2. Literature Review 34 obtain solutions. One of the major difficulties with theoretically-based scatter corrections are the long computation times required to determine their scatter estimates. It is very difficult to compare these techniques, with respect to their computation times, since these times frequently depend on the size of the source and attenuation distribution, the pixel size, the projection matrix-size, and the specific parameters used in each technique. It is also difficult to compare timing-results from computers with different processor speeds. In the following discussion, calculation times and specific details will be provided when they are available. 2.4.1 E F F E C T I V E S C A T T E R S O U R C E ESTIMATION The problems associated with S D S E lead to the development of the effective scatter source estimation ( E S S E ) method, first derived in [72] and later applied in [73, 74, 75, 76, 77]. This technique attempts to find an "effective scatter source" that depends on an estimate of the object activity distribution, the attenuation map, the projection angle and pre-calculated scatter "kernels". The estimated component of scatter in the acquired data or scatter distribution function ( S D F ) is then obtained by simply projecting the effective scatter source using the same projection operation one would use for the unscattered photons. The E S S E method was intended to be used in an iterative reconstruction algo-rithm, as part of the attenuated projection operation (see the discussion in section 1.7.2). In general, the projection operation including the effects of scatter can be written as: Yi = t;C$G>s\il (2.1) 3 = 1 where Xj is the amount of activity present in the discretized object voxel j. The term Y{ is the number of detected counts in the discretized detector bin i with contributions Chapter 2. Literature Review 35 from both scattered and primary photons. The summation over J is over all voxels in the object containing activity. The system matrix element, C ^ ' G ' 5 , contains geometrical (G) considerations such as the collimator response, the effects of photon attenuation (A), and scatter (S). The difficulties associated with a non-sparse (see section 1.7.2) system matrix that includes the effects of scatter are avoided in ESSE by assuming that equation 2.1 can be rewritten as: * = ÂŁ C * Â° (A,-+ A?) (2.2) j=i where C ^ ' G only contains the geometrical and attenuation effects and the term Aj is an effective scatter source. The effective scatter source can be added to the existing source distribution, then a single attenuated projection can be calculated that provides an estimate of the combined primary and scattered projection data. The attenuated projection, described by equation 2.2, can be performed relatively quickly, but the ESSE method to determine A? requires substantially more time. It is rather difficult to discuss the advantages and problems associated with the ESSE method, as well as modifications that have been made to the technique since its inception, without going into detail on its derivation. It is also interesting to discuss this derivation because of the similarities (and differences) between ESSE and the APDI method which is proposed in this thesis. For simplicity, the ESSE method described here [72] is for a detector with a perfect collimator that only accepts photons perpendicular to the detector surface (i.e., photons with incident direction parallel to a unit vector n perpendicular to the detector surface). The calculation of the ESSE scatter point spread function, SPSF(x, ÂŁ), represents the probability that photons emitted from the source position x are scattered inside the object and are then detected at position t on the surface of the detector. This probability is Chapter 2. Literature Review 36 divided into four parts relative to the position x', the location of the last scattering event in the object before the photon is detected. These parts are; (1) the probability that a photon emitted at x, will reach x' by any possible path (including multiple scatter); (2) the probability that the photon will be scattered at position x' in a direction perpendicular to the detector surface; (3) the probability that the scattered photon will reach the detector without being absorbed; and (4) the probability that the scattered photon will be absorbed in the crystal and not be rejected by the energy discrimination of the detector. The probabilities (1), (3) and (4) can be combined to form the product: kn(x, x')pe(x'), where pe(S') is the electron density at position x' and kn{x,x') is the effective scatter source kernel. Using this the SPSF(s, t) can written as: SPSF(x, t) = / / / pe(x') kh{x, x') exp(-Ai,,a(ÂŁ, X')T(X', h)5(x' - x1 â€˘ hn)dx'. (2.3) where pS)n(x,x') is the scatter attenuation coefficient kernel, which is the appropriately averaged attenuation coefficient for all photons that are emitted from x with final scat-tering position x1. The function r(x', n) is the water-equivalent distance from the point x' to the surface of the detector along h. Finally, the delta function in equation 2.3 selects only those final scattering positions along a line perpendicular to the detector surface at The calculation of equation 2.3 would involve three dimensional integrals which would be different for each object and scattering position considered. Therefore, it is assumed that the effective scatter source and scatter attenuation coefficient are spatially invariant. This is not equivalent to assuming that the total scatter response function is spatially invariant. Instead, this means that all the interactions involved in the probabilities: (1), (3) and (4) are assumed to be independent of the position it occurs inside the object. This is only true for single scatter in a convex uniform object [72]. Therefore, it is position t. Chapter 2. Literature Review 37 assumed that all interactions that may occur as the photon travels from x to x' occur in water. This assumption is made in order that the integral in equation 2.3 can be replaced with a convolution. Hence, SPSF(x, i) becomes: SPSF(x, t) = jjj Peix1) kh(x - x') exp(-fi , i f t (x - f )T(X', ti)5(x' - x' â€˘ nn)dx'. (2.4) Using equation 2.3 one can estimate the total scatter distribution function TSDF(i) for the detector pixel at t. The TSDF(i) represents the total contribution from scattered photons to the detector bin at position t, given a particular activity distribution X(x'): TSDF(i) = UJSPSF(x,i) X{x)dx = Iff (HI P e ( & ) k*(* ~ exp(-At,,ft(x - x1 )T(X', n) x 5(x' - x1 â€˘ nn)dx^j\(x)dx. (2.5) Next the order of integration is interchanged and a term is factored out that represents the primary-photon attenuation probability. This probability factor is exp(â€”p^ T x^' â€” n)), where /io is the attenuation coefficient for water. Then equation 2.5 becomes the standard attenuated projection operation acting on an effective scatter source Xs (x) 1: TSDF(i) = IH A 5(f) exp(-/i 0r(f, h))8(x' - x' â€˘ n^dx1, (2.6) where Xs(x1) is the effective scatter source which can be expressed as: Xs(x1) = IH Pe(z') h(x - x')exp(-AfiStn(x - x1 )r(f, h)X{x)dx. (2.7) 1It should be noted that discretization of the object and detector coordinates, as well as geometrical effects, such as the collimator response, can also be included in the attenuated projection operation. In which case, one would use equation 2.2. Chapter 2. Literature Review 38 The term Aps>fl(x) is the relative scatter attenuation kernel, which is given by: The expression for the ESSE effective scatter source can be expressed as a convolution, if one approximates the exponential term in equation 2.7 using three terms of a Taylor series expansion, then Xs(x') is approximately: where <g> represents a three-dimensional convolution. The steps involved in the ESSE method are described by equations 2.9 and 2.8. Since it is assumed that all interactions that occur as the photon travels from x to x' take place in water, the relative scatter attenuation kernel, Aps^(x) and the effective scatter source kernel kn(x') do not depend on the object being imaged. Both kernels can therefore be calculated using computer simulations of water slab phantoms (hence the connection with SDSE). However, the simulation software must be modified to provide the appropriate values. The expression in equation 2.9 is the object-dependent portion of the ESSE method which depends on the attenuation, electron density, and an estimate of the source distribution. The assumptions made to derive ESSE, particularly the assumed spatial invariance of the scatter kernel, can introduce errors in SPR values between 4% to 5% [72] when sources are located close to the edges of an attenuating media. It is also hard to predict how errors these accumulate in different imaging situations. For example, the assumption that all interactions between x to x' occur in water may cause significant errors for highly non-uniform attenuation maps, such as the chest (e.g., coronary perfusion studies). The A/xs,ft(f) = Hsfiix) - po- (2.8) p e ( f ) A(ÂŁ) Â®kh(x') - A(ÂŁ) Â® (kn(x')ApsA(x'))T(x',h) + (2.9) Chapter 2. Literature Review 39 chest contains a large amount of lung and bone tissue, which both have significantly-different attenuation coefficients from that of water. Since the ESSE requires either three image space convolutions or several Fourier transforms followed by frequency domain multiplications, it can be extremely computa-tionally intensive [72, 75]. Using a 64x64x24 activity and attenuation distributions, the ESSE method requires about 7 minutes [75] to generate 64 scatter projections each with a matrix-size of 64x64, on a Sun Ultra computer with dual processors (each 200Mz). This is considered too long for standard clinical reconstructions. Therefore, some effort has been made to accelerate the algorithm [75]. One strategy considered was the use of a larger voxel size (hence smaller matrix size) in equation 2.9, and then expanding back to the full matrix size using interpolation prior to the attenuated projection step (equation 2.6). This method can reduce computation time significantly (from 7 minutes to 45 seconds for the study and computer described above), but leads to an underestimation of the peak regions of the scatter projections which can introduce artifacts in scatter corrected images [75]. Although this method can reduce computation time it is difficult to predict how error will accumulate when a new set of approximations are combined with the existing set of ESSE approximations. Another very promising scatter correction technique which is also uses an effective scatter source but is not based on the ESSE method has recently been described in [78, 79]. This method models the propagation of scattered photons as a slice-by-slice blurring. The parameters of the blurring kernel are calculated "on the fly" and depend on the Klein-Nishina (KN) scattering cross-section, but unlike ESSE this technique can accurately model the propagation of scattered photons through a non-uniform medium from their emission point to their final scattering point. The scattering angle in the Klein-Nishina formula is expressed in terms of the point of emission, the scattering point, and the Chapter 2. Literature Review 40 detection point. A spatial scatter response function (SSRF) is used to propagate scattered photons, starting from the plane containing the emission point, through each scattering slice. At each slice an intermediate scatter source image (ISSI) is convolved with the SSRF scatter kernel, which is based on the KN formula and a projection operation that uses both attenuation and geometric detector response effects. The advantages of this technique include the following: its scatter kernel does not require the generation of complicated pre-calculated tables; it includes the effects of a non-uniform scattering media; and it uses an accurate KN model of scatter. Its disadvantages include fairly long computation times. For example, when an accurate SSRF is used that varies as a function of distance between the source and scattering slice, the method requires 35 minutes to calculate a set of 64 first-order Compton scatter projections with a 64x64 projection matrix-size [78] using a 64x64x64 matrix size for the attenuation and activity maps and a 167 MHz processor. The computation time required to generate the second order Compton scatter distribution is about twice that for the first order calculation [79]. Another problem with this method is its limited flexibility. It can only be applied for the scatter photopeak energy. For example, consider a scan which is performed using 140 keV 99mTc emission source and a lower energy 68.9 to 82.6 keV 2 0 1T1 transmission source to measure the attenuation distribution. Some scattered emission photons will lose enough energy to be detected in the transmission energy window. This distribution of scattered photons cannot be estimated using a direct application of this model [78]. 2.5 MONTE CARLO-BASED SCATTER CORRECTIONS Another way to model the spatial distribution of scattered photons is with the use of Monte Carlo (MC) simulations. These simulators stochastically propagate individual Chapter 2. Literature Review 41 photons through a scattering media and record their interactions. MC is a statistical simulation method that utilize sequences of random numbers to solve problems or model physical systems, which are dependent on particular probability distribution functions (PDFs). These random numbers are used to sample the PDFs that describe the behavior of the system. Due to the statistical nature of MC methods, typically many individual simulations must be performed, often called "trials" or "histories", to obtain a result that accurately reflects the behavior of the system under investigation. In many practical applications one can predict the statistical error ("variance") associated with a given result, and one can therefore estimate the number of Monte Carlo trials that are needed to achieve a given error [22]. Since simulated MC projection data can be easily separated into scattered and un-scattered components, scatter corrections are an obvious application for MC simulations. Some early investigators [58] concluded that accurate MC simulations were too time con-suming for clinical use. A major difficulty with MC methods is that when the "true" probability functions are used the simulation may require an unacceptably long compu-tation times. For example, the geometrical efficiency of a SPECT collimator is on the order of 10 - 4. If an accurate collimator response PDF was used, it would mean that only 1 in every 10,000 simulated photons would pass through the collimator. Therefore, to improve MC performance one must either use more approximate PDFs, degrading the accuracy of the result, and/or employ some form of variance-reduction (VR) technique (also called importance sampling) [80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90]. These meth-ods force the photon to travel along a path in which it is more likely to be detected. The photon is then given a "weight" between one and zero, which represents the probability that it would have traveled along that path. If no VR techniques are used, then each photon would contribute a value of one or zero to the final image depending on whether Chapter 2. Literature Review 42 or not it is detected. Three VR methods commonly used in SPECT are [82]; forced non-absorption, forced detection, and stratification. Forced non-absorption restricts the allowed types of photon interactions to increase the likelihood of detection. For example, one can specify that no scattered photon can lose so much energy that it would no longer be recorded in the energy window of the detector. Similarly the photoelectric effect can be forbidden preventing the complete absorption of photons. Forced detection further restricts the allowed photon interactions by requiring that all photons scatter in a direction that is within the acceptance angle of the collimator. When a photon reaches a site where a number of different interactions are possible the photon history is split. As described in [90], at each possible interaction site a copy of the original photon is created. This copy is forced to escape from the object in a direction that allows it to pass through the collimator and be absorbed in the detector crystal. This direction defines a scattering angle. The copied photon weight is then corrected for the probability of Compton scatter at this angle (using the KN formula), the probability that it is not photo-electrically absorbed, and the probability it will not be attenuated on the path from the interaction site to the detector. This photon weight is then added to the projection bin where it was detected. Then starting from the initial interaction site the original photon is allowed a new interaction, which defines a new direction of propagation. The original photon propagates until it reaches a new interaction site, where another copy is created and forced to interact in such a way to allow its detection. This process is repeated until the original photon escapes the object. In stratification [82] the initial starting position of the photon is not selected at random. Rather the locations that have a higher probability of resulting in photon detection are selected with greater frequency. Stratification can only be be done after Chapter 2. Literature Review 4 3 an initial simulation is performed, generally with fewer photon histories. The results of this initial simulation are then used to determine which initial starting points are more likely to result in detected photons, and therefore should be sampled more frequently in the final simulation. However, MC simulations that employ these VR techniques still require extremely long computation times for large source distributions. Recently convolution-based forced detection (CFD) [91, 92, 88], a method that combines variance reduction with an analyt-ical description of the geometrical collimator response, has permitted the reconsideration of MC for scatter corrections. In standard forced detection the path of a scattered pho-ton is randomly selected from a restricted set of possibilities that would still allow it to pass through the collimator. In CFD each scattered photon is forced to travel in a direction perpendicular to the camera surface and then weighted appropriately. The photon weights are stored in a three-dimensional sub-matrix. In addition to the nor-mal two-dimensional set of detector coordinates, a third coordinate is recorded for each photon, which represents the distance between the detector plane and the plane where the final scatting event took place. After all the photon histories have been calculated, each layer of the sub-matrix is convolved with a distant-dependent collimator response function. Finally, the convolved data is added layer-by-layer to create the final projection data. This method can create smoother and more uniform-looking scatter projections when compared with MC data using standard VR techniques and the same number of photon histories. Since in CFD every photon is projected as a blurring kernel, if too few photon histories are simulated the projections demonstrate correlated noise-textures or "lumpiness" [88]. Some distinct advantages of the CFD technique are: it can be used to very accu-rately account for septal penetration and contamination from X-rays emitted from lead Chapter 2. Literature Review 44 collimators [89]; it can easily be applied to both emission photons [88] and photons that have down-scattered into lower energy windows [92, 87]; and it only requires very simple modifications to any existing MC code. However, the number of histories required to reduce the "lumpiness" effect still make computation times too long for practical clinical use. A single scatter projection calculation requires about 30 minutes, using a computer with a 200 MHz processor, 128x128x128 attenuation and activity maps, a 128x128 pro-jection matrix-size. These calculations simulated 5 orders of scatter and used 107 photon histories. If a 64x64x64 sized object and 64x64 sized matrix is used this time is reduced to 2.3 minutes. Some work has been done to accelerate this method by calculating projection data only for a few projection angles and then using interpolative approximations to estimate the data for the remaining angles [87]. This can decrease computation times substantially. However, this approximation introduces errors to the scatter estimates and applies only to down-scattered photons. These approximations cannot be generalized for photopeak scatter [87]. 2.6 COMPLETE ANALYTICAL SCATTER CORRECTIONS The final category of scatter corrections to be discussed here, are those based on an analytical description of photon scatter. These methods attempt to numerically solve the equations that describe the propagation and interactions of photons from their point of emission through a scattering medium to their point of detection. There are two distinct advantages to this type of technique. The first is that they can model any system with arbitrary accuracy, limited only by the accuracy of the models it uses and the precision of its numerical calculations. The second advantage is the complete absence of stochastic noise in its scatter estimates. This is an extremely favorable property Chapter 2. Literature Review 45 when performing scatter corrections in iterative reconstructions algorithms (e.g., [93, 88]). Two techniques of this type, which use an accurate model of photon scatter and are applicable to non-uniform media, are those first described in [94] and [95]. The first method solves the analytical photon propagation equations "on the fly" using direct numerical integration. The approach is capable of accurately determining all orders of Compton scatter but requires extremely long computation times for anything higher than first order scatter. The calculation of a single projection requires between 16 to 58 hours [94] if one includes primary, first and second order Compton scatter. These calculations were performed using a Sun Sparc 2 computer (no processor speed provided), a small extended source distribution, and a 128 x 128 x 128 attenuation matrix size, and a 128x128 projection matrix size, The second order scatter calculation contributes the majority this computation time but must be considered because it significantly contributes to the detected scatter in SPECT photopeak windows (e.g., [1, 96, 97]). Therefore, this method was not considered practical for accurate scatter correction. The second method is the analytic photon distribution (APD) method [95]; the scatter correction on which the work in this thesis is based. A detailed description of the APD method, as it was originally developed by Wells et al., can be found in chapter 3 section 3.1. 2.7 CONCLUSION Although many scatter correction techniques have been proposed, none of them possess the required accuracy and computational speed that would allow them to be implemented clinically. A number of promising new techniques have been introduced that are based on sound-physics and do not overly simplify the problem of photon scatter in SPECT imaging. Of these the accelerated MC technique, CFD, is probably one of the best methods. The CFD scatter calculation is theoretically-based, uses an accurate physical Chapter 2. Literature Review 46 model, is extremely flexible, and would not be difficult to implement. However, as with all MC techniques it requires long computation times because many photon histories must be generated to obtain low statistical noise. The APDI method proposed in this thesis is also theoretically-based. It uses APD, an analytical description of photon scatter, which is inherently noiseless and extremely accurate. It is hoped that the APDI method can be made fast enough for clinical applications. To achieve this goal, APDI makes some simple but physically reasonable approximations in the context of a SPSF description of scatter. C H A P T E R 3 T H E M E T H O D The following chapter1 describes the APDI scatter estimation method and OSEM-APDI, its implementation as a scatter correction in the OSEM reconstruction algorithm. This chapter is divided into three sections. The first section is an overview of the APD method as it was originally implemented by Wells et al. [3]. APD is an analytical method to determine the spatial distribution for both primary and scattered photons in SPECT projections. This section focuses on the expressions for the distribution of primary photons, first and second order Compton scattered photons, first order Rayleigh scattered photons and approximations made for higher order scattered photons. Also included is a discussion of some of the approxima-tions which were made in the derivation of these expressions. The second section describes APDI, a new scatter estimation method, that determines approximate expressions for scatter point spread functions (SPSFs) using the full APD scatter distributions in a tri-linear interpolation scheme. This approach can substantially reduce computation time and relies on the assumption that the spatial characteristics of *A version of the second and third sections in this chapter have been submitted for publication. E. Vandervoort, A. Celler, R. G. Wells, S. Blinder, K. Dixon and Y. Pang. Implementation of an Analytically Based Scatter Correction in SPECT Reconstructions. IEEE Trans. Nucl. Sci., 2003. 47 Chapter 3. The Method 48 scatter are slowly varying. This section also discusses the validity of the approximations made to derive the APDI method. The final section describes OSEM-APDI, the implementation of the APDI scatter estimation method in the iterative OSEM reconstruction algorithm. Also included is a description of other reconstruction techniques, which were used as reference to evaluate the performance of the OSEM-APDI reconstruction. 3.1 A P D METHOD The APD method numerically evaluates the probability that a photon emitted from a particular source position will undergo Compton or Rayleigh interactions in an attenu-ating medium and be detected by the gamma camera. This calculation involves multi-dimensional integrals of the analytical equations that describe the interaction, propaga-tion and detection of photons in a SPECT imaging system. An important difference between the APD method and other theoretically-based scat-ter estimation techniques is its separation of this calculation into two components. The first component depends on parameters which are independent of the patient, such as the detector and image space geometry, the energy of the emitted photon and the selected energy window. This calculation needs only to be performed once and is then stored as a look-up table for a particular SPECT camera system. The second component depends on patient specific parameters such as the distribution of activity and attenuating material in the patient body. 3 . 1 . 1 PRIMARY P H O T O N D ISTRIBUTION The APD primary or unscattered photon distribution, although not directly used in our scatter correction, plays an important role in OSEM-APDI (see section 3 . 3 ) . Prior to the Chapter 3. The Method 49 s Figure 3.3: Schematic diagram of three different photon paths. The primary or unscattered photon is shown as the dash-dotted line. The primary photon travels directly from the source position: s, to the detector at n. The first order scattered photon, shown as the solid line, originates at source position s, and scatters once at t, and is then detected at n. The second order scattered photon (dashed line) scatters twice, once at position u and again at t, and is then detected at n. reconstruction, the sum of the two distributions for scattered and unscattered photons, calculated using the APDI method, is rescaled to match the experimental projection data. A schematic diagram of the path of an unscattered photon is shown as the dash-dotted line in figure 3.3. In this discussion s, i, and u will refer to continuous vector variables confined to the volume of the source or scattering voxel, while n represents a vector confined to the area of the detector pixel. The vectors s"j, tk, um and n*; are coordinates representing the center of the voxel j, k, m and the center of the detector pixel i, respectively. The APD primary photon distribution PPD(A,nj), for detector bin i can be written in terms of the three-dimensional distribution of activity A as: cos(Q (3.10) where: â€˘ Xj is the source strength (the total number of photons emitted from the source Chapter 3. The Method 50 voxel j), â€˘ F(ÂŁ) is the probability that a photon reaching the collimator surface at an angle ÂŁ with respect to the detector normal will pass through the collimator and be detected, â€˘ PB(0) is the probability that an unscattered photon will be recorded within a 20% energy window centered around the energy of the emission photons, â€˘ rnSj is the distance between the vectors ft and Sj, cost c) â€˘ â€”-â€”d2n is the solid angle of a sphere centered at Sj and subtended by the detector pixel element d 2n and â€˘ p,(x) is the attenuation as a function of distance along the path traveled by the photon. The summation over J is over all voxels containing activity. The expression PJS(0) is a special case of PJE(0), the probability that a photon that has Compton scattered with angle 9 will be detected in the energy window of the detector. The integration in equation 3.10 is performed with respect to the continuous variable n over the area A; of a pixel centered at the position n*. This calculation assumes fixed source position Sj (i.e., the source position it is not a variable over which the integration is performed). To speed up calculation times, the attenuation path of the photons is assumed to be directly from the points Sj to n*j. Certain choices must be made about the exact functional form of the expressions in equation 3.10. For all the work that was done in this thesis, the following models for F(ÂŁ) and P.E(0) were used. Chapter 3. The Method 51 The collimator acceptance function, F ( ÂŁ ) , was modeled as a parallel hole collimator with cylindrical holes and no septal penetration. The functional form of F(f) used was similar to that described in [94, 8] and accurately models the distance-dependent spatial resolution of a S P E C T gamma camera. The energy response of the detector was modeled as a Gaussian with FWHM equal to 10% of the incident photon energy and a 20% energy window was used in the detection process. This form of the energy response is the same as that described in [98]. 3.1.2 F IRST O R D E R C O M P T O N S C A T T E R In figure 3.3 the first order scattered photon is shown as the solid line. The photon originates at the source position s and scatters once at t with Compton scattering angle 9. The photon is then detected at the position n on the detector surface. The A P D first order Compton scatter distribution function S D F ^ ( A , nj), for detector pixel i can be written as: S D F W ( A , fk) = ÂŁ Xj ÂŁ Pe(tk)KÂ§ini exp (- f" p(x)dx - T E(9))dx) , (3.11) j=l k=l V Js"i J t k I where: â€˘ E{9) is the energy of an emission photon that has Compton scattered with angle 9 (equation 1.4), â€˘ ~^:{9,a) is the Klein Nishina scattering cross section (equation 1.3), du â€˘ a is the scattered photon energy E(9), divided by the rest mass of an electron, â€˘ p(x, E(9)) is the attenuation coefficient at position x and photon energy E(9) and Chapter 3. The Method 52 â€˘ pe(tk) is average density of electrons for the voxel at tk-The summations over J and K are over all voxels containing activity and over all vox-els containing scattering material, respectively. The remaining variables are defined in equation 3.10. The factor K^ f c â€ž. is independent of the patient and is precalculated and stored as a look-up table. This calculation also assumes fixed source position ij. The volume integration is performed over all possible scattering positions t within the volume V*. of the voxel centered at t^. The area integration is performed over all positions n within the pixel area Aj centered at position n*. The size of this look-up table could potentially be extremely large if it was necessary to calculate the look-up table element for every combination of source position, scattering site and detector bin. Fortunately, the removal of the patient dependencies allows parameterization of the look-up table into only four distances and one angle. This parameterization is described in figure 3.4. It should be noted that the Klein Nishina scattering cross section (equation 1.3) includes the term S(x, Z), which depends on Z a property of the distribution of matter in the patient. Since, these Z values are unknown at the time of precalculation of the look-up table, the S(x, Z) values for water at the unscattered photon energy are used. The overall error introduced by this approximation to the acquired photon distributions has been estimated to be less than 0.2% [3] for 9 9 m T c photons. It should be noted that maps of pe{tk) values are not usually available in SPECT studies. These values were therefore estimated from the attenuation maps. Linear inter-polation was used to obtain electron densities for each object voxel from the relationship between p and pe for water, lung and bone as given in [101]. To shorten calculation times, an approximation was applied (developed by Glenn Wells after the original publications of the APD method) to estimate the variation of Chapter 3. The Method 53 Figure 3.4: Diagram of the parameterization of the APD look-up tables. The parameters are defined with reference to a plane parallel to the detector surface, shown as the x-y plane in the diagram. The vectors s, t, and n are the positions of the source voxel, the last scattering voxel and the detector pixel, respectively. The detector normal is parallel to the z-axis. The APD parameters are four lengths (a-d) and one angle (e). They are shown as dashed lines in the diagram. The parameters (a) and (b) are the distances from the center of the detector pixel to the center of the final scattering voxel measured perpendicular (a) and parallel (b) to the detector normal. The parameters (c) and (d) are the distance from the center of the source voxel to the center of the final scattering voxel, measured parallel (c) and perpendicular (d) to the detector normal. The parameter (e) is the angle between the vector from s to t and the vector from i to n projected into the x-y plane. Chapter 3. The Method 54 linear attenuation with energy. The attenuation maps used in APD are for the energy corresponding to the emission photons. However, scattered photons have lower energy. Therefore, the p(x,E(9)) values for lower energies had to be appropriately rescaled. Prior to the APD calculation, the attenuation map is read into a separate program that determines the unscattered-photon attenuation-lengths p0ÂŁ, which are given by: where t is the vector coordinate of every voxel that contains scattering material, n is the location of every detector pixel where the photon could be detected, and p(x, E(0)) is the attenuation as a function of distance for the unscattered photon energy, E(0). The number of detector pixels which require these attenuation lengths to be calculated are limited by the acceptance angle of the collimator. This data is stored as an "evaluated at-tenuation map". When these values are later used in the APD calculation for scattered photons with energy E(9) the following approximation is made: where pw(E(0)) and pw(E{9)) are the linear attenuation coefficient for water at the unscattered and scattered photon energy, respectively. For most biological materials, the error made by this approximation should be extremely small. 3.1.3 S ECOND O R D E R C O M P T O N S C A T T E R A possible photon path for a second order scattered photon is shown as the dashed line in figure 3.3. The photon scatters at position u with Compton scattering angle tp, then scatters again at t with Compton scattering angle (p. The photon is then detected at n. (3.13) (3.14) Chapter 3. The Method 55 The APD second order Compton scatter distribution function SDF^(A, nj), for de-tector pixel i can be written as: SDF(2>(A,nj) = J^Xjf^p^K^t fk p(x)dx] x j = i k=i \Js"j J exp (- j** p(x)dx - p(x, E(d))dx^J , (3.15) where: exp (p-w(vSjtk - rUSj - r u t j ) . (3.16) The expression for the second order look-up table (equation 3.16) depends on K^? i / t n., the first order look-up table (equation 3.12), with the voxel position Sj replaced with Uj. The terms fxw and pw are the linear attenuation coefficient (at the unscattered photon energy) and the electron density of water. The expression a 0 is the unscattered photon energy E(9 = 0Â°) divided by the rest mass of an electron. The term PE('4>,(1)) is the probability that a photon that has scattered twice, through angles ip a n d </>, will be detected in the energy window of the detector. The summation over K' is over all voxels in the object. All other variables are the same as those defined in equations 3.10, 3.11 (2) and 3.12. An important feature to note is that <fcn. can be parameterized in the same way as K^| f c n . and is therefore not any larger than the first order look-up table. Unlike the analytical method described by Riauka et al. [94], the APD second order scattered photon distribution (equation 3.15) takes approximately the same amount of time to calculate as the first order distribution (equation 3.11). In order to write the second order calculation in this form, a number of approxima-tions are made in the APD method. A more exact form of the second order scatter Chapter 3. The Method 56 distribution would requires two summations over all potential scattering sites. The first summation would be over all possible first scattering positions Uj and the second would be over all final scattering positions t\.- Both of these summations should be over all voxels containing scattering material (i.e., the summation in equation 3.16, should run from 1 to K, rather than from 1 to K'). In the APD method it is assumed that all potential first scattering sites contribute to second order scatter rather than just those containing scattering material. Since the distribution of attenuating material is unknown at the time of precalculation, it is also further assumed that the attenuation along the path from the point of emission to the first scattering site is the same as that for water. Another significant approximation is made in the calculation of the Klein-Nishina scat-tering cross sections. It is assumed that the energy of the photon is unchanged after the first Compton interaction. The errors introduced by these approximations are minimal. For greater detail on these and further approximations made in the second order scatter calculation, the interested reader is directed to [2, 3, 102]. In general, small errors in the second order scatter distributions will have little effect in scatter corrections since second order accounts for less than 12% of the scattered photons recorded in a 20% symmetric photopeak window for 9 9 m T c [1, 96]. 3.1.4 H IGHER O R D E R C O M P T O N S C A T T E R In 9 9 m T c imaging higher order Compton scatter (greater than second order) does not contribute significantly to the photopeak data. For example, simulations of a small 9 9 m T c line source deep in a homogeneous water phantom [96, 1] show that the Compton scattered photons detected in the photopeak window consist of 87% first order, 12% second order and higher order contributes only 1%. Therefore, the contribution from higher order scatter is approximated by rescaling the second order scattered photon Chapter 3. The Method 57 distribution by 13/12. 3.1.5 R A Y L E I G H S C A T T E R The APD first order Rayleigh scatter distribution function R D F ^ ( A , n*j), is given by: RDF^CA, nj) = ÂŁ Xi ÂŁ Pe(4)Kf^ ni exp (- f" p(x)dx - H p(x, E(6))dx) , (3.17) j=l fc=l V Jgi Jt>' I where: A A PP R O X c1 In equation 3.18 â€”^â€”(8) is an approximate from of equation 1.1, which allows the contribution from Rayleigh scatter to be accounted for with a linear addition to the first order Compton scatter look-up tables. In this way, first order Compton scattered photons and first order Rayleigh scattered photons can be calculated simultaneously by replacing K $ f c n i in equation 3.11 with K $ f c n , + K ^ . The approximate form of the Rayleigh scattering cross section is given by: ^ f p * ) = f (l + cos20) [ n X 2 w Z w ) ) 2 Zw, (3.19) where Zw is the effective atomic number of water, T(xik, Zw), is the Rayleigh scattering form factor for water at the unscattered photon energy and momentum transfer Xik of the scattered photon. A number of approximations are made in the derivation of the APD method to facilitate the combination of the Rayleigh and Compton scatter first order look-up tables. They are expected to have minimal influence on overall the scatter distribution function. More details on these approximations and their significance can be found in [3, 103, 104]. Chapter 3. The Method 58 3.2 THE A P D INTERPOLATION ( A P D I ) METHOD Although the APD method is highly accurate the calculation times required for even a limited sized source map are extremely long. For example, using a small extended source distribution (2388 non-zero voxels) with a 64x64x64 matrix size for the activity and attenuation maps, the APD calculation takes between 50 to 90 minutes per projec-tion angle (866 MHz processor, 1G SDRAM). The time variations are because the APD method requires different computation times depending on the depth of the source distri-bution within the attenuating medium. The time necessary to perform an APD scatter correction would be about 75 hours, for a scan with 64 projections acquired for this sim-ple activity distribution. Since iterative image reconstruction incrementally changes the activity distribution, one might wish to update the scatter correction a number of times during the reconstruction process. This would further increase the computation time. It was therefore beneficial to make a number of approximations to speed up the calculation. The APD method can be used to calculate the scatter point spread function (SPSF) for individual source voxels positioned in a particular attenuation map. The SPSF is a two dimensional distribution of probabilities that photons emitted from a particular source voxel will undergo Compton or Rayleigh scatter interactions in the attenuating medium and are then recorded in the array of detector bins. The APD formula calculating the SPSF for the detector pixel centered around the vector coordinate Hi and the source voxel centered around Sj is given by: exp I- j* p(x)dx - fi(x, E(6))dx\ . (3.20) Chapter 3. The Method 59 In simulation studies, later described in this thesis, only first and second order Comp-ton scatter were considered because the Monte Carlo program which was used to create the data does not model Rayleigh scatter and because higher order Compton scatter does not contribute significantly to 9 9 m T c photopeak data. In this case, the APD formula for the SPSF simplifies to: SPSF(55,ni) = ^pe(tk)(K^tkni+K^tknft\(x)dx) x k=l \ Js"i ) exp ^ - J*k p(x)dx - p(x, ÂŁ(0))dxj . (3.21) The SPSFs defined by equations 3.20 and 3.21, can be used to determine the total distribution of detected scattered photons given a particular activity map. The APD total distribution of scattered photons is obtained by summing the product of the ac-tivity intensity and the appropriate SPSF over all voxels containing activity. The APD expression for the total scatter distribution function T S D F / â€ž H ( A , Hi), for the detector pixel at rii is given by: j T S D F / U â€ž ( A , Hi) = ÂŁ A jSPSF(5,-, Hi), (3.22) J ' = I where SPSF (<fj, Hi) is defined by equation 3.20 or 3.21. The summation over J is over all voxels containing activity and Aj is the total number of photons emitted from activity voxel j. When an analysis was performed using a digital phantom of the distribution of atten-uating medium in the thorax, very little variation was observed in the APD-calculated SPSFs for individual activity voxels that were separated by small spatial distances. Therefore, in order to speed up computation times the APD interpolation (APDI) method was developed. The procedure for the APDI can be described in three steps. First, the Chapter 3. The Method 60 full A P D calculation is performed for a smaller subset of activity voxels arranged in an evenly spaced cubic grid. Secondly, the SPSF values for all remaining voxels are esti-mated using tri-linear interpolation. At this step two approximate corrections are made. Preferential photon scattering at low angles results in peaked scatter projections and so in order to improve the A P D I scatter estimate the sampled A P D probability distributions are shifted. The interpolated SPSFs are also approximately corrected for the amount of attenuating material that the scattered photons have passed through. The final step is to use the A P D I approximate SPSFs to determine the estimated total distribution of scattered photons for a particular activity map (as in equation 3 . 2 2 ) . 3 . 2 . 1 D E F I N I N G T H E A C T I V I T Y M A P S U B S E T In the A P D I method the two dimensional SPSF(s}, ni) matrices are calculated using the original A P D method only for activity voxels on the corners of a three dimensional cubic grid. The spacing between the corners of this cubic grid could be set to any number of pixel lengths. However, for the work done in this thesis the size of the grid used was always 4x4x4. This value was selected because it gave the best compromise between the decrease in calculation time and the loss of accuracy associated with using a larger grid size. Smaller grid spacings were not investigated because the calculation times would still be too long for implementation in image reconstruction. Larger grid spacing would introduce unacceptably high errors to the interpolated SPSF values. To save time the SPSFs are only calculated for regions of the grid that contain activity. Figure 3 . 5 shows a transaxial slice through a digital phantom representing an activity distribution in the heart. Overlaid upon the activity distribution in black is the interpolation grid with 4x4x4 spacing. Chapter 3. The Method 61 a a â€˘ a a a a a a a a a i i i i i i 10 20 30 40 50 60 Figure 3.5: Transaxial slice through a digital activity distribution for the heart. Activity intensity values are shown on a linear scale from white to gray. Overlaid upon the activity distribution in black is the interpolation grid with 4 x 4 x 4 spacing. 3 . 2 . 2 E S T I M A T I N G T H E S P S F s U S I N G T R I - L I N E A R I N T E R P O L A T I O N The procedure for tri-linear interpolation that was used in A P D I is extremely simple and was derived from [51]. A more complicated interpolation scheme could easily be implemented. For example, one could use a method, such as tri-cubic spline interpola-tion [106, 51], that ensures higher order smoothness by estimating gradients. However, these methods would increase calculation time. A more complex interpolation method is probably not necessary because the equations that describe the propagation of scat-tered photons in the full A P D method are generally quite smooth, with the exception of discontinuities between voxels belonging to different tissue types. As a result projections calculated using A P D or any linear combination of A P D projections (such as those used in a tri-linear interpolation scheme) are generally quite smooth. 20 30 h 40 ^ 50 Chapter 3. The Method 62 In the following discussion it is convenient to define a separate coordinate system for the activity voxels and detector pixels. In figure 3.6 the activity voxel coordinates are shown as the x-y-z axes and the detector coordinates are shown as the y'-z' axes. A n activity voxel centered at s has Cartesian coordinates (sx,sy,sz) and the detector pixel centered at n has coordinates (nyi,nzi). In the A P D calculation the position of the de-tector surface is fixed in a particular y-z plane relative to the activity voxel coordinate system. We can therefore define the y'-z' plane as x â€” 0 in the activity voxel coordinate system. Hence, the detector position n has coordinates (ny>,nzi) in the y'-z' detector co-ordinate system and coordinates (0,nyi,nzi) in the x-y-z activity voxel coordinate system. a Detector I Plane Figure 3.6: Schematic diagram of to illustrate tri-linear interpolation as it was used in the A P D I method to estimate the SPSF, represented here as f(s), for an activity voxel centered at s with spatial coordinates sx,sy, and sz. Points for which the value of the two dimensional SPSF(s}, rf,) array has been determined are labeled as open circles. Points for which the function is estimated using linear interpolation are shown as solid circles. The detector position is fixed in the y-z plane at x=0. The scatter point spread function SPSF(s , rf;) can be thought of as a function that depends on the spatial coordinates of the activity voxel s. Consider the function f which depends on the spatial coordinates x,y, and z. The value of f(s) = i(sx,sy,sz) can be estimated using linear interpolation i f the value of f is known for f(s"0), / ( s i ) , / (s^), Chapter 3. The Method 63 (3.23) Given that: /(ss), f{s~l), f(s*5), /(sJ), and /(Â«}) , where: / ( 5 o ) = f(sx0i s y o > 5 z o ) ' = / ( s m svo> s z o ) > / ( S _ 2 ) = / ( s x i , Sj / i , S 2 o ) , / ( S 3 ) = / ( S i 0 , S y i , S 2 Q ) , f{si) = f(sxo> syo> szi)' f(s$) ~ fisxi> s y o > szi)> f(sh) â€” f(Sxi, Syi, s 2 i ) i / ( s _ 7 ) â€” f(sx0, S y i , SZ1). $xo Sx <i SXl, Sy0 < <C S y x , S Z o <C S 2 <C SZl $xi $XQ ~ $yi Syo â€” Â§zi $zo â€” & sx &xo â€” b, Sy Syo â€” ^ and Sz SZQ â€” d, we can use linear interpolation to estimate the values of the function! at the points go, gl, g2 and g3 as: 00 = f(sx,syo,szo) Â« (1 - ÂŁ) / (s5) + ÂŁ1 = f(sx,sVl,sZ0) Â« (1 - J)/(si) + 52 = / ( s x , S y o , S Z l ) Â« (1 - + J/(se), 03 = f(sx,SyiiSXl) Â« (1 - f ) / ( 5 4 ) + J / O S ) . Similarly, = / ( Â« * , Sâ€ž , S z o ) Â« (1 - f )0O + ^01, / i l = / ( S x > Sy, SZ1) Â« (1 - f )03 + ^02, and /(*) = /(ss> sV) sz) w (1 - - ) / i â€ž + - / i i (3.24) o o We can express equation 3.24 in terms of the known values f(s~o), /(si), / ( S 2 ) , / ( S 3 ) , / ( S 4 ) . /(ss), /(s"e) and f(s"7) as: 7 /(s) = f(sx, sy, sz) Â« ÂŁ w m / (4) , (3.25) m=0 where: Chapter 3. The Method 64 fo = ( i - ÂŁ - i - d _ i _ be bed , a ^ a2 a 3 " t " 6d , c d \ a 2 a 2 J ' = (" i feed a 3 be a2 6d\ a2' w2 = / be bed \ \a2 a 3 / ' w 3 = ( ' i 6cd + a 3 be a 2 cd\ a2> w* = (i + 6cd 6d cd \ a 3 a 2 o 2 / ' w5 _ /6d â€” U2 bcd\ a 3 )â€˘> w6 (bcd\ ~ V a 3 >' / cd bcd\ a 3 / ' (3.26) The coefficients of interpolation, wm in equation 3.26, only depend on the size of the interpolation grid and can be stored and reused for a given grid size. In this case, tri-linear interpolation is to be used to obtain an approximate expression for SPSF(s, ni), using a subset of APD-calculated SPSFs (given by equation 3.20 or 3.21). These SPSFs are a function both of the source voxel coordinates s and the detector coordinates n. The limited angular extent of the collimator acceptance function F(ÂŁ') and preferential photon scattering at low angles result in peaked scatter projections. In figure 3.6 the x coordinate is defined as the direction perpendicular to the detector surface. Therefore, the interpolated SPSF for the activity voxel s with coordinates (sx,sy,sz) should have a peak at the detector coordinates: ny< = sy and nz> â€” sz. To obtain a better estimate of the interpolated SPSFs the APD-calculated SPSFs must be shifted appropriately. Equation 3.25 can be modified to define the shifted-interpolated estimate of the scatter point spread function SPSFj n t e r p(s, ni) as: 7 SPSF i n f e r p(5,ni) = YI W m S P S F ^ , ni - A m ) , (3.27) m=0 where A m are the two dimensional detector coordinate shifts for each of the APD-calculated SPSF(s^,ni) matrices. The shifts for many of the sampled source positions Sm, are the same due to symmetry (i.e., A 0 = A \ , A 2 = A 3 , A 4 = A 5 and A 6 = A 7 ) . This occurs because two SPSF(s^,ni) matrices will have peaks in the same detector co-ordinate postions if their sy and sz source coordinates are the same. These shifts can be Chapter 3. The Method 65 written as: = Sy ~ syo' Ay = sz - S z 0 = Sy ~ syo> Ay = sz - Szo Ay'2 sy â€” syi, Ay z2 = 5 Z ~~ SZ0 \'3 = Sy â€” syi, Ay Z3 = Sz - S z 0 Ay'4 = Sy ~ 5 yo> Ay = sz ~ s * i â€” Sy â€” s yo> Ay = Â«* - szr Av'e = Sy â€” s y i > Ay ze = ~ Szi Ay'7 sy - syi, Ay = sz - sZl (3.28) Another feature that needs to be considered is the influence of attenuating material between the source voxel and the detector. Using a digital thorax phantom, it was observed that when highly attenuating material was located between the source voxel and detector surface, in the direction parallel to the detector normal, the size of the peak was reduced in the SPSF. The peak region primarily consists of very low angle scatter. However, pixels that are further away from the peak region are generally unchanged. This can be explained by noting that large variations in attenuation values in human tissue are due primarily to the presence of bones. Bones that could typically influence the scatter probability distribution in a SPECT cardiac study are the ribs and spine. These bones are small and localized and only affect a narrow region of the probability projection. In an attempt to approximately correct the interpolated SPSFs for the amount of attenuating material that the scattered photons have passed through, equation 3.27 is modified in the following way: S P S F i n t e r p ( s , fU) = exp (- r" p(x)dx] wm S P S F ( y A â„˘ ) ( 3 . 2 9 ) V J ' J m=o exp ( - n{x)dx) , Chapter 3. The Method 66 where Hvm and np are the coordinates of the final position of a photon that has traveled directly from the source positions sm and s, respectively, in the direction parallel to the detector normal to the detector surface. In the coordinate system defined in figure 3.6, the factor exp(â€” jpP u.(x)dx) is the probability that a photon is attenuated when trav-eling directly from a source voxel at the coordinates (sx,sy,sz) to the detector surface at coordinates (0,sy,sz). This will inevitably be an underestimation of the attenuation effect in the peak region since all photons that have scattered by non-zero angle will have traveled some distance further than has been approximated. This underestimation was a compromise which prevented an over-correction in the regions of the SPSF further from the peak. These regions correspond to larger scattering angles and are not strongly influenced by the attenuation between s and Hp. 3.2.3 T H E I N T E R P O L A T E D T O T A L S C A T T E R D I S T R I B U T I O N F U N C T I O N Equation 3.29 can be used to define the interpolated total scatter distribution function TSDFj n i e r p (A, Hi) for detector bin i (which shall also sometimes be referred to simply as the APDI scatter estimate) as: j TSDF i n t e r p (A, Hi) = ÂŁ A J - S P S F i n t e r p ( 5 ' - , fU), (3.30) where SPSF interp(sj, Hi) is given in equation 3.29 and Xj is the total number of photons emitted from activity voxel j. 3.2.4 I T E R A T I V E C A L C U L A T I O N O F T H E S C A T T E R D I S T R I B U T I O N It should be noted that equation 3.29 does not explicitly depend on the activity intensity Aj. This is a distinct advantage in a reconstruction algorithm where the activity dis-tribution is an unknown quantity that changes iteratively. The interpolation algorithm, Chapter 3. The Method 67 given in equations 3.30, proceeds relatively quickly compared to the APD calculation of SPSFs, described by equation 3.20 or 3.21. As a result it is convenient to store the APD-calculated SPSFs for each of the voxels on the corners of the interpolation grid. As the activity distribution is iterated in the reconstruction algorithm these SPSFs can then be reused. 3.3 IMPLEMENTATION OF SCATTER CORRECTIONS IN S P E C T R E -CONSTRUCTIONS This section describes OSEM-APDI, the implementation of the APDI scatter estima-tion method into an iterative OSEM reconstruction algorithm. In order to evaluate the performance of this scatter corrected reconstruction technique, we have compared a number of different reconstruction protocols. All methods (including OSEM-APDI) use the OSEM algorithm with attenuation correction (AC) and three-dimensional detec-tor response compensation (3D-DRC) [107] but differ in their implementation of scatter compensation techniques. 3.3.1 T H E OSEM-BB M E T H O D As described in chapter 2, a simple and well known holistic scatter correction method, frequently used to partially compensate for scatter, is to rescale the attenuation maps used in the OSEM algorithm from narrow to broad-beam attenuation coefficients. This method, which shall be referred to as OSEM-BB, employs an attenuation correction in which the original narrow-beam attenuation maps has been scaled by a factor of 0.8. This factor was selected to rescale the attenuation coefficient of water from 0.15 cm - 1 to 0.12 cm""1 for 140 keV photons, as in [27]. Chapter 3. The Method 68 3.3.2 T H E OSEM-ISR M E T H O D A second reconstruction method employed OSEM with ideal scatter rejection (OSEM-ISR). In this method images are reconstructed using only the primary photons. Therefore, this method is only possible for simulated data. This is equivalent to an ideal data acquisition in which one can perfectly identify and reject all scattered photons. Although OSEM-ISR is "ideal" in the sense that it represents a scatter-free reconstruction, it is still limited by noise in the data, problems associated with having only a finite number of projection views, and approximations included in the AC and 3D-DRC. Since there was no scattered photons in the acquired data, only 3D-DRC and narrow beam AC are modeled in the reconstruction system matrix. 3.3.3 T H E O S E M - A P D I M E T H O D The OSEM-APDI technique requires three distinct steps for its scatter correction. First, an initial estimate of the activity distribution is made using the OSEM-BB method. Secondly, the original attenuation map and the OSEM-BB image are used to calculate the APDI scatter estimates, as in equation 3.30. Finally, a second OSEM reconstruction (OSEM-APDI) is performed that uses dual matrix ordered subsets (DMOS) expectation maximization (see section 1.7.2). In this reconstruction 3D-DRC, narrow beam AC, and the APDI scatter estimate are used in the forward projection step, but the back projec-tion step only uses 3D-DRC and narrow beam AC. If necessary this whole process can be repeated iteratively, with new APDI scatter estimates calculated using subsequent OSEM-APDI activity distributions. The procedure for OSEM-APDI is shown schemati-cally in figure 3.7. The first OSEM-APDI reconstruction can be written as: Chapter 3. The Method Broad-Beam Experimental u-map Projections OSEM - B B Narrow-Beam Narrow-Beam (i-map Experimental Projections Figure 3.7: Schematic diagram of a reconstruction procedure that includes an APDI scatter correction. Chapter 3. The Method 70 A " + 1 = ÂŁ M â€” > ( 3 - 3 1 ) thn 13 ^ ÂŁ C ^ k + bTSDFinterp(XBB, nt) fc=i where A" is the estimate of the activity for voxel j at the nth iteration of the OSEM algorithm, Ctk is the system matrix element that includes AC and 3D-DRC, Sn represents the ordered subset of projection angles for iteration n, and Yt are the number of counts recorded in detector bin t. The APDI scatter estimate T S D F j â€ž f e r p ( A B B , nt) can be written as: T S D F i n t e r p ( A B B nt) = ÂŁ Af B S P S F i n t e r p ( 5 j , nt), (3.32) where A B B is the OSEM-BB estimate of the activity intensity for voxel j. In equation 3.31 a scaling factor b is applied to the T S D F i n t e r p ( A B B , n*) term in the denominator. The use of broad-beam attenuation maps in image reconstruction does not reliably produce quantitatively accurate images. As a result these reconstructions could be biased. Therefore, the scaling factor b is introduced. In the APDI algorithm the primary photon distribution given by equation 3.10 can be generated at the same time as the scatter distribution function. The scaling factor b is estimated to match the experimentally measured number of counts in each detector bin Yt with YApm(XBB,nt), the sum of the APDI primary and scattered photon distributions. An expression for y A P D i (A B B ,n t ) can be written as: r A P D I ( A B B , nt) = PPD(A B B , nt) + T S D F j n i e r p ( A B B , nt), (3.33) where PPD(A B B ,n t ) is the APD primary photon distribution (given by equation 3.10) using the A B B activity distribution. Chapter 3. The Method 71 The scaling factor b is determined by taking a weighted average of the ratio of Yt over Y A P D I ( \ B B , nt). The experimental measurements of YT, are governed by Poisson statistics. For counts much greater than one, the statistical uncertainty oyt, can be approximated as yfYt. The weighted average scaling factor can then be written as: ÂŁ -5=-^-^-, (3.34) E ^ 2 where: Y h = Y ^ ( \ B B , n t ) ' ( 3 3 5 ) dbt dYt < = , F L . ^ Â» . (3-36) ( F A P D I ( A B B , n ( ) ) and the summation from t = 1 to T" is over all detector bins with recorded counts greater than five. It is assumed in the calculation of o\t that the uncertainty in Y A P D 1 ( X B B , n t ) , is much smaller than the statistical uncertainty in Yt. Perhaps a more obvious choice for this scaling factor would be the total ratio of the sum of the counts (i.e., b = (X) Yt)/(Y1 YAT>m(\BB, nt)), where the summation is over all detector bins). In this implementation a weighted average scaling factor was used instead, because it takes into account the statistical uncertainty in each detector bin measurement. This is achieved by giving a smaller weight to those ratios determined for detector bins with low count values than those with higher count values. The lower portion of figure 3.7, schematically shows the procedure to iteratively up-date the APDI scatter estimate. The APDI-calculated SPSFs are reused to iteratively generate new estimates of the scatter distribution function as the activity distribution Chapter 3. The Method 72 function changes in the reconstruction process. We can define the scatter corrected activ-ity distribution as A A P D I p , where p is an iterative variable which represents the number of times the scatter estimate has been updated. We can define this new estimate of the scatter projection matrices as T S D F i n t e r p ( A A P D I p , nt). Subsequent OSEM-APDI recon-structions that use the updated scatter distribution function can be written as: \ " + ' = v & T - E S ^ . (3-37) <Sl 3 teS" ÂŁ Ctk^k + bpTSDFinterp(\APT)lp, nt) k=i where: T S D F j n ( e r p ( A A P D I p , nt) = ÂŁ A / P i 3 / p S P S F i n ( e r . p ( 5 j , nt). (3.38) j=i 3.3.4 T H E O S E M - A P D T M E T H O D In simulation studies the true activity and attenuation distributions are known. This a-priori knowledge allows for additional investigation of the OSEM-APDI reconstruction algorithm that is not possible with experimental data. In the simulation studies (to be discussed in further detail in chapter 4) both Monte Carlo and noiseless projection data were generated. The noiseless data were simulated by inputting the true activity and attenuation distributions into the full APD calculation. These projections represent a perfectly noiseless data set. The generation of the noiseless projections allowed for modification of the OSEM-APDI algorithm, in which the true APD-calculated scatter projections are used in place of those calculated with the approximate APDI technique. This method will be referred to as OSEM-APDt. The OSEM-ISR and OSEM-APDt reconstructions were compared to investigate the limitations the DMOS approximation, in which scatter is only included in the forward projection step. The OSEM-APDt and Chapter 3. The Method 73 OSEM-APDI reconstructions were compared to investigate inaccuracies and biases that may be introduced by the APDI approximation relative to the full A P D method. C H A P T E R 4 VALIDATION OF A P D I AND O S E M - A P D I USING COMPUTER SIMULATIONS The purpose of the work described in this chapter1 was to validate the APDI scatter es-timation technique and the OSEM-APDI reconstruction using computer simulated data. This goal was achieved in three distinct stages. Section 4.1 describes the anatomically realistic digital phantoms that were used to model the distribution of activity and attenuation in the thorax for a SPECT coronary perfusion study. This section also describes the specific acquisition parameters which were used to simulate both noiseless and Monte Carlo projection data. In the first stage the errors which may be introduced by the APDI interpolation were 1A version of this chapter has been published. E. Vandervoort, A. Celler, R. G. Wells, S. Blinder, K. Dixon and Y. Pang. Implementation of an Analytically Based Scatter Correction in SPECT Recon-structions. In 2003 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, pages M10-256, October, 2003. A version of this chapter has also been submitted for publication. E. Vandervoort, A. Celler, R. G. Wells, S. Blinder, K. Dixon and Y. Pang. Implementation of an Analytically Based Scatter Correction in SPECT Reconstructions. IEEE Trans. Nucl. Sci., 2003. 74 Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 75 evaluated. In section 4.2 scatter projections generated with the full APD method are compared with those calculated with the approximate APDI technique using the same activity and attenuation maps. The comparison was performed separately with activity distributions for the heart, liver and body using a complete set of 64 projections evenly spaced over 360Â°. The second step, described in section 4.3, was to determine the optimal way to im-plement the OSEM-APDI reconstruction technique. This was achieved by comparing two different versions of the OSEM-APDI algorithm used to reconstruct both noiseless and Monte Carlo simulated projection data. The first version requires long computation times, but is meant to provide an extremely accurate scatter correction by performing many iterative updates of the APDI scatter estimate. The second method is meant to minimize the reconstruction time by performing only one APDI scatter estimation step with a thresholded version of the OSEM-BB reconstruction. In simulations both the true activity map and the true distribution of scattered photons in the projection data is known. The different versions of OSEM-APDI were therefore evaluated in two ways. The first consisted of a comparison between the APDI scatter estimates, generated at each step of the OSEM-APDI reconstruction, with the true scatter projections. The second evaluation came from a comparison of the image quality of the OSEM-APDI re-constructions to that of the true activity distribution and that of the other reconstruction techniques: OSEM-ISR, OSEM-APDt, and OSEM-BB. The final validation, discussed in section 4.4, evaluates the consistency of the opti-mized OSEM-APDI algorithm using the noiseless and MC simulated projection data for phantoms with coronary perfusion defects in a variety of positions. In this analysis the OSEM-APDI reconstructions were also compared to the true activity distributions and the OSEM-ISR, OSEM-APDt, and OSEM-BB reconstructions. The goal of this work Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 76 was to verify if the OSEM-APDI algorithm behaves consistently regardless of the defect position. 4.1 DIGITAL PHANTOMS AND SIMULATION ACQUISITION PARAME-TERS Two sets of simulated projection data were considered in this chapter. The first data set was obtained by using digital activity and attenuation maps with the Monte Carlo (MC) simulation software SimSET [99], which models photon interactions in the body and the acquisition process to generate realistic SPECT projections with noise characteristics similar to those of a clinical scan. For the second set, the same activity and attenuation maps were input in the full APD method to generate noiseless projection data. The long calculation times required to perform MC simulations and to perform the full APD calculation, made it advantageous to reuse projection data in the various sections of this chapter. The acquisition parameters for APD and SimSET were made as similar as possible so that the reconstructions methods which used either dataset could be compared to each other more easily. In order to do this, the same activity maps, attenuation map and acquisitions parameters must be used to generate both sets of data. In all the studies discussed in the following chapter, the dynamic mathematical cardiac-torso phantom (MCAT) [105] was used to generate anatomically realistic ac-tivity and attenuation maps for the thorax of a standard sized female. Separate activity maps were generated representing different tissue types. Since the detection of emitted photons from radioactive decay is a linear process, the simulated projection data for each tissue could be generated separately and added together to create the final datasets. Activity maps were created for the healthy myocardium and myocardium with defects in the anterior, apical, inferior, lateral and septal walls of the myocardium, as well as the Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 77 liver and body (the lungs and all other remaining tissue). The MCAT program models the beating heart and was used to produce an activity distribution averaged over several time-frames of a complete heart cycle. The anterior, inferior, lateral and septal defects were elliptical and positioned at the midpoint along the length of the left ventricle. Each of these defects occupied 14% of the total volume of the left ventricle. The apical defect was circular and located in the center of apex with a volume equal to 11% of the total left ventricle. The activity maps corresponded to 9 9 m T c activity concentrations of 4.0, 0.0, 4.0, 0.5 and 1.0 ^Ci/cm 3 in the heart, defect, liver, lungs and all other remaining tissue respectively. These values were selected to give clinically realistic concentrations and total injected activities for each organ. The total injection was 1.0, 7.0, 2.0 and 30.0 mCi for the heart, liver, lungs and all other remaining tissue respectively. Although the MCAT program could be used to generate a pixelized map of the linear attenuation values in cm - 1 , this map was not used because it was incompatible with the SimSET simulation software. SimSET uses internal tables that specify p(E), the variation of linear attenuation coefficients with photon energy for specific materials and tissue types. The attenuation file required for SimSET does not consist of linear atten-uation coefficients, but rather is a map of integer indices. Each integer index specifies a particular material or tissue type which is then used to obtain the appropriate tab-ulated p(E) values2. The MCAT attenuation map could not be used directly because the MCAT program models organs as continuous geometric objects. This leads to the partial volume effect, in which a given voxel may contain a number of different tissue types each with different attenuation values. The motion of the heart leads to a further 2 SimSET also uses an activity file input that consists of a a map of integer indices. These indices are used to translate each integer value to a unique /uCi/cm3 activity value. The tables are user-defined and limited to a maximum length of 256 unique /iCi/cm 3 values. The MCAT program generated activity maps with fewer than 256 unique /iCi/cm 3 values for all cases considered. As a result, there was no loss of information when importing MCAT activity files into SimSET. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 78 blurring of the attenuation values in the region near the heart. These effects made it difficult to definitively classify different voxels as belonging exclusively to one or another tissue type. To allow for the same attenuation map to be used by both SimSET and APD, the MCAT program was used to generate maps that represented tissue percentage values. A voxel that consisted of more than 50% of a particular tissue type was assigned the corresponding SimSET index for that tissue type. In the unlikely case where a particular voxel was on the boundary of three tissue types, with tissue percentages each less than 50%, this voxel was assigned the SimSET index for muscle. Each voxel was classified as either air, fat, blood, bone, heart, lung, muscle, gastro-intestinal tract, or liver tissue. The APD and APDI scatter projection calculation requires the distribution of linear attenuation values for the unscattered photon energy (140 keV). For consistency, the attenuation map used in the APD and APDI methods were obtained by translating this map of integer indices into linear attenuation values using SimSET's internal tables for 140 keV photons. Both APD and SimSET used a geometric collimator response function for a circular parallel-hole, Low Energy Ultra High Resolution (LEUHR), collimator with no septal penetration. The energy response of the detector was modeled as a Gaussian function with FWHM equal to 10% of the incident photon energy, a 20% energy window was used in the detection process, and 100% efficiency was assumed for all photon energies. The detector crystal thickness was modeled as zero in the image formation process. The detector was rotated in a circular orbit with a radius of rotation of 23.8 cm. A total of 64 views evenly spaced over 360Â° were acquired. The projection matrix was 64x64 with 0.72 cm bins. For the MC simulations, forced non-absorption variance reduction (see chapter 2) was utilized to speed up computation times. An adequate number of photon Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 79 histories, as recommended by the SimSET documentation, were simulated to obtain the correct noise characteristics for an experimental scan with a 20 second acquisition for each view ( Â« 2xl0 5 counts per projection). The APD projection data represents perfectly noiseless data corresponding to the specific activity concentrations and length of scan described above. This data could represent for example, the average result of an infinite number of scans under identical conditions. For all simulations in this chapter, third and higher order Compton and Rayleigh scatter were not considered. This was because SimSET does not model Rayleigh scatter for SPECT and because higher order Compton scatter does not contribute significantly to the 9 9 m T c photopeak data. As a result, the simplified version of the scatter point spread function SPSF(s}, rTj), given in equation 3.21, was used throughout this chapter. It should also be noted that SimSET does not model either Rayleigh or Compton scattering form factors. Therefore, equations 1.1 and 1.3 did not include the terms T(x,Z) or S(x,Z) in the APD look-up table calculations performed for this chapter. The method used to generate APD projections did differ from the SimSET calculation in a few significant ways. The steps involved in the generation of APD (or APDI) primary and scatter projections are shown schematically in figure 4.8. The inputs to the APD calculation are an activity map and an attenuation map. A Gaussian based diffusion rotation algorithm [108] is used to simulate the acquisition of projections from different angles. This is necessary because the location of the detector is fixed in the APD calculation. The rotated attention map is then used to estimate the distribution of electron density values. The rotated attention map is also used to create an "evaluated attenuation map" that consists of attenuation lengths (as given by equation 3.13). The evaluated attenuation map, density map, activity map and the original attenuation map are all used to generate the APD primary and scattered photon projections. In SimSET Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 80 the activity map is not rotated. The detector is modeled as though it is simultaneously at every projection angle position. The electron density values for each tissue type are obtained from SimSET's internal tables for specific tissue types. Similarly the variation of linear attenuation with energy is approximately determined in APD (using equation 3.14) while SimSET uses its internal tables for the p,(E) values for specific tissue types. Inputs: Activity Map Attenuation Map ^u(x)d). Evaluated AttenuatioB-Map APD or APDI Pi oicilion C.iUuKilion Output: Primary and Scatter Projections Figure 4.8: Schematic diagram to illustrate the steps involved in the APD (or APDI) projection calculation for the data used in this chapter. The figure shows all the intermediate steps which are required to accurately model the simulated data acquisition. In this chapter, all figures that are shown as a function of projection angle are defined according to the following convention. The start angle of all acquisitions was defined as 90Â° with the camera at the left hand side of the phantom. The rotation proceeded underneath the phantom to the posterior position at 180Â°, through to the right hand side of the patient (270Â°), to the anterior position (360Â°) and ended just short of its starting position (444.375Â°). This convention corresponds to the angles which were used in the rotation program. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 81 4.2 C O M P A R I S O N O F T H E F U L L APD M E T H O D T O APDI In order to quantify the errors introduced by the interpolation approximation, a com-parison was made between the scatter projections generated using the same activity and attenuation distribution for both APD and APDI. The comparison was performed sepa-rately with activity distribution maps for the healthy myocardium (with no defect), the liver and the body (which including the lungs and remaining tissue) for the complete set of 64 projections. 4.2.1 R E S U L T S To obtain an estimate of the errors introduced by APDI, two figures of merit were con-sidered. The first was x 2 defined as: v 2 = v (TSDF / U â€ž(A, fU) - TSDF i n t e r p (A, n*))2 X 2-*/ _2 ' (4.cjyj - \ \2 where the summation over i is over all detector bins in a given projection, A is the true MCAT activity distribution and <7j is the expected standard deviation of detector bin i. Since photon detection is a statistical process that follows a Poisson probability distribution, o-j is equal to ^TSDFfuu(X, fij). Hence x2 divided by the number of degrees of freedom (d.f.) provides a measure of the average size of the differences between APDI and APD relative to the statistical uncertainty expected in that projection data for a given activity distribution and length of scan. In this case, the number of degrees of freedom was taken to be the total number of detector bins with TSDF/ UÂ»(A, n*) greater than one. Since the APD calculation involves integration of scatter probability functions, there is the possibility of obtaining values between one and zero. All these values were set to zero probability of photon detection and were neglected in the analysis. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 82 If the difference between TSDFj n t e r p (A, Hi) and TSDF/ UÂ»(A, Hi) is on average compa-2 rable or larger than the statistical variation in the projection data, the -fj value will be 2 one or greater. If -Xy is considerably less than one, it is unlikely that these differences will be significant in a reconstruction based scatter correction. Figure 4.9 (a) shows fj as a function of projection angle for each of the activity distributions considered in this section. Although useful for quantifying the size of the errors relative to statistical uncertainty X2 gives no indication of bias that may be introduced by the interpolation approximation. Therefore, the second figure of merit considered was bias defined as: bias Sj=i TSDFj n t er P(A, Hi) E f = i T S D F / l l â€ž ( A , n i ) (4.40) The bias as a function of projection angle was determined for each set of projections and is summarized in figure 4.9 (b). A comparison of mean and bias values, averaged over all 64 projections, for each activity distribution is summarized in table 4.1. 0.6 \ to .5 u "8 ÂŁ 0,4 (a) Heart - - Liver - Body 100 150 200 250 300 350 Projactlon Angle (degrees) 200 250 300 Projection Angle (degrees) (b) Figure 4.9: Comparison of the (a) and bias value (b) as a function of projection angle, using the activity distributions for the heart, liver and body. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 83 Table 4.1: Comparison of mean fj and the mean bias, averaged over all 64 projections and their standard deviations for the heart, liver and body activity distributions. Standard Deviation Bias Standard Deviation Heart Liver Body 0.03 0.04 0.05 0.07 0.17 0.15 1.02 0.06 1.02 0.05 1.03 0.03 Examples of projections calculated using the APD and APDI methods are shown in figures 4.11 and 4.10 for detector angles 180Â° (posterior view) and 360Â° (anterior view) respectively. Displayed are the sums of the distributions obtained for the MCAT heart, liver and body activity distributions. The lower parts of these figures show horizontal profiles drawn through the pixel with maximum count value on linear and logarithmic scales. These particular views were selected to represent "best" and "worst-case scenar-ios". In the 180Â° view the spine is positioned between the source distribution and the detector. In the 360Â° view it is primarily soft tissue between the source and detector. Therefore, errors due to the approximate correction made for attenuation in the APDI method are expected to be more significant for the 180Â° view than for the 360Â° view. The amounts of CPU time (1.7 GHz processor, 1024 Mb of RDRAM) that were necessary to calculate a single scatter projection for each of the activity maps considered in this section, using the APDI and APD method, are summarized in table 4.2. A range of values are given because the computation time for both methods depends on the amount of attenuating material, between the activity distribution and the detector, which changes when the data is acquired from different angular views. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 84 0 20 40 60 80 0 20 40 60 Detector Bin Detector Bin Figure 4.10: Comparison of full APD and approximate APDI scatter distribution functions for the MCAT phantom (sum of heart, liver and body activity distributions) with the detector at the anterior view (360Â°). Horizontal profiles through the maximum count value are shown using linear and logarithmic scales in (a) and (b) respectively. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 8 5 Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 86 Table 4.2: Timing results for one scatter projection calculation, for the APD and APDI methods, using the heart, liver and body activity distributions. The calculations were performed on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM. Activity Map CPU Time for APD Per Projection CPU Time for APDI Per Projection Number of Non-Zero Activity Voxels Heart 12-31 min. 1-4 min. 2960 Liver 12-33 min. 2-5 min. 3316 Body 4.1-4.7 hours 16-21 min. 37579 4.2.2 D I S C U S S I O N Qualitatively the APDI projection data shown in figures 4.10 and 4.11 are smoother than the full APD data. This is due to the averaging effects associated with tri-linear interpolation. These projections also show that there is a slight overestimation of APDI data near the edges of the scatter projections, which is due to detector coordinate shifts applied in equation 3.29. Due to the approximate correction for attenuation of scattered photons in APDI, the largest differences between the APD and APDI data occurs for detector coordinates in which highly attenuating bone material (e.g., spine, ribs, and sternum) was situated between the source and the detector. As discussed previously, this can produce APDI probability projections with overestimated peak values. In figure 4.11 these differences were as high as 20% of the maximum count in a given projection bin. In figure 4.9 there were a small number of projection angles for which 2 2 the -fy approaches the value of 1.0 and bias values were as high as 1.23. A jfj Â« 1.0 implies that the errors introduced by the APDI method were on average comparable to the expected statistical uncertainty in these projection values. A bias value of 1.23 means that there is a 23% overestimation of the total number counts in the APDI data relative to APD. As significant as these problems seem to be, there are two reasons why they may not affect the reconstruction of a typical cardiac scan. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 87 The first reason is that for all activity distributions the highest values of x2 occurred for angles in which the detector was positioned underneath the phantom. In these views, the spine was located between the activity distribution and the detector. Therefore, the approximate correction for attenuation in APDI was more significant. However, in most cardiac acquisition protocols the camera rotates over 180Â° between the left posterior oblique angle and the right anterior oblique angle. This corresponds to the subset of projection angles in figure 4.9 from 90 to 135Â° and from 315 to 444.375Â°. These angles had much smaller x2 a n d l e s s variation in bias values than the remaining angles. Secondly, the attenuation map used here was selected to be compatible with the MC simulation software SimSET, which used an unusually high value for the linear attenuation coefficient of bone (0.34 cm - 1 at 140 keV). The overestimation in scatter projection data, associated with the approximate attenuation correction made in the APDI method, is highly sensitive to large variations in the attenuation maps. To evaluate the significance of this effect, a comparison was made between the APD and APDI data using another attenuation map with a slightly lower estimate of the linear attenuation coefficients for bone (spinal bone: 0.17 cm - 1 and skeletal bone: 0.21 c m - 1 at 140 keV). This map was the direct output of the MCAT program (not modified to be compatible with SimSET) and uses attenuation coefficients from the tabulated Medical Internal Radiation Dose (MIRD) radiological data [112]. Unfortunately, the long computation times required for the APD method prevented an analysis from being performed on a full set of 64 projection angles. However, the 2 and bias values were determined for the three projection angles that demonstrated the highest x2 values in figure 4.9. The results are summarized in table 4.3. The values were reduced by a factor of ten and a substantial improvement in bias values ( Â« 10% to 15%) was observed using this new attenuation map. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 88 Table 4.3: Comparison of APD and APDI projections using the SimSET p values and p values from the MCAT program. SimSET /i-map d.f. MCAT p-map d.f. Activity-Distribution Projection Angle Bias Bias Heart Liver Body 219.375Â° 135.0Â° 258.750Â° 0.33 0.86 0.81 1.22 1.23 1.13 0.03 0.08 0.08 1.04 1.05 1.04 The significance of the errors associated with the interpolated APD method might be better evaluated using the mean values summarized in table 4.1, since image reconstruc-tion uses data for all projection angles. 'Although the projection to projection variation in 2 2 fj values was high, the mean value was much lower than one. The average jfj- overlaps with zero plus or minus its standard deviation for all data sets considered except the body. Similarly in figure 4.9 (b) a large degree of variation in bias values was observed from projection to projection. These values ranged from over 20% overestimation to about a 10% underestimation of the total number of counts. However, the average ratio over all projection angles plus or minus its standard deviation overlapped with one for all cases considered. An important practical question that was addressed is under what conditions would the differences between APD and APDI become significant. If one increases the data acquisition time or increases the activity concentrations, the resulting projection data would have lower statistical noise. If one applies a linear scaling factor to the APDI and APDI scatter projection data, used in this section, one can determine the level of statistics where the differences between APD and APDI would become significant. The X2 for a scan with the same relative ratios of activity as the MCAT phantom but c times longer acquisition time or c times greater absolute activity can be written as: Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 89 2 ^ ( c x TSDF / t l Â»(A,nQ - e x TSDF i n < e r p (A, nQ)2 Xc h (cx^)2 c2 (TSDF / U â€ž(A, nQ - TSDF i n f e r p (A, n,))2 ^ c x T S D F / u â€ž ( A J n i ) = e x f i , (4.41) i=l Using the average values over all projection angles (table 4.1), the difference between the APD and APDI data would be comparable to the statistical uncertainty in the measured data for a scaling factor of c = 33 for the heart, c = 20 for the liver, and c = 5.9 for the body. 4.2.3 C O N C L U S I O N S The analysis of calculation times demonstrate that depending on the size and complex-ity of the activity distribution, the APDI approximation with the 4x4x4 interpolation grid spacing reduces computation time by a factor of 10 to 20 when compared with the full APD method. The differences between the APD and APDI scatter estimates are small relative to the statistical uncertainty in SPECT projection data if one considers the x2 averaged over all projection angles. Similarly no significant bias is introduced when averaged over all projection angles. It is unlikely that large differences observed in individual pixels in a small subset of projection angles will influence the OSEM-APDI reconstruction, since the entire set of projection data is used to reconstruct the image. 2 The average *j values in table 4.1 indicate that the differences between APD and APDI would only become significant (i.e., they would be comparable to the statistical uncer-tainty in myocardial projection data) for acquisition times between 33, 20, and 6 times longer than a typical clinical scan for the heart, liver, and body activity distributions, respectively. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 9 0 4.3 OPTIMIZING THE O S E M - A P D I RECONSTRUCTION ALGORITHM The following section investigates the best way to implement the OSEM-APDI image reconstruction algorithm. This goal was achieved by comparing two different versions of the OSEM-APDI method, using noiseless and MC simulated projection data for the MCAT coronary phantom with a defect in the inferior wall of the myocardium. The first version: OSEM-APDI(0%)P, iteratively updates the APDI scatter estimate several times (represented by the variable p), during the reconstruction process. It is anticipated that this method will provide an extremely accurate scatter correction, but will require long computation times. The second version: OSEM-APDI(15%), only performs one APDI scatter estimation using a thresholded activity distribution to reduce computation times. There were three main goals for the analysis performed in this section. The first was to determine if it is reasonable to use thresholded activity maps in terms of the accuracy of the APDI scatter projections and the image quality of the resulting OSEM-APDI reconstruction. The second goal was to determine if it is necessary to iteratively update the APDI scatter estimates and if so, how many scatter iterations are optimal. The final goal was to determine how many iterations should be used within the OSEM reconstruction algorithm itself (not be confused with the scatter iterations). To avoid any possible confusion, the following discussion will always distinguish between OSEM iterations and scatter iterations. 4.3.1 D A T A A N A L Y S I S Since the true distribution of scattered photons in the simulated projection data was known, the OSEM-APDI reconstructions were evaluated was by comparing the APDI scatter estimates, generated at each step of the OSEM-APDI reconstruction, with the true component of scatter in the data (which could be either noiseless or MC simulated Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 91 data). This evaluation was carried out using two figures of merit: and bias, previously defined in equations 4.39 and 4.40. Since the OSEM-APDI method uses the weighted average scaling factor b to match the experimentally data with the APDI projections, the degrees of freedom (d.f.) are reduced by one in this case. For the SimSET data, one would expect that both figures of merit should be approximately equal to one if the simulated data has the correct Poisson distributed statistical error and if the APDI scatter projections accurately reflect the underlying distribution of scattered photons. In the case of the noiseless projection data, these values should be compared to the results of the previous section where APD and APDI were compared using the same activity distribution for the true MCAT phantom. In this section we can now evaluate the accuracy of the APDI scatter estimates when a reconstructed image is used instead of the true activity distribution. The OSEM-APDI reconstructions were also evaluated by computing two different image quality measures. This analysis was performed in a similar manner to [93]. The first measure was the defect contrast, defined as: p^jj, where a and b are the mean count values in a volumes of interest (VOI) for the healthy myocardium and defect tissue respectively. For this analysis, all VOIs were groups of voxels obtained from the true MCAT phantom. It should be noted that although the defect contrast calculated from the known activity concentrations is 1.00 partial volume effect and blurring due to heart motion reduced this contrast to 0.71 in the true MCAT phantom. The second measure was the normalized standard deviation (NSD) defined as the standard deviation divided by the mean count value in a VOI (in this case all liver voxels) that contain a constant activity concentration in the true MCAT phantom. The image quality measures for the two versions of OSEM-APDI were compared to the true MCAT activity map and the reconstructions: OSEM-ISR, OSEM-APDt, and OSEM-BB. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 92 Another important consideration in image reconstruction is the time required to per-form the reconstruction. As a result all reconstructions were also evaluated with respect to the reconstruction CPU time. 4.3.2 A N A L Y S I S O F T H E O S E M - A P D I ( 0 % ) P M E T H O D OSEM-APDI(0%) O S E M - B B Activity Map Narrow-Beam Attenuation Map Experimental Projections Experimental Projections APDI S i n n e r I'miectmn Calculation n o APDI -Scatter Projections O S I I M - A P D K U ' . i j Reconstruction Ai'ni spinel Pn>|uuion I'.ili'ulaliiin A P D I -Scatter Projections OSI1M-APDI.01.;) . 2 Reconstruction Analysis: â€” and Bias d.f. Analysis: â€˘ Contrast and N S D Analysis: â€” and Bias d.f. Analysis: - Contrast and N S D First Scatter Iteration .Second Scatter Iteration Figure 4.12: Schematic diagram to illustrate the steps involved in OSEM-APDI(0%)P and its analysis. Figure 4.12 is a schematic diagram to help illustrate the analysis of the OSEM-APDI(0%)p algorithm. This flow chart shows the input data, the steps involved in the reconstruction, the data output at each step, and how that data was analyzed. In this Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 93 approach the eighth OSEM iteration of the OSEM-BB reconstruction is used without any threshold to generate the first APDI scatter estimate: TSDF i n t e r P ( A B B , n t ) . This APDI data is then used to perform the OSEM-APDI(0%)i reconstruction. Then a second scat-ter iteration begins. The OSEM-APDI(0%)i (again the eighth OSEM iteration with no threshold) is used to calculate a new APDI scatter estimate: T S D F i n t e r p ( A A P D l 1 , nt). Then a second OSEM-APDI reconstruction: OSEM-APDI(0%)2 is performed. This pro-cess is repeated 16 times to evaluate how many scatter iterations are necessary. The variable p refers to the number of times the scatter estimate has been updated. Each APDI scatter estimate was analyzed using the and bias figures of merit and each OSEM-APDI(0%)P reconstruction was analyzed using the contrast and NSD image qual-ity measures. 4.3.3 A N A L Y S I S O F T H E OSEM-APDI(15%) M E T H O D O S E M - A P D I ( 1 5 % ) Thresholded OSEM-BB Activity Map Narrow-Beam Attenuation Map Experimental Projections :tlon plictiiSBi OSLM-M'DMV. I S ! Analysis: APDI -Scatter Projections Analysis: - Contrast and NSD Reconstruction Figure 4.13: Schematic diagram to illustrate the steps involved in the analysis of the OSEM-APDI(15%) reconstruction technique. Figure 4.13 schematically shows the steps involved in OSEM-APDI(15%) and its anal-ysis. In the first and only scatter iteration, the activity map used in the APDI calculation Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 94 is the eighth iteration of the OSEM-BB reconstruction. This activity map is thresholded to only include those voxels with activity greater than 15% of the maximum. This par-ticular threshold was selected because it includes all of the voxels in the heart and liver, very few of the voxels for the lungs, but the majority of the voxels for the remaining tissue ( Â« 64%). In total the threshold included 56% of all non-zero voxels in this particular OSEM-BB reconstruction. This threshold can reduce APDI computation time signifi-cantly. However, it can also lead to biased APDI scatter projections. The APDI scatter 2 estimate was analyzed using the fj and bias figures of merit and the OSEM-APDI(15%) reconstruction was evaluated using the contrast and NSD image quality measures. 4.3.4 A N A L Y S I S O F T H E O S E M - I S R , O S E M - A P D T , A N D O S E M - B B R E C O N -S T R U C T I O N S Figure 4.14 schematically shows the steps involved in the OSEM-ISR, OSEM-APDt and OSEM-BB reconstructions and their analysis. These images were evaluated using the contrast and NSD image quality measures. Both the OSEM-ISR and OSEM-APDt re-constructions represent a type of ideal scatter correction in which the true component of scatter is known. The OSEM-APDt reconstruction was not performed for the SimSET data because the true component of scatter in this case contains statistical noise3. 4.3.5 R E S U L T S 2 Figure 4.15 shows the fj and bias values as a function of projection angle for the noiseless and SimSET projection data. In figures 4.15 (a) and (b) the x and o markers are for the 3The OSEM-APDt using MC data could be thought of as an ideal scatter estimation method in which the true scatter component in the projection data is known but that this component is corrupted by Poisson noise [93]. However, this type of scatter correction is not particularly relevant here, since the topic of this thesis is an analytical scatter correction which by definition does not contain any statistical noise. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 95 Narrow-Beam Attenuation Map 0Â«O Primary Photon Projections (No Scatter) OSEM-ISR OSLM-ISR <2IS> Analysis: - Contrast and NSD Reconstruction OSEM-APDt The True MCAT Activity Map Narrow-Beam Attenuation Map Experimental Projections l i l t luh APD Pro i l m ( j I U L I i n O s l M - M ' l l i The True Scatter Projections Analysis: - Contrast and NSD Reconstruction OSEM-BB Broad-Beam Attenuation Map Experimental Projections O N I M - B U Analysis: - Contrast and NSD Reconstruction Figure 4.14: Schematic diagram to illustrate the steps involved in the analysis of the OSEM-ISR, OSEM-APDt, and OSEM-BB reconstruction techniques. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 96 first scatter iteration of the 0SEM-APDI(15%) and OSEM-APDI(0%) reconstructions, respectively. These data represent the APDI scatter distribution generated using OSEM-BB as the activity distribution (i.e., TSDFj n t e r p (A B B , nk)). The cross and square markers are for the second and third scatter iterations of the OSEM-APDI(0%)P reconstruction (i.e., TSDF i n ( e r p (A A P D I l , n f c ) and TSDF i n < e r p (A A P D l 2 ,n k ) , respectively). There were min-imal differences between the APDI data beyond the third scatter iteration, therefore all these data are plotted using dashed lines. Figure 4.16 shows the contrast and NSD as a function of OSEM iteration for all the reconstructions considered using the noiseless and SimSET projection data. In the case of the noiseless data, these reconstructions are: OSEM-ISR, OSEM-APDt, OSEM-APDI(15%), OSEM-APDI(0%)P and OSEM-BB. The. same reconstructions are shown for the SimSET data with the exception of OSEM-APDt, which was not performed in this case. A total of 16 OSEM iterations were performed for each reconstruction method. The data at each scatter iteration of the OSEM-APDI(15%) and OSEM-APDI(0%)P reconstructions are shown using the same data markers, where applicable, as figure 4.15: -^markers for OSEM-APDI(15%), o-markers for OSEM-APDI(0%)i, cross-markers for OSEM-APDI(0%)2, square-markers for OSEM-APDI(0%)3, and dashed lines for higher scatter iterations of OSEM-APDI(0%)P. The timing results for the OSEM-APDI(0%)P and OSEM-APDI(15%) reconstruction methods are summarized in table 4.4. The CPU time for OSEM-APDI(0%)P is split up into four parts. The first part is the time necessary to perform the initial OSEM-BB reconstruction, the second part is the time needed to compute all 64 APDI scatter projections using the OSEM-BB activity map, the third part was the time necessary to perform the OSEM-APDI reconstruction, and the final part was the time necessary to Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 97 (a) (b) E 6 1.4 - â€” 0 S E M - A P D I ( 1 5 % ) -e- O S E M - A P D I ( 0 % ) 1st Iteration - t - O S E M - A P D I ( 0 % ) 2nd Iteration - B - O S E M - A P D I ( 0 % ) 3rd Iteration - - O S E M - A P D I ( 0 % ) Higher Iterations 200 250 300 Projection Angle (degrees) 350 400 - Â» - O SEM-APDI (15% ) -o- O S E M - A P D I ( 0 % ) 1st Iteration - t - O S E M - A P D I ( 0 % ) 2nd Iteration - a - O S E M - A P D I ( 0 % ) 3rd Iteration O S E M - A P D I [ 0 % ) Higher Iterations 200 250 300 Projection Angle (degrees) (c) (d) Figure 4.15: Comparison of the %j and bias values as a function of projection angle for each scatter iteration of the OSEM-APDI reconstructions. The APDI data for both the OSEM-APDI(15%) and iterative OSEM-APDI(0%) P reconstructions are shown. Fig-ures (a) and (b) are for the noiseless projection data, while (c) and (d) are for the SimSET projection data. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 98 - * - OSEM-ISR - Â» - OSEM-APDt OSEM-BB OSEM-APDI(1S%) 1st Iteration -e- OSEM-APDI(0%) 1st Iteration - i - OSEM-APDI(0%) 2nd Iteration -a- OSEM-APDI(0%) 3rd Iteration - - OSEM-APDAPDI(0%) Higher Iterations OSEM-ISR - 9 - OSEM-APDt OSEM-BB - Â« - OSEM-APDI(15%) 1st Iteration -e- OSEM-APDI(0%) 1st Iteration OSEM-APDI(0%) 2nd Iteration â€”Bâ€” OSEM-APDI(0%) 3rd Iteration - - OSEM-APDAPDl(0%) Higher Iterations 10 OSEM Iteration 12 14 OSEM Iteration (a) - * - OSEM-ISR - A - OSEM-BB â€” OSEM-APDI(15%) 1st Iteration -e- OSEM-APDI(0%) 1st Iteration - i - OSEM-APDI(0%) 2nd Iteration -e- OSEM-APDI(0%) 3rd Iteration - - OSEM-APDAPDI(0%) Higher Iterations OSEM Iteration 10 12 14 (b) OSEM Iteration (c) (d) Figure 4.16: Comparison of the contrast and NSD image quality measurements as a function of OSEM iteration, for the OSEM-ISR, OSEM-APDt (where applicable), OSEM-APDI(15%), OSEM-APDI(0%) P and OSEM-BB reconstructions. The contrast between the healthy my-ocardium and the defect are shown for the noiseless projection data in figure (a) and for the SimSET projection data in figure (b). Figure (b) and (d) show the liver NSD values for the noiseless and SimSET data, respectively. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 99 perform each additional APDI scatter calculation using subsequent iterations of OSEM-APDI(0%)p. The full APD-calculated SPSFs are only determined in the first APDI scatter estimation step and then reused in later scatter iterations. Therefore, the first APDI calculation requires the most time. Since the OSEM-APDI(15%) reconstruction only performs one APDI scatter estimation, only the first three entries in the table are provided. Table 4.4: Timing results for each step of the OSEM-APDI(15%) and OSEM-APDI(0%)P reconstructions, for a 64x64x24 reconstruction matrix running on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM. CPU Time OSEM-APDI(0%)P OSEM-APDI(15%) OSEM-BB 280 s 280 s APDI First Scatter Calculation 26.7 hr 9.2 hr OSEM-APDI 280 s 280 s Each Subsequent APDI Calculation 2.6 hr N/A Voxels Included in APDI Threshold 45856 25551 Figure 4.17 shows reoriented short axis slices through the true MCAT activity distri-bution and the OSEM-ISR, OSEM-BB, OSEM-APDt, OSEM-APDI(15%) and OSEM-APDI(0%)3 reconstructions for the noiseless data. The third scatter iteration of the OSEM-APDI(0%)P was selected because there were minimal differences, in terms of both the reconstructed image itself and the APDI scatter estimates, between this scatter itera-tion and subsequent ones. The data shown here is for the eighth OSEM iteration of each reconstruction. The true MCAT activity distribution was normalized to have the same total number of counts as the OSEM-ISR reconstruction. All images are shown using cubic-spline interpolation to reduce pixelization. Each image is displayed using a gray Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 100 scale normalized to the maximum count in that reconstruction. Figure 4.18 shows ver-tical profiles through the original data (i.e.; the non-interpolated reconstructed data) in the inferior wall defect region. The position of the profiles are shown on figure 4.17 as the dashed line on the reconstructed slices. Similarly, figures 4.19 and 4.20 show reoriented short axis slices and profiles respectively, through the true MCAT activity map, and the OSEM-ISR, OSEM-BB, OSEM-APDI(15%) and OSEM-APDI(0%)3 reconstructions for the SimSET projection data. 4.3.6 D I S C U S S I O N 2 The and bias values for the noiseless projection data, shown in figure 4.15 (a), closely resembles the data from the previous section (figure 4.9) for the true MCAT body activity distribution. The body makes up the majority of the activity in the MCAT torso. In the previous section the true activity maps were used in both the APD and APDI method. In this section data reconstructed from noiseless projection data were used to generate the APDI scatter estimates. Therefore, the errors due to using these approximate ac-tivity maps (i.e., the reconstructed data) are small relative to the errors associated with the APDI interpolation procedure. In particular, the jfj values were well below one at almost every projection angle for both the OSEM-APDI(15%) and OSEM-APDI(0%)P reconstructions. This indicates that the difference between the APDI scatter estimate and true component of scatter in the noiseless data is much smaller on average than the expected statistical error in the data. The largest were for the OSEM-APDI(15%) reconstruction and the first scatter iteration of OSEM-APDI(0%)P. This was not unex-pected, since these reconstructions use an activity distribution (OSEM-BB) which could contain scatter-related artifacts to generate their APDI scatter estimates. The for the iterative OSEM-APDI(0%)P reconstruction converges by the third scatter iteration. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 101 Figure 4.17: Comparison of short axis slices for the true MCAT activity distribution with each of the reconstructions considered for the noiseless simulated data: OSEM-ISR, OSEM-APDt, OSEM-APDI(0%) 3, OSEM-APDI(15%) and OSEM-BB. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 102 Y-Axis Cross Section 60 r 1 1 1â€” Figure 4.18: Comparison of vertical profiles through inferior wall defect region of the true MCAT activity map and each of the reconstructions considered for the noiseless simulated data. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 103 Figure 4.19: Comparison of short axis slices for the true MCAT activity distribution with each of the reconstructions considered for the SimSET simulated data: OSEM-ISR, OSEM-APDI(0%) 3 ) OSEM-APDI(15%) and OSEM-BB. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 104 Y-Axis Cross Section 501 1 1 1 â€” Figure 4.20: Comparison of vertical profiles through inferior wall defect region of the true MCAT activity map and each of the reconstructions considered for the SimSET simulated data. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 105 The maximum percent difference in values for any projection angle, when compared to the next scatter iteration, was less than 1% for OSEM-APDI(0%)3. The same trends were observed in the plots of bias as a function of projection angle. The accuracy of the APDI scatter estimates are improved using the iterative OSEM-APDI(0%)P relative to the more approximate OSEM-APDI(15%). Similar trends, in terms of convergence, were observed in figures 4.15 (c) and (d) when the SimSET projection data was used. However, for this data all values were greater than one and the bias values were all substantially lower than one. If the APDI data accurately represents the distribution of scattered photons and the MC data has the correct statistical variance, it was expected that these values should both be close to one. However, there were systematic differences between the APD (and hence APDI) and SimSET scattered photon distributions. In appendix A these differences are discussed 2 in greater detail. However, the fj values for the MC data were still fairly low on average about 1.8 for OSEM-APDI(15%). This indicates that the difference between the APDI and SimSET scatter projections were on average about 80% larger than the statistical error in the MC data. For OSEM-APDI(0%)P, the average ÂŁj was 1.6 to 1.4 (i.e., 60% to 40% larger than the statistical noise). The bias values demonstrate a consistent underestimation of the scattered photon component for the APDI data relative to SimSET. The bias values for both OSEM-APDI(15%) and OSEM-APDI(0%)P were on average 0.95 to 0.9, corresponding to a 5% to 10% underestimation of the total number of counts in the APDI data. The plots of contrast as a function of OSEM iteration, shown in figure 4.16 (a) and (c), demonstrate that the ideal OSEM-ISR reconstruction resulted in higher contrast values than any other method considered for both noiseless and SimSET data at almost every Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 106 iteration. This result was anticipated, since a simulation study [93] found that ideal scat-ter rejection provides better contrast recovery than OSEM reconstructions that include scatter in the system matrix. The one exception was the OSEM-APDI(15%) reconstruc-tion. At very high OSEM iteration, its contrast exceeded the OSEM-ISR result. It is hypothesized that this occurred because of the less-accurate, thresholded activity map produces biased APDI scatter estimates. Since these biased APDI scatter projections are rescaled to match the data, certain regions of the scatter projections are over-estimated. These regions correspond to pixels that are more strongly influenced by activity voxels included in the 15% threshold. Similarly regions that would be influenced by voxels which were not included in the threshold are under-estimated. For the noiseless data, the OSEM-APDI(0%)P contrast fell somewhere between the OSEM-ISR and OSEM-APDt values. Since OSEM-APDt uses the true scatter projec-tions, this implies that the APDI component of scatter is over-estimated in detector pixels which are significant to the defect contrast. In the first scatter iteration of OSEM-APDI(0%)p, the APDI scatter estimate is generated using the OSEM-BB activity dis-tribution. This OSEM-BB image has lower defect contrast which corresponds to higher levels of activity in the defect region. Since scatter is "forward-peaked" (i.e., photons are more likely to undergo small angle scatter) scatter is over-corrected in detector pixels influenced by this false activity in the OSEM-BB defect region. This over-correction was apparently only slightly reduced by performing further scatter iterations. The liver NSD values for the noiseless projection data shown in figure 4.16 (b) demon-strate that all reconstructions reach a minimum level of variation for uniform regions at seven or eight OSEM iterations. The OSEM-APDI(0%)P reconstructions demonstrated the lowest minimum NSD while the OSEM-ISR reconstruction shows the largest mini-mum NSD. The OSEM-APDt, OSEM-BB, OSEM-APDI(15%) and OSEM-APDI(0%)P Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 107 reconstructions all have minimum NSD values within 1% of each other. It is unclear why OSEM-ISR fared more poorly than the other reconstructions. However, all recon-structions demonstrated very similar minimum NSD values. All NSD values fell between 0.072 and 0.066 indicating very similar image uniformity for all reconstructions. This means the activity in the liver VOI for all reconstructions had a standard deviation that was about 7% of the mean. It is perhaps more meaningful to compare NSD values for projection data that con-tained statistical noise. In this case, the NSD values provide an indication of how sta-tistical noise in the projection data propagates in the reconstruction. The poorest NSD values were for the OSEM-ISR reconstruction. This is because the projection data used by OSEM-ISR only contained unscattered photons and therefore had higher statistical noise than the combined scattered and unscattered data. The two versions of the OSEM-APDI reconstruction had lower or the same NSD values as OSEM-ISR, due to the way in which the scatter correction was included in these reconstruction. If the APDI scatter projections were directly subtracted from the measured data, the reconstructed image noise would increase above the level in OSEM-ISR. Although an explicit subtraction is not performed, the inclusion of scatter in the OSEM-APDI reconstruction still con-tributes to the noise in the image. The OSEM-APDI(0%)P reconstructions had slightly lower NSD values than OSEM-APDI(15%), probably due to the use of the thresholded activity map and hence biased APDI scatter estimates. Increasing the number of scatter iterations in OSEM-APDI(0%)P did not significantly change the NSD values. In all cases the OSEM-BB reconstruction resulted in the lowest NSDs. However, the NSD values for both versions of OSEM-APDI were not any higher than the ideal case of OSEM-ISR. The reduced contrast of the OSEM-BB reconstruction is clearly visible in the short axis slices shown in figures 4.17 and 4.20. Another indication of image contrast can be Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 108 calculated using the image profiles shown in figures 4.18 and 4.20. The ratio of the min-imum depth of the profile in the defect region relative to the maximum height of the myocardium region was compared. A higher ratio indicates poorer contrast between the two regions. The ratio for the OSEM-BB reconstruction for both the noiseless and Sim-SET projection data was about 200% greater than either the true MCAT phantom or any of the other reconstruction techniques. This ratio differed by at most 20% between each of the reconstructions that used an accurate scatter correction (i.e., OSEM-APDI(15%) OSEM-APDI(0%)P, OSEM-ISR and OSEM-APDt). Visually it is difficult to distinguish between the OSEM-APDI(15%) OSEM-APDI(0%)P, OSEM-ISR and OSEM-APDt short axis slices. However, for both sets of projection data the OSEM-BB short axis slices ap-pear different from the other reconstructions. The OSEM-BB image shows a false increase in uptake in the right ventricle, inferior wall and the tissue between the heart and the liver. Scatter from the viable heart tissue into the defect region also contributed signifi-cantly to the observed loss in contrast. These effects were reduced greatly in the scatter corrected data. 4.3.7 C O N C L U S I O N S As described in the introduction, the presence of statistical noise in the projection data causes image noise in MLEM based reconstructions (such as OSEM) to increase with the number iterations. This effect can be observed in the constantly increasing NSD values in figure 4.16 (d). However, for very low OSEM iterations the image contrast is extremely poor. The use of ordered subsets also requires that the number of OSEM iterations be evenly divisible by the the number of angular views acquired. This requirement suggested that the best compromise between image contrast and image noise is eight OSEM iterations for the all the reconstruction techniques considered. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 109 There seems to be little need to perform additional scatter calculations beyond the initial APDI calculation using the OSEM-BB estimate of the activity distribution. For the image quality measures considered here, the first scatter iteration of the OSEM-APDI(0%)p algorithm was not very different from the sixteenth scatter iteration. The timing results in table 4.4 demonstrate that thresholding the activity map reduced the CPU time required to perform the first APDI calculation by about a factor of three. Thresholding appears to introduce some systematic errors in the APDI scatter estimates when compared with the true scatter projections. However, these errors did not signif-icantly degrade the image quality of the OSEM-APDI(15%) reconstruction relative to the ideal case of OSEM-ISR. Based on these observations, it was concluded that OSEM-APDI(15%) provided the best compromise between image quality and the required CPU time in this particular case. If the CPU time is ignored, the best reconstruction was OSEM-APDI(0%)3 in terms of convergence and accuracy of the APDI scatter estimates. The OSEM-APDI(15%) reconstruction was selected for the studies discussed in the fol-lowing section because it provided results that were closer to the ideal case of OSEM-ISR in terms of image quality and required much shorter computation times than OSEM-APDI(0%)3. 4.4 CORONARY PERFUSION SIMULATIONS The following section evaluates the consistency of the OSEM-APDI(15%) algorithm using the noiseless and MC simulated projection data for simulated phantoms with perfusion defects in a variety of positions. The simulated activity distributions used to generate the projections in this section were the MCAT thorax phantoms with defects in the anterior, apical, inferior, lateral and septal walls of the myocardium. The OSEM-APDI(15%) reconstructions were compared to the OSEM-ISR, OSEM-BB and for noiseless projection Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 110 data the OSEM-APDt reconstructions. The goal of this work was to verify that the OSEM-APDI(15%) consistantly provided accurate reconstructions for perfusion defect studies regardless of the defect position. 4.4.1 R E S U L T S The reconstructed images were compared with respect to the image quality measures described previously. In this analysis the contrast between the coronary perfusion defect and the healthy portion of the myocardium was determined for each defect position. The contrast for each of the defect positions are displayed for the SimSET projection data in figure 4.21 (a) and for the noiseless projection data in figure 4.21 (b). In this analysis NSD values for both the liver and body VOIs were considered. The body VOI was defined as only those voxels in the original MCAT phantom that were assigned the activity concentrations of 1.0 juCi/cc. Therefore, this VOI did not include any lung tissue or voxels influenced by the partial volume effect. The noiseless and MC simulated projection data used in this section consisted of the sum of the three separate simulations generated using the components of the MCAT phantom for the defective myocardium, liver, and body (lungs and remaining tissue). The same set of simulated projection data for the liver and body component of the MCAT phantom was used for each different defect position. It was found that there were negligible differences in the NSD values for activity distributions with different coronary defect positions. This suggests that the activity distribution in the myocardium has minimal influence on the liver and body uniformity. These NSD values are summarized for the noiseless and SimSET data in table 4.5. Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 111 Anterior Apex Inferior Lateral Septal Anterior Apex Inferior Lateral Septal Defect Position Defect Position (a) (b) Figure 4.21: Comparison of contrast between the healthy myocardium and defect VOIs for each of the five coronary defect positions using the (a) the noiseless and (b) the SimSET projection data. 4.4.2 D I S C U S S I O N The contrast values for the noiseless projection data shown in figure 4.21 (a) were also consistent with the results of the previous section, regardless of position of the coronary defect position with only one exception. The OSEM-APDI(15%) reconstruction had lower contrast than O S E M - A P D t for the anterior wall defect. In all other cases, the OSEM-APDI(15%) had contrast between the O SEM- A PD t and OSEM-ISR results. The maximum difference between the OSEM-APDI(15%) and OSEM-ISR contrast was 6% for the anterior wall defect. The average percent difference between the OSEM-APDI(15%) and OSEM-ISR contrast was 2%. In comparison, the average percent difference between the O S E M - B B and OSEM-ISR contrast was 22% with a maximum difference of 29%. Similar results were obtained for the SimSET projection data and are shown in fig-ure 4.21 (b). These results were also consistent with the previous section with a single exception. In the case of the septal wall defect, the OSEM-APDI(15%) was somewhat Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 112 Table 4.5: Comparison of NSD image quality measures for the liver and body VOIs for each reconstruction technique using the noiseless and SimSET projection data. Noiseless Liver Body NSD NSD SimSET Liver Body NSD NSD OSEM-ISR OSEM-APDI(15%) OSEM-APDt OSEM-BB 0.072 0.242 0.069 0.242 0.067 0.242 0.067 0.223 0.20 0.30 0.19 0.28 N/A N/A 0.17 0.26 higher than the OSEM-ISR reconstruction. For each of the other defect positions, the contrast values for OSEM-APDI(15%) were slightly lower than that of the OSEM-ISR reconstruction. The average percent difference between OSEM-APDI(15%) and OSEM-ISR was 3% with a maximum difference of 8%. The average percent difference between OSEM-BB and OSEM-ISR contrast was 25% with a maximum difference of 35%. The NSD values for both the liver and body VOI summarized in table 4.5 demonstrate that for the noiseless projection data the OSEM-BB, OSEM-APDt, OSEM-APDI(15%) and OSEM-APDI(0%) reconstructions all resulted in similar NSD values in both VOIs. The OSEM-ISR resulted in the highest NSD, as in the previous section. The NSD values for the SimSET projection data were also consistent with the results of the previous section. The OSEM-BB reconstruction had the lowest NSD values, OSEM-ISR had the highest, and OSEM-APDI(15%) fell somewhere between the two extremes. 4.4.3 C O N C L U S I O N S The OSEM-APDI(15%) reconstruction algorithm was shown to behave consistently for most of the defect positions considered in this analysis. This reconstruction does not sig-nificantly degrade image uniformity relative to the ideal case of OSEM-ISR and greatly Chapter 4. Validation of APDI and OSEM-APDI Using Computer Simulations 113 improves image contrast. The poor contrast associated with the OSEM-BB could have serious repercussions in the quantification of the size and level of perfusion in cardiac defects. The OSEM-APDI(15%) reconstructions had significantly improved image con-trast relative to OSEM-BB. The OSEM-APDI(15%) reconstruction consistently produced results which were qualitatively and quantitatively very close to the ideal case of OSEM-ISR. C H A P T E R 5 EXPERIMENTS The following chapter1 describes the experimental validation of the OSEM-APDI recon-struction algorithm. This chapter is divided into three sections. Before the OSEM-APDI reconstruction method can be applied to experimentally acquired projection data, one must generate the APD look-up tables used in the APDI scatter projection calculation. These look-up tables depend on the patient-independent characteristics of the SPECT imaging system being used. Section 5.1 describes how the necessary characteristics were determined for a Siemens E-Cam camera system. Once these look-up tables were generated, it was then possible to perform the OSEM-APDI reconstruction algorithm using data acquired on this camera. Section 5.2 describes an experiment that employed a physical phantom to realistically model the activity and attenuation distributions in the thorax of a patient with a myocardial perfusion defect. The phantom was injected with known quantities of radioactivity and was then used X A version of this chapter has been published. E. Vandervoort, A. Celler, R. G. Wells, S. Blinder, K. Dixon and Y. Pang. Implementation of an Analytically Based Scatter Correction in SPECT Recon-structions. In 2003 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, pages M10-256, October, 2003. A version of this chapter has also been submitted for publication. E. Vandervoort, A. Celler, R. G. Wells, S. Blinder, K. Dixon and Y. Pang. Implementation of an Analytically Based Scatter Correction in SPECT Reconstructions. IEEE Trans. Nucl. Sci, 2003. 114 Chapter 5. Experiments 115 to acquire experimental SPECT projections. These projections were employed in the OSEM-APDI reconstruction algorithm. The accuracy of OSEM-APDI was evaluated by comparing the APDI data with the experimental projections. The properties of the OSEM-APDI and OSEM-BB reconstructions were compared to the known activity con-centrations in the phantom. Section 5.3 discusses the application of the OSEM-APDI scatter correction to clini-cally acquired data. A myocardial perfusion study was selected for a patient that demon-strated an extremely elevated liver. In this analysis the APDI data were also compared to the measured data. However, since one cannot know the true distribution of activity in a patient study, the OSEM-BB and OSEM-APDI images could only be compared qualitatively in this case. 5.1 GENERATION OF THE EXPERIMENTAL A P D L O O K - U P TABLES The purpose of the work performed in this section was to determine the characteristics of the Siemens E-Cam gamma camera, which are necessary for the calculation of the look-up tables used in the APDI and APD scatter calculations. The patient-independent parameters on which these look-up tables depend include the energy resolution of the detector, the energy of the emitted photon, the Klein-Nishina scattering cross section, the distance dependent spatial resolution of the collimator, and the intrinsic resolution of the detector. The energy response of the camera system is modeled in the look-up table calculation with the function PE{8) and the geometrical collimator response is accounted for by F(ÂŁ'). These functions can be determined either experimentally or modeled theoretically. It should be stressed that the only variables that change between the APD look-up table calculations for different SPECT imaging systems are geometrical variables; such as Chapter 5. Experiments 116 the size of the detector pixels and object voxels, and the two functions: F(ÂŁ') and P.E(0). The APD method was extensively validated for experimental data (e.g., [104]) and is extremely accurate providing that the F(f') and P (^f9) functions used in the look-up table calculation realistically model the collimator and energy response of the detector under investigation. The function F(ÂŁ') represents the probability that a photon is transmitted through the collimator if the angle between the photons direction of propagation and a vector perpendicular to the detector surface is In the original work performed by Wells et al. [104, 3], the F(ÂŁ') function was determined experimentally through the analysis of projections acquired from a point-like sources positioned at various distances from the camera surface. These experimental F(ÂŁ') values were compared to a theoretical model [94, 8], but the final form of the F(ÂŁ') function used in the APD look-up tables was a seventh order polynomial fit to the experimental values. Since the differences between the theoretical F(f') and the function fit to the exper-imental F(ÂŁ') values were extremely small2, the theoretical form of F(ÂŁ') was used here. This theoretical model assumes a parallel hole collimator with circular holes, however the physical collimator had hexagonal holes (shown schematically in figure 5.22). The hole diameter was selected to be the average of the width and height of the hexagonal hole. From the manufacturers specifications, a circular-hole diameter of 0.118 cm, a collimator thickness of 2.36 cm, and a septal thickness of 0.016 cm were used. As in the original work by Wells et al., the P.e(0) was modeled as a Gaussian with a FWHM equal to 10% of the incident photon energy and a 20% energy window. In order to evaluate the accuracy of the F(ÂŁ') and PB(0) functions, which were later used in the generation of the APD look-up tables, experimental projection data were 2See for example figure 1 in [104] Chapter 5. Experiments 117 Septal Thickness Figure 5.22: Schematic diagram of a hexagonal collimator hole. analyzed for line sources in air (i.e., with negligible photon scatter and attenuation). In this case, one cannot directly determine the F(f') probability values, however its accuracy can be indirectly evaluated through its influence on the primary photon distribution. In the absence of an attenuating medium, the APD (and hence APDI) primary photon distribution for detector bin i, PPD(A, n,) is given by: PPD(A,n,) = ÂŁ W O ) f d 2 n F ( 0 ^ , (5.42) where all variables are defined as in chapter 3. In this case, PPD(A,n;) only depends on: the position and intensity of the activity source, P.E(0) the energy response of the detector for the unscattered photon energy, and the integral of F(f'). Throughout this discussion the distribution described by equation 5.42 shall be referred to as the APDI primary photon distribution. However, it should be noted that in the APDI method no approximations were made to the original APD primary calculation. As a result, there is no difference between the APD and APDI primary photon distributions. Another topic that was investigated here, was the accuracy of a method which was used to increase the matrix size of the generated APDI projections. In clinical SPECT Chapter 5. Experiments 118 scans, projection data is typical acquired in a 128x128 matrix size. In order to perform the OSEM-APDI scatter correction, one must use scatter projections which have the same dimensions as the acquired data. The APDI projection calculation used in the previous chapter was only capable of generating 64x64 projections. It was therefore beneficial to investigate the use of expanded 128x128 projection matrices. The larger projection matrix size was obtained by performing the original APDI calculation to generate 64x64 projections and then using bi-linear interpolation to increase the size to 128x128. The line source experiments used here were performed by Blinder et al. and were presented in the 2001 IEEE Nuclear Science Symposium and Medical Imaging Conference Record [107]. These experiments used capillary tubes that were each filled with the same amount of 9 9 m T c . The tubes were 1.11 mm in diameter and contained activity over a length of 60.0 Â± 0.5 mm. The scan was performed using an dual head Siemens E-Cam camera equipped with a Low Energy High Resolution (LEHR) collimator. The camera heads were positioned at 180Â° to each other and a total of six tubes were placed in Styrofoam supports at various distances between the two detector heads. Two line sources were placed in each support, one source oriented along the vertical and another along the horizontal axis of the detector surface. Projections were collected using a 1024x1024 matrix size with a 0.599 mm pixel length. The projections were acquired for 120 s resulting in approximately 2.2xlO6 counts per projection. In order to generate the APDI projections, it was necessary to create an activity map that approximated the line source distributions which where used in the experiment. The APDI calculation was performed using a 64x64x64 activity map and a 64x64 projection matrix size with 9.590 mm pixel and voxel lengths. The line source length in all cases was 57.54 mm (6 pixel lengths). The linear attenuation coefficient of air was assumed to be zero. In the APDI calculation, all photons are assumed to be emitted from the center Chapter 5. Experiments 119 of each voxel. The position of each activity voxel is therefore an integral number of pixel lengths from the surface of the collimator. As a result, the distance from the surface of the collimator to each line source in the APDI calculation could only be matched to the experimental distance to within Â± 4.759 mm (i.e., half a pixel length). The APDI projection data was then expanded using bi-linear interpolation to a 128x128 matrix size with 4.759 mm sized pixels. It was also necessary to correct the APDI projection data to account for the intrinsic resolution of the gamma camera. As described in chapter 1, the spatial localization of a detected photon in the scintillating crystal is achieved by measuring the output signal from an array of photomultiplier tubes (PMT) that amplify the scintillation light produced in the crystal. Both the conversion of photon energy to light and the PMT amplification are statistical processes that cause an intrinsic limitation in the camera's spatial resolution. The intrinsic spatial resolution can be accounted for by convolving the two dimensional array of PPD(A, Hi) values with a function, frequently a Gaussian, that describes the intrinsic resolution of the camera (see for example [8, 63, 104, 107]). Previous quality control studies have shown that the intrinsic resolution for this imaging system can be accurately modeled as a Gaussian function with FWHM equal to 4.51 mm. The convolution was carried out on the 64x64 APDI projection matrix prior to expanding to the 128x128 size, in order not to enhance a blurring effect caused by the bi-linear interpolation. This effect will be discussed in more detail subsequently. 5 . 1 . 1 R E S U L T S Figures 5.23 (a) and 5.24 (a) show the experimental projections for the first and second camera heads, respectively. The projections have been reduced to a 128x128 projection matrix size to allow them to be compared to the APDI projections. The reduction Chapter 5. Experiments 120 was performed by adding 8x8 pixels in the 1024x1024 matrix to form one pixel in the 128x128 matrix. Figures 5.23 (b) and 5.24 (b) show the expanded 128x128 APDI projections projections for the first and second camera head respectively. Approximately 5 to 7 profiles were extracted for each line source from the reduced 128x128 experimental projections. Each profile was fit with a Gaussian function using unconstrained nonlinear optimization based on the Nelder-Mead Simplex search method [110]. In the same way, 5 to 7 profiles through the 128x128 APDI projections were fit for each line source. Table 5.6 summarizes the average FWHM for each source distance, determined for both the vertical (FWrIM s) and horizontal (FWrlMj) profiles. The error quoted for each FWHM value is the standard deviation over all the profiles fit. Table 5.6: Comparison of mean FWHMs fit to profiles through both experimental and APDI projection data. The APDI FWHMs were fit to projection data which was generated using 64x64 APDI projection matrix and expanded to 128x128 using linear interpolation. The experimental FWHMs were fit to the 128x128 projection matrix size, which had been reduced from 1024x1024, to allow for its comparison with the APDI projections. Line Source APDI F W H M X APDI FWHMy Experimental Experimental Distance (cm) (mm) (mm) FWHM x(mm) FWHMy (mm) 0.00 9.21 Â± 0.00 9.21 Â± 0.00 5.56 Â± 0.08 5.69 Â± 0.04 10.15 9.40 Â± 0.00 9.40 Â± 0.00 8.54 Â± 0.14 7.93 Â± 0.04 20.38 12.42 Â± 0.04 12.42 Â± 0.03 12.34 Â± 0.13 12.21 Â± 0.08 39.22 21.04 Â± 0.06 21.03 Â± 0.06 20.57 Â± 0.20 20.07 Â± 0.14 49.45 25.24 Â± 0.12 25.26 Â± 0.13 25.38 Â± 0.40 24.92 Â± 0.15 59.60 29.41 Â± 0.11 29.45 Â± 0.11 30.34 Â± 0.41 29.59 Â± 0.25 5.1.2 D I S C U S S I O N The average FWHM of the Gaussian functions that were fit to the A P D I projection profiles were within 3% of the those fit to the experimental data for all but the line sources located at 10.15 cm and 0.00 cm. For these line sources distances, the profiles Chapter 5. Experiments 121 20 40 60 80 100 120 (a) 20 40 (b) 80 100 120 Figure 5.23: Comparison of the experimentally acquired projection data (a) and the expanded 128 x 128 APDI photon distribution (b) for the line source experiment with sources located at 0, 10.15 and 20.38 cm from the collimator surface. The experimental projections have been reduced to a 128x128 projection matrix size to allow them to be compared to the APDI projections. 20 40 60 OO 120 (a) 6000 20 5000 40 4000 60 3000 80 2000 100 1000 120 20 40 60 80 100 120 (b) Figure 5.24: Comparison of the experimentally acquired projection data (a) and the expanded 128x128 APDI photon distribution (b) for the line source experiment with sources located at 39.22, 49.45 and 59.60 cm from the collimator surface. The experimental projections have been reduced to a 128x128 projection matrix size to allow them to be compared to the APDI projections. Chapter 5. Experiments 122 Figure 5.25: Schematic diagram to illustrate an error introduced when bi-linear interpolation is used to expand from a 64x64 to a 128x128 matrix size. Figure (a) shows a profile extracted from the original 64x 64 matrix with pixel length Ax. This profile represents an impulse function of height h. Figure (b) shows the resulting function when linear interpolation is used to increase the grid sampling by a factor of two (i.e., pixel length = ^p). The resulting triangular function has a FWHM equal to Ax. through the original 64x64 APDI projection matrices were essentially equal to an impulse function (i.e., these profiles contained only one pixel with a significant number of detected photons). The FWHMs for these distances reflect inaccuracies introduced by the bi-linear interpolation method that was used to expand the projection matrix size. This effect is shown schematically in figure 5.25. The minimum FWHM of any profile through a 128x128 matrix obtained using linear interpolation on a 64x64 matrix is equal to the original pixel length, in this case 0.959 cm. The FWHMs fit to APDI data, for the 10.15 cm and 0.00 cm source distances, were both approximately equal to 0.959 cm. The differences between the FWHMs and the pixel length are probably due to difficulties associated with fitting a Gaussian to a triangular function. To evaluate the significance of errors introduced in the expansion of the APDI ma-trix size, the FWHMs of 3 to 4 profiles through the original 64x64 APDI projections were fit for each line source. It was expected that profiles through the 64x64 projection data should have approximately the same FWHM as the expanded 128x128 data. The 64x64 Matrix Size 128x128 Matrix Size Linear Interpolation >-Ax Minimum Ax FWHM y = Ax Z Chapter 5. Experiments 123 FWHMs for line sources at 10.15 cm and 0.00 cm were not considered in this analysis because the profiles extracted from the 64x64 matrices were essentially impulse func-tions and therefore could not be reliably fit with Gaussian functions. The results are summarized in table 5.7. Table 5.7: Comparison of mean FWHM fit to profiles through APDI line source projection data for the original 64x64 and the expanded 128x128 matrix size. Line Source Distance (cm) APDI 64x64 matrix FWHM X FWHM y (mm) (mm) APDI 128x128 matrix FWHM X FWHMy (mm) (mm) 20.38 39.22 49.45 59.60 12.51 Â± 0.08 12.50 Â± 0.07 19.78 Â± 0.08 19.78 Â± 0.08 24.08 Â± 0.15 24.08 Â± 0.15 28.53 Â± 0.16 28.56 Â± 0.14 12.42 Â± 0.04 12.42 Â± 0.03 21.04 Â± 0.06 21.03 Â± 0.06 25.24 Â± 0.12 25.26 Â± 0.13 29.41 Â± 0.11 29.45 Â± 0.11 The differences between the FWHMs fit to profiles through the original APDI 64x64 sized matrices were within 6% of the FWHMs of the expanded 128x128 matrices. These differences may be due to difficulties associated with fitting a Gaussian to the profiles from the smaller 64x64 sized matrix. In this case there were fewer profiles and fewer sampled data points to fit. Another way to evaluate the significance of the differences observed between the var-ious FWHMs measured in this section is to compare those differences with the length of a pixel. The pixel length is a limiting factor in how well the features of the projection data can be measured. Therefore, any differences between the APDI and experimental FWHMs that are considerably less than a pixel length are not significant. For the com-parison of the 128x128 experimental projections with the 128x128 APDI projections, the maximum difference between the experimental and APDI FWHMs was 0.97 mm, ex-cluding the 0.0 cm and 10.15 cm line sources distances. The pixel length for the 128x128 Chapter 5. Experiments 124 matrix size is 4.795 mm, hence the maximum difference is only 21% of the pixel length. In the comparison of the 128x128 and 64x64 APDI projections the limiting length is the 64x64 pixel size (9.59 mm). The maximum difference observed here was 1.3 mm, which is only 21% of the 64x64 pixel length. 5.1.3 C O N C L U S I O N S The small size of the differences between the experimental and APDI FWHMs, indicate that P^(0) and F(ÂŁ') functions accurately reflect the energy resolution and geometrical response of the Siemens E-Cam LEHR camera. As a result, the APD look-up tables generated using these functions are expected to be appropriate for this SPECT imaging system. The APD (or APDI) projections that are generated using these look-up tables should therefore accurately reflect the true distribution of scattered and unscattered photons. The expansion of the original 64x64 APDI projection matrix to the 128x128 matrix size does not introduce significant errors to the projection matrices, except in the case where the data contains features that change rapidly over a length of about 1 cm. This effect was observed for sources at small distances from the detector (less than about 10 cm). This is a fundamental limitation connected to the pixel length used in the APDI calculation (0.959 cm) and the bi-linear interpolation method used to expand the matrices. These small differences are not expected to have a significant affect on the calculation of APDI scatter projections, since these distributions are extremely smooth and therefore do not change very rapidly over small spatial distances. 5.2 P H A N T O M EXPERIMENTS This section describes the application of the OSEM-APDI reconstruction algorithm using experimental SPECT projection data. The data was acquired using a physical phantom Chapter 5. Experiments 125 that realistically models the distribution of activity and attenuation in the thorax of a patient wi th a myocardial perfusion defect. The thorax phantom used in this experiment consisted of an ell iptical cylinder wi th inserts for the lungs, spine and heart. The heart insert had a fillable defect that was positioned in the inferior wall of the myocardium. The phantom was injected wi th measured quantities of 9 9 m T c and experimental projection data was then acquired. The experimental projection data was reconstructed using the O S E M - A P D I reconstruction algorithm. This analysis could not be performed in exactly the same manner as chapter 4. The accuracy of the A P D I scatter estimate cannot be evaluated directly, since the true com-ponent of scatter in the acquired projection data is not known. Therefore, the sum of the primary and scattered A P D I projections was compared to the experimental data, which contain the contribution from both scattered and unscattered photons. Similar ly the O S E M - A P D t and O S E M - I S R reconstructions are not possible wi th experimental data, since they both require the true scattered photon data. Therefore, only the O S E M - A P D I and O S E M - B B reconstructions were compared in this analysis. Since there was no ideal scatter corrections wi th which to compare the O S E M - A P D I reconstruction, the image quality was more difficult to evaluate using experimentally ac-quired data. However, each of the simulated organs in the thorax phantom had known volumes and were injected wi th measured amounts of 9 9 m T c . Therefore, the image con-trast of the O S E M - A P D I and O S E M - B B reconstructions could be compared to the true contrast calculated from the known activity concentrations in each of the phantoms or-gans. Another difficulty arises because the true 3D distribution of activity is not available to be used to define VOIs . Errors are introduced because the VOIs can only be estimated using the reconstructed images. Chapter 5. Experiments 126 5.2.1 T H E E X P E R I M E N T The thorax phantom used in this experiment was an Elliptical Lung-Spine Body Phantomâ„˘ with a Cardiac Insertâ„˘ and a Fillable Defect Setâ„˘ from Data Spectrum Corporation. The combined cardiac and defect insert create three isolated chambers that represent the defect, left ventricle and heart chamber and each chamber could be filled with different activity concentrations. The defect was positioned in the inferior wall of the myocardium. The phantom contained left and right lung inserts filled with Styrofoam beads to simulate the lower attenuation of lung tissue in the human body. The phantom also contained a cylindrical TeflonÂ® rod insert that modeled the higher attenuation of the spine. A l l other organs were filled with water to simulate the attenuation of soft tissue. Two water bags (each Â« 2 liters) containing no activity, were strapped to the anterior of the phan-tom to simulate the presence of mammary tissue. Mammary tissue enhances the effects of both attenuation and scatter. The phantom was injected with measured quantities of 9 9 m T c to approximate 5 times the activity of a typical clinical injection. In a typical coronary perfusion study, the patient is injected with 50 MBq of 9 9 m T c , of which approximately 10% is taken up by the myocardium tissue. Two SPECT acquisitions were performed. In the first scan the data was acquired for 20s per projection angle (high statistics) while the second scan had a 4s acquisition (normal statistics). The volumes, total injected activity and activity concentrations for each of the phantoms organs are summarized in table 5.8. The concentrations have been corrected for the radioactive decay that occurred during the time between measuring the activity and midway through each scan. The activity injected in the remaining volume of the cylinder and the heart chamber were selected to give similar relative concentrations as those observed clinically. The defect, lungs and spine did not contain any activity. Chapter 5. Experiments 127 Table 5.8: Summary of volumes, injected activities, and activity concentrations corrected for decay for both high and normal statistics scans for each of the phantoms simulated organs. Organ Injected Volume High Statistics Normal Statistics Activity (mL) Concentrations Concentrations (MBq) (kBq/mL) (kBq/mL) Healthy Myocardium 26 88.0 261.5 257.3 Defect 0 11.4 0 0 Heart Chamber 4.2 62.0 60.6 59.6 Remaining Cylinder 274 5863 42.9 42.3 The data was acquired using Siemens' standard cardiac acquisition protocol, on a dual-head Siemens E-cam camera equiped with LEHR collimators and a Profile trans-mission system. The two heads were set at ninety degrees to each other. Emission and transmission data were acquired simultaneously. A total of 64 projection angles evenly spaced over 180Â° were acquired. The detectors were rotated counter-clockwise over the anterior side of the phantom between the left posterior oblique (-45.00Â°) to the right anterior oblique position (135Â°). The data was acquired in "stop and shoot" mode with a 128x128 projection matrix size and 0.4795 cm bins. The detectors were rotated using a non-circular orbit (NCO), in which the acquisition software shifts the detector as close as possible to the phantom at each projection angle. 5.2.2 OSEM-APDI RECONSTRUCTION SPECIFICATIONS The method used to generate APDI projections in the OSEM-APDI reconstruction dif-fered slightly from that used in chapter 4. In experimental studies third and higher order Compton and Rayleigh scattered photons must be considered. As a result the complete version of the total scatter distribution function TSDF(s"j, rfj), given in equation 3.20, was used throughout this chapter. Rayleigh and Compton scattering form factors were also Chapter 5. Experiments 128 included in the APD look-up table calculations. Also new to this chapter is the inclusion of the effects of the intrinsic resolution of the detector system and the expansion of the APDI projection matrices from their original size of 64x64 to 128x128. The procedure for the generation of APDI scatter estimates, as part of the OSEM-APDI reconstruction, is shown schematically in figure 5.26. The APDI calculation re-quires an activity map and an attenuation map to generate scatter projections. In this case, the activity map used was the OSEM-BB reconstruction. The attenuation map used was reconstructed from the SPECT transmission scan using Siemens' software. Previous experiments have shown that maps acquired from larger patients require segmentation to improve image quality. In this study all maps were segmented for consistency. The attenuation and activity maps both consisted of a series of 128 x 128 sized transax-ial slices. A Gaussian based diffusion rotation algorithm [108] was used to rotate the OSEM-BB reconstruction and the attenuation map. The maps were then shifted toward or away from the detector to account for the NCO that was used to acquire the pro-jections. The APDI calculation requires 64x64x64 sized attenuation and activity input files. Slices containing only zeros were added to the activity and attenuation maps to make them both 128x128x128 sized volumes. The activity map was then reduced to 64x64x64 by adding 2x2x2 voxels in the 128x128x128 matrix to form one voxel in the new activity map. The average linear attenuation coefficient of the 2x2x2 group of voxels was used to determine the value for a single voxel in the new attenuation map. The rotated, shifted and smaller sized attention map was then used to generate a map of electron density values and the evaluated attenuation map. The OSEM-BB reconstruction was thresholded to include all voxels greater than 1% of the maximum count. This threshold was selected based on visual inspection of the image. When this threshold was applied to the OSEM-BB data, nearly all voxels that Chapter 5. Experiments 129 appear to be within the elliptical cylinder were included. This threshold removed those voxels outside the cylinder and a number of voxels inside the lung regions. It was known a-priori that these voxels did not contain activity, and were therefore an artifact of the OSEM-BB reconstruction. This lower threshold should give better accuracy than the 15% threshold used in the simulations and did not significantly increase computation time. The thresholded OSEM-BB activity map, attenuation map, electron density map, and evaluated attenuation map were input into APDI to generate the primary and scat-tered photon projections. The 64x64 sized APDI projection matrices were convolved with a Gaussian function to account for the intrinsic resolution of the detector and then expanded using bi-linear interpolation. The final output was 128x128 APDI projection matrices that were then used in the OSEM-APDI reconstruction algorithm. This recon-struction and the scatter estimation step shall be referred to as OSEM-APDI(l%) and APDI(1%) respectively, throughout the following discussion. 5.2.3 R E S U L T S 2 The -fy and bias measurements could not be evaluated directly using the scatter pro-jections generated in the APDI(1%) step of OSEM-APDI(l%). In order to get some estimate of the accuracy of the APDI scatter projections, the fy and bias measure-ments were performed using the experimental projections, which obviously contain the contribution from both scattered and unscattered photons, and the sum of the APDI primary and scattered photon distributions. Figure 5.27 shows the fy and bias values as a function of projection angle for these experiments. Figures (a) and (b) are for the high statistics projection data while (c) and (d) are for the normal statistics data. Table 5.9 2 summarizes the mean fy and bias values as well as their standard deviations over all Chapter 5. Experiments 130 Inputs: OSEM-BB Activity Map 128x128 Narrow-Beam Attenuation Map 128x128 Rotate Shift for in 2D n.c.o. Reduce to 64x64 ^uOOdx Evaluated Attenuation Map Create Electron Density Map APDI Scatter Pi ejection Calculation Convolve by Intrinsic Resolution Expand to 128x128 64x64 Projections Output: Primary and Scatter 128x128 Projections Figure 5.26: Schematic diagram to illustrate the steps involved in the APDI projection calcu-lation for the data used in this chapter. The figure shows all the intermediate steps which are required to accurately model clinical data acquisition. projection angles. Table 5.9: Comparison of mean fj and the mean bias, averaged over all 64 projections, and their standard deviations for the high and normal statistics phantom experiments. The experimentally acquired projection data was compared to the sum of the primary and scattered APDI(1%) photon distributions. Standard Deviation Bias Standard Deviation High Statistics Normal Statistics 1.75 0.32 0.77 0.16 0.96 0.02 0.98 0.02 To evaluate the contrast and NSD image quality measurements, VOIs that accurately reflect the position of the different organs in the phantom must be determined. In the simulation studies the true pixelized activity distribution was available to create these VOIs. This is not possible in experimental studies. In particular, it is difficult to Chapter 5. Experiments 131 1.2 -50 j â€” Inferior Wall Defect High Statistics | fl : 1 Â» 1 i i Â« i i i '\ ' i . ' â€˘ " n .' 1 i . " " i 1 ' â€˘ . " > ' i 1 ' â€˘ , ' ' . ' 1 i i . , 1 ' / * ' â€˘ â€˘ i ' ' (a) 50 100 Projection Angle (degrees) 1 0.99 0.98 0.97 Â« ÂŁ 0.96 CO 0.95 0.94 0.93 150 0.92 | Inferior Wall Defect High Statistics | i i f \ i i* i i t ! 'K' i , i i i > I #* i ft 1 Â» - - ' i * *â€˘ i / * 14 ' '1 ' >' -50 0 50 100 Projection Angle (degrees) 150 (b) 0 50 100 Projection Angle (degrees) (c) 150 0 50 100 Projection Angle (degrees) 150 (d) Figure 5.27: Plot of the fy and bias values as a function of projection angle for the phantom 2 experiments. The fy values are shown for the high and normal statistics projections in figures (a) and (c), respectively. The bias values are shown for the high and normal statistics projec-tion data in figures (b) and (d), respectively. The experimentally measured projections were compared to the sum of the primary and scattered APDI(1%) projection data. Chapter 5. Experiments 132 classify voxels in the reconstruction as belonging to one organ or another. Therefore, VOIs representing the healthy myocardium and defect were determined directly from SPECT images using the image segmentation software iQuant [111]. iQuant is an image analysis toolkit created on a Matlab platform. Thresholds, defined as a percentage of the maximum count in the myocardium, are used to define viable myocardial tissue and to provide a rough outline of the complete myocardium. This outline together with the visible thickness of the viable myocardium is used as a guide to determine the location and size of defective tissue. The iQuant VOIs were determined from the OSEM-APDI reconstruction. A VOI representing the background region in the phantom was obtained by manually selecting polygonal regions from the reconstructed slices. The background VOI represents what has been referred to as the remaining cylinder. In table 5.10 the "true" image quality measures, that were calculated using the known activity concentrations, are compared with values measured using the reconstructions. Reconstructions using both the high and normal statistics projection data are considered. The true calculated contrast values, determined using activity concentrations that were corrected for radioactive decay, were negligibly different for the high and normal statistics scans. The amounts of CPU time that were necessary to perform the OSEM-APDI(l%) reconstruction, using a computer with a 1.7 GHz processor and 1024 Mb of RDRAM, are summarized in table 5.11. The CPU time was split up into three parts. The first part was the amount of time necessary to perform the initial OSEM-BB reconstruction. The second part was the total time required to perform all 64 APDI(1%) scatter projection calculations. The final part was the amount of time necessary to perform the OSEM-APDI(1%) reconstruction. Chapter 5. Experiments 133 Table 5.10: Summary of image quality measures calculated from the known activity con-centrations in the phantom with those measured for the OSEM-APDI(l%) and OSEM-BB reconstructions using both high statistics and normal statistics projection data. Image Quality Measure and VOIs Used Calculated ("True") High Statistics OSEM- OSEM-APDI(1%) BB Normal Statistics OSEM- OSEM-APDI(1%) BB Contrast between Healthy Heart and Defect 1.000 0.453 0.383 0.458 0.391 Contrast between Healthy Heart and Background 0.718 0.619 0.562 0.634 0.591 NSD in Background N/A 0.23 0.18 0.30 0.25 Figures 5.28 and 5.29 shows reoriented vertical long axis slices through the OSEM-APDI(1%) and OSEM-BB reconstructions for the high statistics and normal statistics data, respectively. This slice was selected to clearly show the inferior wall defect. Each reconstructed image is shown using cubic spline interpolation to reduce pixelization and are displayed using different gray scales, each normalized to the maximum count in that reconstruction. The lower portion of each figure shows horizontal profiles through the original data in the inferior defect region. The position of the profiles is shown as the dashed line on the reconstructed slices. 5.2.4 D I S C U S S I O N The values for the normal statistics projection data, shown as a function of projection angle in figure 5.27, were all close to or less than one. The average values for this data, recorded in table 5.9, was 0.77 with a standard deviation (s.d.) of 0.16. This means that the differences between the experimental projection data and the APDI estimate of the Chapter 5. Experiments 134 0 S E M - A P D I ( 1 % ) O S E M - B B Max Horizontal Profile 70 OSEM-APDI(1%) - - OSEM-BB 90 Fi gure 5.28: Comparison of reoriented vertical long axis slices through the OSEM-APDI(l%), and OSEM-BB reconstructions which used the high statistics projection data. The two images are displayed using different gray scales normalized to the maximum count value in each re-construction. The position of the horizontal profiles shown in the lower figure is shown as the dashed line on the reconstructed slices. Chapter 5. Experiments 135 0 S E M -APDI ( 1 % ) OSEM-BB l u n H a w i t i i i u u n H B n H 4|$W am a s s SSB Max 0 Horizontal Profile 40 45 50 55 60 65 70 75 80 85 90 Figure 5.29: Comparison of reoriented vertical long axis slices through the OSEM-APDI(l%), and OSEM-BB reconstructions which used the normal statistics projection data. The two images are displayed using different gray scales normalized to the maximum count value in each reconstruction. The position of the horizontal profiles shown in the lower figure is shown as the dashed line on the reconstructed slices. Chapter 5. Experiments 136 Table 5.11: Timing results for APDI(1%) and 0SEM-APDI(1%) running on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM. OSEM reconstructions were per-formed using 128x128x38 matrix size, while the APDI calculation was performed on a 64x64x19 sized matrix and expanded using linear interpolation. There were negligible differences in the CPU time required for high and normal statistics projection data. CPU Time OSEM-BB 47 m APDI Scatter Calculation 2.8 h OSEM-APDI 47 m sum of the primary and scattered photon distributions were small relative to the statis-tical uncertainty in these SPECT projection measurements. The values for the high statistics projection data were higher than one for many individual projection angles and was on average 1.75 with a s.d. of 0.25. It should be reiterated that unlike the simulation 2 studies of chapter 4 the measurements for this chapter were necessarily made using both primary and scattered photon data. It is probable that the inclusion of the primary 2 photons would increase the values, since they contain higher detected count values (and hence lower signal to noise ratio) and greater spatial variation than the smoother scattered photon distribution. The errors introduced by using the inaccurate OSEM-BB activity distribution to generate APDI(1%) projections were therefore more significant when the sum of the primary and scattered photon distributions were considered. These differences may be less significant if it was possible to only compare the scattered photon distributions. The average bias values for the normal statistics projection data, that are summarized in table 5.9, show that the mean bias was one standard deviation below unity. The value for the high statistics data was two standard deviations below unity. The lowest bias values occured between -45Â° to 0Â° and from 45Â° to 100Â° in both figures 5.27 (b) and (d). Chapter 5. Experiments 137 These angles corresponded to the camera being located on the left and right hand side of the phantom when the lungs were between the detector and heart insert. These periodic fluctuations in bias values may be connected to the segmented atten-uation map used for these studies. The segmentation routine takes the reconstructed transmission scans and assigns all voxels above a particular threshold the attenuation coefficient for water and all remaining voxels are set to zero. The image is then processed using a smoothing filter. This segmentation and smoothing can reduce the volume of the lungs in the reconstructed attenuation map. An underestimation of the size of the lungs would cause an overestimation of photon attenuation for projection angles where the lungs are between the detector and heart insert. However, for projection angles in which the portion of the phantom between the heart and detector does consist of water (e.g., the anterior view) there would be no overestimation of photon attenuation. Un-fortunately, this segmentation is necessary due to the poor quality of the non-segmented attenuation maps. This hypothesis could be tested if a more accurate attenuation map was available. Another possible cause for these lower bias values is the poor contrast between the heart and background in the original OSEM-BB image which was used to generate the APDI(1%) projections. This reduced contrast was clearly visible in the images and profiles shown in figures 5.28 and 5.29. The technique used to match the experimental data to the APDI(1%) projection estimates does not necessarily preserve counts. In this particular case, it may be difficult to determine an appropriate scaling factor because the ratio of activity in the heart and background is so much further from the true values for the OSEM-BB image. The reader will recall that scaling factors are calculated on a pixel by pixel basis and then a weighted average is used to find the best factor to match the experimental and APDI(1%) distributions. The small size of the phantom's heart, Chapter 5. Experiments 138 relative to the background region of the cylinder, would cause more detector pixels to be influenced by the background region than by the heart. If the detector pixels influenced by the background have sufficiently high count values, a weighted average may favour a lower scaling factor rather than one that would match higher counts in the heart region. This could lead to an overall underestimation of the count values in the scaled APDI(1%) projections. This effect can be seen in figure 5.30, which shows the experimental and APDI(1%) estimate of the projection data at an angle of 50.625Â° for the high statistics projection data. At this angle the bias reached a particularly low value of 0.938. Profiles through the heart region of the projection data are shown on linear and logarithmic scales in figures 5.30 (a) and (b), respectively. The location of the profiles are indicated in the top images as dashed lines. It is clear from these figures that the APDI(1%) estimate was less than the experimental measurements primarily in the heart region. The NSD values for the background VOI, summarized in table 5.10, were consistent with the results of chapter 4. The OSEM-BB reconstruction resulted in lower NSD values than the OSEM-APDI(l%) reconstructions for images reconstructed using both high and normal statistics data. The contrast values in table 5.10 demonstrate that OSEM-APDI(1%) had improved image contrast relative to OSEM-BB. The percent differences between the OSEM-APDI(l%) and OSEM-BB contrast values for the heart and defect was 17% and 16%, while the heart and background contrast was 10% and 7% for the high and normal statistics data, respectively. In all cases, the OSEM-APDI(l%) reconstructions have contrast values closer to those calculated from the known activity concentrations than the OSEM-BB reconstruction. The improvements are slightly less dramatic for the myocardium and defect contrast. The thickness of the myocardium insert (1 cm) was comparable to the pixel length (0.475 cm). The defect VOI generated by the iQuant program was approximately three times Chapter 5. Experiments 139 Horizontal Profile Horizontal Profile Dfitfirtnr Rin Dptertnr Rin Figure 5.30: Comparison of the APDI(1%) estimate of the combined primary and scattered photon distribution with the experimental projections acquired for the high statistics phantom study. The image displayed in the top left is the experimental projection and the top right is the APDI(1%) estimate. Horizontal profiles through the heart region are shown using linear and logarithmic scales in (a) and (b), respectively. The detector was positioned at projection angle 50.625Â°. Chapter 5. Experiments 140 larger than the measured defect insert volume. The partial volume effect and difficulties associated with defining a VOI which is only a few pixels in length may be the reasons the heart and defect contrast was further from the calculated values than the contrast between the heart and background. It should also be noted that the iQuant method assumes that the concentration of activity present in the defect VOI is higher than the surrounding tissue. The method assumes that the user can employ two thresholds. The first threshold is used to create a VOI for the healthy myocardium. A second lower threshold is used to define the entire myocardium including the defect. These two VOIs are then used as a guide to define the defect VOI. If the concentration of activity in the defect is too low, it is not possible to define a threshold that would include the entire myocardium and defect but exclude the surrounding background voxels. In this case, one must define the defect VOI without the use of the second VOI as a guide. The defect VOI is therefore much more difficult to create. The CPU times summarized in table 5.11 indicate that the time necessary to perform the APDI(1%) scatter estimate for the experimental projection data is 2.8 hours, about one third that required for the simulation study. The longer APDI(1%) computation time required for the simulation data was a result of the smaller pixel length and hence larger activity and attenuation maps that were used in these studies. It should also be noted that the initial OSEM-BB and final OSEM-APDI(l%) reconstruction required substantially more time for the experimental study. This was a result of the difference in matrix size in these two studies. In the simulations all reconstructions were performed using a 64x64x24 matrix size while these phantom studies used a 128x128x38 matrix size. Chapter 5. Experiments 141 5 . 2 . 5 CONCLUSIONS 2 The fy values for the normal statistics projection data demonstrate that for a SPECT scan with clinically realistic noise properties the difference between the APDI(1%) es-timate of the scattered and unscattered photon distribution and the experimentally ac-quired data was small relative to the statistical uncertainty in the projection measure-ments. However, for the scan with five times the normal activity concentrations, the difference between the APDI(1%) estimate and the experimental data were significant relative to the projection measurement uncertainty. The bias values for both sets of data were somewhat lower than one, by 1 to 2 standard deviations, indicating that the magni-tude of the APDI(1%) estimate was slightly biased toward lower values. It is hypothesized that these effects may be connected to the use of segmented attenuation maps and the poor contrast of the OSEM-BB image used to generate the APDI(1%) scatter estimate. The NSD values for the background VOI indicated that OSEM-APDI(l%) resulted in higher image noise than OSEM-BB. The contrast values for this study demonstrate that the OSEM-APDI(l%) reconstruction more accurately reflected the true activity con-centrations injected in the phantom than OSEM-BB. The CPU time required for the experimental APDI(1%) scatter estimation was substantially less than that for the simu-lation studies. This was a result of the smaller pixel length and hence larger activity and attenuation maps used in the simulations. However, the experimental reconstructions required more CPU time because the projection matrix size was larger (128x128) in this case. 5.3 PATIENT EXPERIMENTS This section describes the application of the OSEM-APDI algorithm to clinical data. A myocardial perfusion study was selected for a patient that demonstrated an extremely Chapter 5. Experiments 142 elevated liver with very high activity levels, relative to the heart, in the OSEM-BB reconstruction. The data was acquired using the standard cardiac acquisition protocol, described in the previous section, on a Siemens LEHR E-cam camera with the Profile transmission system. The procedure for the generation of APDI scatter estimates, as part of the OSEM-APDI reconstruction, were the same as in the previous section (see figure 5.26). For this study, the OSEM-BB reconstruction was thresholded to include all voxels greater than 0.1% of the maximum count. The reconstruction and scatter estimation step shall be re-ferred to as OSEM-APDI(0.1%) and APDI(0.1%), respectively, throughout the following discussion. This threshold was selected to include all activity voxels that were inside the contours of the attenuation map. Since the injected activity is located inside the patients body, any reconstructed activity outside the attenuation map is presumably an artifact of the OSEM-BB reconstruction. 2 Although the jfy and bias values can be determined as in the previous section, the image quality measurements cannot be applied to patient data. The contrast values cannot be compared to the measured injection because it is impossible to know how much radio-isotope is taken up by a particular organ. The calculation of NSD values is not possible in patient studies because one cannot assume that the concentration of radio-isotope in any particular tissue type or volume is constant. 5.3.1 RESULTS 2 Figure 5.31 shows the the and bias values as a function of projection angle for the patient data. As in the previous section, the sum of the APDI(0.1%) distributions for the primary and scattered photon distributions are compared to the experimental projection 2 measurements. Table 5.12 summarizes the mean jÂŁj and bias values and their standard Chapter 5. Experiments deviations over all projection angles. 143 1.04 1.02 0.98 0.96 0.94 0 50 100 Projection Angle (degrees) 150 9-io r Patient Study | s J l J/ ' I l l ' >h I: '.'V' â€˘ ' a â€˘â€˘ n' i i 1 1 1 ' I â€˘ ' l 11 " l 1 I , 11 l l i i ' * i i i ' i * ' i * l l ' l I I' < 1 ' ' â€˘ 1 " \' / . i 1 1 i ' 1 i ' i â€˘' 'â€˘: !.' "> . \ M ij i ' i i i 1 ; i â€˘ 0 50 100 Projection Angle (degrees) (a) (b) Figure 5.31: Plot of the fj (a) and bias values (b) as a function of projection angle for the patient data. Table 5.12: Comparison of mean fy and the mean bias averaged over all 64 projections and their standard deviations for the patient data. fy- Standard Deviation Bias Standard Deviation Patient Study 1.08 0.19 0.989 0.03 Figure 5.32 shows reoriented short axis slices and profiles through the OSEM-APDI(0.1%) and O S E M - B B reconstructions. The images are displayed in the same manner as in fig-ures 5.29 and 5.28. The amounts of C P U time necessary to perform the initial O S E M - B B reconstruction, the APDI(0.1%) scatter estimation step, and the OSEM-APDI(0.1%) reconstruction, are summarized in table 5.13. Chapter 5. Experiments 144 Figure 5.32: Comparison of reoriented short axis slices through the OSEM-APDI(0.1%) and OSEM-BB reconstructions for the patient data. The plots in the lower portion of the figure show horizontal and vertical profiles through the slices, indicating regions of the plot corresponding to the left ventricle (L.V.), right ventricle (R.V.) and gall bladder (G.B.). Chapter 5. Experiments 145 Table 5.13: Timing results for the patient data for each of the OSEM-BB, APDI(0.1%) and OSEM-APDI(0.1%) calculations running on a computer with a 1.7 GHz processor and 1024 Mb of RDRAM. OSEM reconstructions were performed using 128x128x38 ma-trix size, while the APDI(0.1%) calculation was performed on a 64x64x19 sized matrix and expanded using linear interpolation. CPU Time OSEM-BB 43 m APDI Scatter Calculation 3.3 h OSEM-APDI 43 m 5.3.2 D I S C U S S I O N 2 For most projection angles, the and bias values shown in figure 5.31 fluctuated around 2 one. The average and bias values, summarized in table 5.12, were both within one standard deviation of unity. As in the phantom experiments, the lowest bias values occurred in the range of projection angles between -45Â° to 0Â° and between 45Â° to 100Â°. This may have been due to the reduced volume of the lungs in the attenuation maps that was described previously. The short axis slices and profiles shown in figure 5.32 shows that the OSEM-APDI(0.1%) reconstruction had improved contrast and better separation of the heart from nearby organs (i.e., the liver and gall bladder) relative to OSEM-BB. Figure 5.32 also shows increased uptake, probably due to scatter from the liver, in the inferior wall region of the OSEM-BB image which was reduced in the OSEM-APDI(0.1%) image. Table 5.13 indicates that the CPU time required to perform APDI(0.1%) scatter estimation step is about 3.3 hours. This was slightly longer than the phantom study due to the larger physical size of the patient relative to the phantom, and larger activity and attenuation map used in this study. Chapter 5. Experiments 146 5.3.3 C O N C L U S I O N S 2 The fj values for this section indicate that the average difference between the sum of the APDI(0.1%) primary and scattered photon distributions and the clinical data was small relative to the statistical noise in the data. Although the average bias for this data set was within a standard deviation of unity, there does seem to be a systematic underestimation associated with the APDI calculation for specific projection angles which was observed in both phantom and patient data sets. It is hypothesized that this error is due to the poor quality of the experimentally determined attenuation maps. The increased uptake in the inferior wall region observed in the OSEM-BB image, relative to OSEM-APDI(0.1%), could potentially obscure the presence of an inferior wall defect and would certainly result in an overestimation of the level of perfusion in this region. Although the CPU time required for the APDI(0.1%) scatter estimation was slightly longer than the phantom study, it was still substantially shorter than that required for the simulations. ' C H A P T E R 6 C O N C L U S I O N S A N D F U T U R E W O R K In this work APDI, a scatter estimation technique based on the APD method, and OSEM-APDI, its implementation as a scatter correction in an iterative reconstruction algorithm, were presented. The correction and reconstruction were applied to computer simulated projection data, experimental data obtained from physical phantoms, and patient data. The following chapter provides a brief overview of the work performed thus far and some general comments about the applicability of the method proposed here. The chapter concludes with a description of future work which could improve the method's performance and accuracy. 6.1 OVERVIEW It has been demonstrated that the approximate APDI scatter estimation method gen-erates scatter projections that are qualitatively and quantitatively very similar to those generated using the full APD method. Using clinically realistic digital activity and at-tenuation maps, it was estimated that the average differences between APD and APDI would become significant (i.e., they would be comparable to the statistical uncertainty in SPECT projection measurements) for acquisition times between 20 to 30 times longer 147 Chapter 6. Conclusions and Future Work 148 than a typical clinical scan for the liver and heart activity distributions, respectively. These errors were slightly higher using the body activity distribution, and would become significant for scans that are about 6 times longer than a clinical scan. The optimal number of times the APDI scatter estimate should be updated in the OSEM-APDI reconstruction was then investigated. Monte-Carlo and noiseless simulated projection data was generated for a digital phantom with a perfusion defect in the inferior wall of the myocardium. It was concluded that there was little need to perform more than one scatter estimation step in OSEM-APDI. The APDI scatter projections generated us-ing the OSEM-BB reconstruction were not significantly different from those generated using subsequent OSEM-APDI iterations. The image quality measures considered also did not change substantially when additional scatter estimation steps were performed. Therefore, the OSEM-APDI reconstruction (which included only one APDI scatter esti-mation step) was performed using simulated SPECT data for a series of phantoms with perfusion defects in various positions in the myocardium. In all cases, the OSEM-APDI images had substantially better contrast than the OSEM-BB images and were extremely close to OSEM-ISR, the ideal scatter-free reconstruction, in terms of image uniformity and contrast. The OSEM-APDI reconstruction was then applied to two sets of experimental data each with different noise characteristics. Both employed a physical phantom to model the activity and attenuation distribution for a patient with an inferior wall defect in a myocardial perfusion study. The OSEM-APDI reconstructions demonstrated contrast values which were closer to those calculated from the known activity concentrations than the OSEM-BB reconstructions. Finally, the OSEM-APDI reconstruction was applied to patient data. It was observed that the OSEM-BB image demonstrated increased uptake in the inferior wall, probably due to scatter from the liver, relative to OSEM-APDI. Chapter 6. Conclusions and Future Work 149 6.2 CURRENT APPLICABILITY Although the average differences between APD and APDI were small, it should be stressed that this was the case for a specific set of acquisition parameters such as the length of scan, the SPECT detector model, a LEUHR collimator geometry and the par-ticular activity and attenuation distributions used. If any of these parameters were to change these results may not be applicable. There also may be situations where the dif-ferences between APD and APDI observed in individual scatter projection pixels may be significant. The APDI method is sensitive to large variations in the attenuation distribu-tion and is therefore less accurate for activity voxels at the boundaries between different attenuating media. The APDI method also tends to produce scatter projections that are more "blurred" than the APD method. One can conceive of situations where these problems could become significant. A specific example would be lung tumor oncology studies, in which a small hot source could be located near the boundary of the lungs and soft tissue. In situations such as this, the author recommends the use of the full APD method which can be performed relatively quickly for small sources and would be much more accurate in this case. The CPU time required for the APDI calculation depends strongly on the number of object voxels containing activity and attenuating material. In the simulation studies, a fairly severe threshold was applied to the OSEM-BB activity distribution prior to its use in the APDI scatter estimation step. The application of this threshold was necessary to reduce the CPU time to a reasonable amount for the particular phantoms used in this study. However, the author does not recommend the use of extreme thresholds unless absolutely necessary. In general, the APDI threshold should be selected to include all activity voxels in the OSEM-BB reconstruction that are within the contours of the attenuation map, based on visual inspection of the two images. Chapter 6. Conclusions and Future Work 150 In the noiseless simulation studies the OSEM-APDI images had contrast values that consistently fell somewhere between that of OSEM-ISR and OSEM-APDt. It was con-cluded that an overestimation had occurred in the APDI scatter projection for pixels which are significant to myocardial contrast. It was hypothesized that this effect was connected to the poorer contrast in the OSEM-BB activity distribution, which can lead to an over-correction for scatter when cold spots are present in the object being imaged. In this particular case, these inaccuracies did not result in contrast values greater than the ideal correction, OSEM-ISR. However, it is obvious that poor contrast in the OSEM-BB image can lead to less accurate APDI scatter projections. This effect manifested itself in a different way in the phantom experiments. In this case, the APDI scatter projections were biased toward lower values for pixels influenced by the heart. In some cases, it may therefore be preferable to perform more than one scatter it-eration in the OSEM-APDI reconstruction, depending on the quality of the OSEM-BB image used in the first APDI scatter estimation step. Even when the exact truth is not known, one can predict when an OSEM-BB image will be particularly poor. For example, from visual inspection of the OSEM-BB reconstruction, one could observe if a patient has an extremely hot liver in close proximity to the heart. In this situation, the OSEM-BB image would have reduced contrast in the heart region, which could result in less accurate scatter projections in the first scatter iteration of OSEM-APDI. This effect would become more severe as the level of activity in the liver increases relative to that in the heart. Under these circumstances, the author recommends that additional scatter iterations should be performed in OSEM-APDI. Chapter 6. Conclusions and Future Work 151 6.3 FUTURE WORK One of the major obstacles to the development of quantitative SPECT imaging is the computation time required for accurate scatter corrections. A comparison of the amounts of CPU time that were necessary to perform each APDI corrected reconstruction for all of the data sets considered in this work are summarized in Table 6.1. The APDI scatter estimation method can be performed in about 3 hours for the experimental data but required about 9 hours for the simulation studies. In general, a reconstruction technique is considered to be clinically viable if it can be performed in approximately the same time it takes to scan a patient 30 minutes). This is a practical limitation which prevents backlogs in the clinics. If this limitation is met, while a new patient is being scanned, the previous patient data can be reconstructed, essentially allowing data to be acquired continuously and thereby maximizing the number of patients that can be scanned in a work day. Table 6.1: Timing results (1.7 GHz Processor, 1024 Mb RDRAM) for each step of the scatter correction and reconstruction for the simulation, phantom and patient data sets. CPU Time Simulation Phantom Patient OSEM-BB 280s 47m 43m APDI Calculation 9.2h 2.8h 3.3h OSEM-APDI 280s 47m 43m Threshold in APDI 15% 1% 0.1% No. Voxels in APDI Â«26000 Â«12000 ^15000 Matrix Size in APDI 64x64x24 64x64x19 64x64x19 Matrix Size in OSEM 64x64x24 128x128x38 128x128x38 In order reduce the computation time required for the APDI scatter calculation, a number of changes must be made to the algorithm. A large amount of CPU time in both the APD and APDI calculation (Â«65% and Â«30%, respectively) is spent determining Chapter 6. Conclusions and Future Work 152 the attenuation length: f*~ p(x)dx, from the source position: S j , to the final scattering position: t^. If some of these attenuation lengths could be reused, a great deal of time could be saved. Since a separate APDI calculation is performed for each projection angle, many attenuation lengths through the object are calculated repeatedly from one projection angle to the next. One obvious solution is to determine all projection angles simultaneously. However, there are some practical problems with this approach. The first problem is the size of the APD look-up table. If this method were imple-mented naively, the look-up table size would increase by the number of projection angles. Currently the size of the APD look-up tables are approximately 25 Mb. An APD look-up table which included separate elements for 64 different projection angles, would be 1.6 Gb. This is larger than the random access memory (RAM) of the computer used for most of the timing results in this thesis, which had 1 Gb RDRAM (most of the computers used in clinics have even less RAM). Therefore, this computer would not be able allocate the entire look-up table in RAM, which would increase the computation time substantially. The size of the look-up tables that include all projection angles simultaneously can be decreased by restricting the number of projections based on the energy of the scattered photon. A schematic diagram to illustrate how this could be done using a simplified 2D model is shown in figure 6.1. If we ignore the energy resolution of the detector, 126 keV is the minimum energy of a scattered 9 9 m T c photon which would be detected in a symmetric 20% energy window around 140 keV. This corresponds to a maximum first order Compton scattering angle of 53Â°. If we assume a perfect collimator, only those scattered photons with directions perpendicular to the the detector surface can be detected. If data is acquired over a 360Â° rotation and the problem is restricted to 2D, there would always be the same number of projection angles for which look-up table calculations would need to be performed (23x650f). Clearly the exact position of this subset of projection angles Chapter 6. Conclusions and Future Work 153 would depend on the incident photon direction (see figure 6.1). However, the size of this subset would increase if one includes the imperfect energy response of the detector and the collimator spatial resolution. The problem also becomes more complex if one includes second order scatter and extends the problem to 3D. Figure 6.1: Schematic diagram of the subset of projection angles for which APD look-up table elements would need to be calculated for two different incident photon directions. The problem here is restricted to 2D, the detector energy response and collimator spatial resolution are assumed to be perfect, and data is collected in a symmetric 20% energy window around 140 The second problem is connected to the flexibility of the calculation. In the method described here, one can easily account for a variable angular interval between projections and non-circular detector orbit by rotating and shifting the activity and attenuation map before they are input in the APDI scatter calculation. If look-up table elements were calculated for all projection angles, the angular interval and radius of rotation would have to be fixed. Typically the NCO orbit changes based on the object being imaged. This would cause difficulties because the look-up table calculation involves F(ÂŁ), the geometrical collimator response, which is depth dependent. Therefore, there is no simple way to account for a variable radius of rotation in an APD look-up table calculation with multiple projection angles. One possible solution would be to adopt an approach similar Detector "Path""' Detector ""Path""--.. keV. Chapter 6. Conclusions and Future Work 154 to that used in [88], in which all scattered photons are assumed to travel along paths perpendicular to the detector surface. The distance from the last scattering site to the detector is used to perform a convolution with a depth dependent Gaussian function, to account for the collimator response, in the final stage of the calculation. In order to further verify the accuracy of OSEM-APDI, it may be necessary to test the algorithm using more diverse and therefore more realistic simulation data. In the work performed in this thesis, all the simulation studies essentially used the same phantom with only one free variable; the location of the perfusion defect. In the MC studies the same statistical realization of the data was used for each of the components of the body and liver activity distributions. To validate more thoroughly this scatter correction it may be useful to perform separate MC simulations for a population of digital phantoms, as in [77], each with different body types, organ sizes and orientations. In this way, the OSEM-APDI method could be tested using data that more accurately reflects the type of variation observed in a clinical patient population. Another area that may require further validation is the use of the OSEM-BB activity distribution to generate the APDI scatter projections. For all the cases considered in this thesis, the results were reasonable. However, for reasons discussed in the previous section, this may not always be the case. A possible solution would be a more incremental updating of the scatter estimation in the reconstruction routine. The programs used to generate the APDI scatter estimate and the OSEM-APDI reconstruction were completely separate. If the two programs could be combined efficiently then updating the APDI scatter estimate could be performed after each OSEM iteration. This would avoid the direct use of the OSEM-BB reconstruction in APDI, which can contain artifacts precisely because scatter has not been accurately corrected for. A P P E N D I X A C O M P A R I S O N O F A P D W I T H S I M S E T In section 4.3 the spatial distribution of scattered photons in the SimSET-simulated projection data was compared to the APDI scatter estimates generated in the OSEM-APDI(15%) and OSEM-APDI(0%)P reconstructions. In this comparison, it was observed 2 that all fy values were greater than one and the bias values were substantially lower than one. If the APDI data accurately reflects the underlying spatial distribution of a scattered photons and the SimSET data has the correct statistical variance, the fy and bias values should both be approximately equal to one. The purpose of this appendix is to demonstrate that systematic differences between the SimSET and APD (and hence APDI) data caused this discrepancy, rather than any problem with the OSEM-APDI reconstruction method. A comparison was made between the SimSET and APD data (i.e., without the APDI interpolation approximation) using the true MCAT activity distribution for the phan-tom with an inferior wall defect. The APD and SimSET data was matched using b, the weighted average scaling factor defined by equation 3.34. To investigate if this dis-crepancy was caused by this scaling factor, another analysis was performed using APD data rescaled to have the same total number of counts (i.e., the sum of both scattered 155 Appendix A. Comparison of APD with SimSET 156 and primary photons) as the SimSET data. A comparison was also made between APD and SimSET data which had been separated into three components corresponding to the contributions from primary photons, first order Compton scattered photons, and second order Compton scattered photons. A . l RESULTS A comparison of the projections calculated using the full APD method and the SimSET simulated data are shown in figures A.2 and A. l for detector angles 180Â° (posterior view) and 360Â° (anterior view), respectively. This data represents the total photon distribu-tions including scattered and unscattered photons. The lower parts of these figures show horizontal profiles drawn through the pixel with maximum count value on linear and log-arithmic scales. These graphs show the total from all photons and the contribution from each component (primary photons, first and second order Compton scattered photons). The SimSET data is shown using various markers, while the APD data is shown using different line styles. To verify that the differences observed between the APD and SimSET scattered pho-ton projection data was not a result of using the weighted average scaling factor, b, an analysis was also performed with APD data which had been rescaled using the total sum of the counts. This scaling factor can be written as: T y S i m S E T &TSC = ~jr-^ , (A.l) ÂŁ y A P D ( A T r " e , n t ) t = i where the summation from t = 1 to T is over all projection pixels, Y t S i m S E T is the SimSET count value (including both scattered and unscattered photon contributions) for detector Appendix A. Comparison of APD with SimSET 157 Max SimSET fess Max Horizontal Profile Horizontal Profile 400 40 Detector Bin (a) 20 40 60 Detector Bin (b) Figure A . l : Comparison of the full APD and SimSET projection data. The projection images are shown for the sum of scattered and unscattered photons, for the MCAT phantom with an inferior wall defect and the detector at the anterior view (360Â°). Horizontal profiles through the maximum count value are shown using linear and logarithmic scales in (a) and (b), respectively. Data in the plots are shown separately for all photons, unscattered photons, first order and second order Compton scattered photons. Appendix A. Comparison of APD with SimSET 158 Max 400 Max Horizontal Profile Horizontal Profile o SimSET Total â€˘ SimSET Primary A SimSET First Order â€˘ SimSET Second Order APD Total - - APD Primary - â€˘ APD First Order APD Second Order 20 40 60 Detector Bin 20 40 60 Detector Bin 80 (a) (b) Figure A.2: Same as figure 4.10, except that the projections shown are for detector at the posterior view (180Â°). Appendix A. Comparison of APD with SimSET 159 pixel *, and Y A P D ( \ T r v e , nt) is the sum of the full APD primary and scattered photon distribution generated using the true MCAT activity distribution. Table A . l compares the total sum of the counts and the counts from each of the individual scatter components using both scaling factors. Table A. l : Comparison of the total sum of the counts for the SimSET and APD projection data. Also included is the contribution from primary photons, first and second order Compton Scattered photons. The APD data using both scaling factors, the weighted average (WA) and the total sum of the counts (TSC), are included. Dataset Total Number Primary Compton First Compton Second of Counts Order Order APD (WA) 8380561 6514647 1702302 163612 APD (TSC) 8022800 6236541 1629631 156628 SimSET 8022800 5930365 1782068 310367 The SimSET data and the APD scattered photon distributions (using both scaling factors) are compared in figure A.3 using the fy and bias figures of merit as a function of projection angle. The mean and standard deviation of the fy and bias values over all projection angles are summarized in table A.2. Table A.2: Comparison of the SimSET and APD scattered photon data (with both scaling factors) using the mean fy and the mean bias over all 64 projections and their standard deviations Scaling Factor fy Standard Deviation Bias Standard Deviation Weighted Average Sum of the Counts 1.42 0.06 1.69 0.08 0.89 0.01 0.85 0.01 Appendix A. Comparison of APD with SimSET 160 â€” Weighted Average - - Sum of the Counts 1.9h I Projection Angle (degrees) I â€” Weighted Average - - - Sum of the Counts A â€˘ M i * A " â€˘ 1 B 1 Â» x1 \ i * . â€˘ ' l . * ' f \ l Â» -Â» â€˘ " l ' Â« 1 V l l , ! Â» " , ;, % ' V * â€˘ > 1 x .. V, . n H , l / Â» ' 1 1 1 ' Â» u 1 1 ' .'â€˘ t 50 100 150 200 250 300 350 400 450 Projection Angle (degrees) (a) (b) 2 Figure A.3: Plot of (a) and bias (b) as a function of projection angle to compare the SimSET and APD scattered photon data. The solid and dashed lines are for the APD data using scaling factors determined using the weighted average and total sum of the counts, respectively. A . 2 DISCUSSION The total APD and SimSET photon distributions (sum of the scattered and unscattered photons) shown in figures A.2 and A . l appears to agree fairly well. However, there are discrepancies in the individual components from primary photons, first and second order Compton scattered photons. In particular, there is a large systematic difference between SimSET and APD in the relative magnitude of the second order Compton scatter distribution. This is also apparent in table A . l in terms of the ratio of each photon component to the total sum of the counts. For SimSET 74% of the total counts come from primary photons, 22% from first order Compton scatter, and 4% from second order. For APD the percentages are 78% primary, 20% first order, and 2% second order. When the scaling factor &TSC was used the APD total number of counts from primary photons became closer to SimSET result. However, the APD total counts from Compton scatter for both orders were further from the SimSET values using &TSC (see table A.l). The use Appendix A. Comparison of APD with SimSET 161 of &TSC a l s o resulted in poorer fy and bias values for the scattered photon distributions (see figure A.3 and table A.2). This seems to indicate that the weighted average scaling factor, b, is preferable when the photon distributions are not exactly matched in terms of the ratio of each photon component. This indicates that systematic differences exist between the APD and SimSET data. The ratio of scattered photons to primary photons was different in SimSET than in APD. This makes it difficult to match the SimSET-simulated projections to the sum of the APD scattered and primary photon distributions. The plots of fy and bias in this section (figure A.3) closely resemble those in section 4.3 (figures 4.15 (c) and (d)). Therefore, these problems also manifested themselves in the APDI data generated in the OSEM-APDI reconstruction technique which used the SimSET data. It is not clear where these systematic differences came from. One possible explanation considered was that assumptions made in the derivation of the original APD method introduced some large inaccuracies. In particular, an assumption in the second order lookup table calculation could introduce a bias toward lower values. It is assumed that the attenuation between the point of emission and the first scattering site is the same as that for water at the unscattered photon energy. If an attenuation map contained a large number of voxels with lower attenuation than water (such as lung tissue), the APD calculation could produce scatter projections biased toward lower values (i.e., APD over-estimates the photon attenuation). This hypothesis was tested by comparing the full APD and SimSET projections calculated using an attenuation map that consisted entirely of water. If the hypothesis was correct no bias would be introduced by the second order approximations in this case. However, the bias was still observed. A great deal of work was done to attempt to sort out this problem. What was particularly puzzling was that the original APD method agreed quite well with an earlier version of Appendix A. Comparison of APD with SimSET 162 the SimSET program [102, 3, 2]. Another possible explanation considered was that the variance reduction employed in this work was somehow used incorrectly. If the weight (see section 2.5) of the second order Compton scattered photons was overestimated, it could explain these differences. If time had permitted it, this hypothesis could have been tested by comparing these results with a different MC program. A . 3 CONCLUSIONS The only conclusions that could be made was that there are systematic differences be-tween APD and SimSET data, with respect to the contributions from unscattered, first 2 and second order Compton scattered photons. Lower fy and bias values for the scattered photon distributions were achieved using the weighted average scaling factor as opposed to the total sum of the counts. The cause of the systematic differences between APD and SimSET is still unknown. B I B L I O G R A P H Y [1] C. Floyd, R. Jaszczak, C. Harris, and R. Coleman. Energy and spatial distribution of multiple order Compton scatter in SPECT: a Monte Carlo investigation. Phys. Med. Biol., 29:1217-1230, 1984. [2] R.G. Wells, A. Celler, and R. Harrop. Analytical calculation of photon distributions in SPECT projections. IEEE Trans. Nucl. Sci, 45(6):3202-3214, 1998. [3] R.G. Wells. Analytical calculation of photon distributions in SPECT projections. PhD thesis, University of British Columbia, 1997. [4] L. Freeman, editor. Freeman and Johnson's Clinical Radionuclide Imaging. Grune and Stratton, Inc., Orlando Florida, third edition, 1984. [5] E.L. Kramer and J.J. Sanger. Clinical SPECT Imaging. Raven Press, 1995. [6] P.J Early and D.B Sodee. Principles and Practice of Nuclear Medicine. Mosby-Year Book, Inc., St. Louis, 1995. [7] G.N. Ramachandran and A.V. Lakshminarayanan. Three-dimensional reconstruc-tion from radiographs and electron micrographs: application of convolution instead of fourier transforms. Proc. Nat. Acad. Sci., 68(9)-.2236-2240, 1971. [8] C E . Metz, F.B. Atkins, and R.N. Beck. The geometric transfer function component for scintillation camera collimators with straight parallel holes. Phys. Med. Biol., 25(6):1059-1070, 1980. [9] B. Hutton, H. Hudson, and F. Beekman. A clinical perspective of accelerated statistical reconstruction. Eur. J. Nucl. Med., 24:797-808, 1997. [10] L.A. Shepp and Y. Vardi. Maximum liklihood reconstruction for emission tomog-raphy. IEEE Trans. Med. Imag., MI-L113-122, 1982. [11] H.M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag., 13:601-609, 1994. [12] G.L. Zeng and G.T. Gullberg. Frequency domain implementation of the three-dimensional geometric point response correction in SPECT imaging. IEEE Trans. Nucl. Sci., 39:1444-1453, 1992. [13] G. Zeng and G. Gullberg. Valid back projection matrices which are not the trans-pose of the projection matrix (abstract). J. Nucl. Med., 37:206p, 1996. 163 Bibliography 164 [14] C Kamphuis, F. Beekman, M. Viergeve, and P. Rijk. Accelerated fully 3D SPECT reconstructions using dual matrix ordered subsets (abstract). J. Nucl. Med., 37:62p, 1996. [15] A. Welch and G.T. Gullberg. Implementaion of a model-based nonuniform scatter correction scheme for SPECT. IEEE Trans. Med. Imag., 16:717-726, 1997. [16] R.J. Jaszczak, K.L. Greer, C E . Floyd, C C . Harris, and R.E. Coleman. Improved SPECT quantification using compensation for scattered photons. J. Nucl. Med., 25(8):893-900, 1984. [17] M. A. King, W. Xia, D. J. deVries, T-S. Pan, B. J. Villegas, S. Dahlberg, B. M. W. Tsui, M. H. Ljungberg, and H. T. Morgan. A monte carlo investigation of artifacts caused by liver uptake in single-photon emission computed tomography perfusion imaging with technetium 99m-labeled agents. J. Nucl. Cardiol, 3:18-29, 1996. [18] Y.S.Gur, T.H. Farncombe, P.H. Pretorius, H . C Gilford, M.V. Narayanan, E . C Frey, D. Gagnon, and M.A. King. Comparison of scatter compensation strategies for myocardial perfusion imaging using tc-99m labeled sestamibi. IEEE Trans. Nucl. Sci, 49(5):2309-2314, October 2002. [19] M. O'Connor, C. Caiati, T. Christian, and R. Gibbons. Effects of scatter correction on the measurement of infarct size from SPECT cardiac phantom studies. J. Nucl. Med., 36:2080-2086, 1995. [20] I. Buvat, H. Benali, A. Todd-Pokropek, and R. Di Paola. Scatter correction in scintigraphy: the state of the art. Eur. J. Nucl. Med., 21:675-694, 1994. [21] I. Buvat, M. Rodriguez-Villafuerte, A. Todd-Pokropek, H. Benali, and R. Di Paola. Comparative assessment of nine scatter correction methods based on spectral anal-ysis using Monte Carlo simulations. J. Nucl. Med., 36:1476-1488, 1995. [22] M. Ljungberg. Development and evaluation of attenuation and scatter correction techniques for SPECT using the Monte Carlo method. PhD thesis, University of Lund, Sweden, 1990. [23] M.A. King, B.C. Penney, and S.J. Glick. An image-dependent metz filter for nuclear medicine images. J. Nucl. Med., 29:1980-1989, 1988. [24] B.C. Penney, M.A. King, and S.J. Glick. Restoration of combined conjugate images in spect: comparison of a new wiener filter and the image-dependent metz filter. IEEE Trans. Nucl. Sci., 37:707-712, 1990. Bibliography 165 C. C. Harris, K. L. Greer, R. J. Jaszczak, C. E. Floyd, E. C. Fearnow, and R. E. Coleman. Tc-99m attenuation coefficients in water-filled phantoms determined with gamma cameras. Med. Phys., 11:681-685, 1984. S.A. Larsson. Gamma Camera Emission Tomography: development and properties of a multi-sectional emission computed tomography system. PhD thesis, University of Stockholm, Sweden, 1980. R.J. Jaszczak, R.E. Coleman, and F.R. Whitehead. Physical factors affecting quan-titative measurements using camera-based single photon emission computed tomog-raphy (SPECT). IEEE Trans. Nucl. Sci, NS-28:69-79, 1981. D. R. Gilland, R.J. Jaszczak, K.L. Greer, and R.E. Coleman. Quantitative SPECT reconstruction of iodine-123 data. J. Nucl. Med., 32:527-533, 1991. J.A. Siegel, R.K. Wu, and A.H. Maurer. The buildup factor: Effect of scatter on absolute volume determination. J. Nucl. Med., 26:390-394, 1985. R.K. Wu and J.A. Siegel. Absolute quantification of radioactivity using the buildup factor. Med. Phys., 11:189-192, 1984. M. Ljungberg and S-E. Strand. Scatter and attenuation correction in SPECT using density maps and Monte Carlo simulated scatter functions. J. Nucl. Med., 31:1579-1567, 1990. M. Ljungberg and S-E. Strand. Attenuation corrections in SPECT based on trans-mission studies and Monte Carlo simulations of buildup functions. J. Nucl. Med., 31:493-500, 1990. W.L. Rogers, N.H. Clinthorne, J. Stamos, K.F. Koral, R. Mayans, G.F. Knoll, J. Juni, J.W. Keyes, and B.A. Harkness. Performance evaluation of sprint, a single photon ring tomography for brain imaging. J. Nucl. Med., 25:1013-1018, 1984. R. Jaszczak. Scatter compensation techniques for SPECT. IEEE Trans. Nucl. Sci, NS-32:786-793, 1985. G. Hademenos, M. Ljungberg, M. King, and S. Glick. A Monte Carlo investigation of the dual photopeak window scatter correction method. IEEE Trans. Nucl. Sci., 40:179-185, 1993. M. Ljungberg, P. Msaki, and S-E. Strand. Comparison of dual-window and con-volution scatter correction techniques using the Monte Carlo method. Phys. Med. Biol, 35:1099-1110, 1990. Bibliography 166 K.F. Koral, F.M. Swailem, S. Buchbinder, N.H. Clinthorne, W.L. Rogers, and B.M.W. Tsui. SPECT dual-energy-window compton correction: scatter multiplier required for quantitation. J. Nucl. Med., 31(l):90-98, 1990. P.H Pretorius, A.J. van Rensburg, A. van Aswegen, M.G. Lotter, D.E. Serfontein, and CP. Herbst. The channel ratio method of scatter correction for radionuclide image quantitation. J. Nucl. Med., 34(2):330-335, 1993. K.W. Logan and W.D. McFarland. Single photon scatter compensation by photo-peak energy distribution analysis. IEEE Trans. Med. Imaging, 11:161-164, 1992. M.A. King, G. Hademeos, and S.J. Glick. A dual photopeak window method for scatter correction. J. Nucl. Med., 33:606-612, 1992. K. Ogawa, Y. Harata, T. Ichihara, A. Kubo, and S. Hashimoto. A practical method for position-dependent Compton-scatter correction in single photon emission CT. IEEE Trans. Med. Imag., 10:408-412, 1991. K. Ogawa and N. Nishizaki. Accurate scatter compensation using neural networks in radionuclide imaging. IEEE Trans. Nucl. Sci., 40(4): 1020-1025, 1993. R.N. Beck, L.T. Zimmer, D.B. Charleston, and P.B. Hoffer. Aspects of imaging and counting in nuclear medicine using scintillation and semiconductor detectors. IEEE Trans. Nucl. Sci, 19:173-178, 1972. J.R. Halama, R.E. Henkin RE, and L.E. Friend. Gamma camera radionuclide images: Improved contrast with energy-weighted acquisition. Radiology, 169:533-538, 1988. L.V. East, R.L. Phillips, and A.R. Strong. A fresh approach to Nal scintillation detector spectrum analysis. Nucl. Instrum. Methods, 193:147-155, 1982. K. Koral, X. Wang, L. Rogers, N. Clinthorne, and X. Wang. SPECT Compton-scattering correction by analysis of energy spectra. J. Nucl. Med., 29:195-202, 1988. X. Wang and K. Koral. A regularized deconvolution-fitting method for Compton-scatter correction in SPECT. IEEE Trans. Med. Imag., ll(3):351-360, 1992. F. Cavailloles, D. Morvan, F. Boudet, J.P. Bazin, and R. Di Paola. Factor anal-ysis of dynamic structures as an aid for vesicoureteral reflux diagnosis. Contrib. Nephrol, 56:238-242, 1987. Bibliography 167 [49] D. Gagnon, A. Todd-Pokropek, A. Arsenault, and G. Dupras. Introduction to holospectral imaging in nuclear medicine for scatter subtraction. IEEE Trans. Med. Imag., 8(3):245-250, 1989. [50] D. Gagnon, N. Pouliot, L. Laperriere, A.Arsenault J. Gregoire, and G. Dupras. Post acquisition linearity correction for holospectral imaging in nuclear medicine. IEEE Trans. Nucl. Sci, 37(2):621-626, 1990. [51] W.H. Press, B.R. Flannert, S.A. Teukolsky, and W.T. Vetterling. Numerical Recipes. Press Syndicate of the University of Cambridge, New York, 1986. [52] D.C. Barber. The use of principal components in the quantitative analysis of gamma camera dynamic studies. Phys. Med. Biol, 25:283-292, 1980. [53] R. Di Paola, J.P. Bazin, F. Aubry, A. Aurengo, F. Cavailloles, J.Y. Herry, and E. Khan. Handling of dynamic sequences in nuclear medicine. IEEE Trans. Nucl. Sci, 29:1310-1321, 1982. [54] P. Hannequin, J.C. Liehn, and J. Valeyre. Correction de la diffusion compton par analyse factorielle des structures dynamiques. J. Med. Nucl. Biophys., 12:460, 1989. [55] I. Buvat. Correction de la diffusion Compton en imagerie scintigraphique par anal-yse factorielle sous contraintes des structures spectrales. Master's thesis, Institut Gustave Roussy, 1990. [56] I. Buvat. Correction de la diffusion en imagerie scintigraphique. PhD thesis, Uni-versite de Paris XI, 1992. [57] B. Axelsson, P. Msaki, and A. Israelsson. Subtraction of Compton-scattered pho-tons in single-photon emission computerized tomography. J. Nucl. Med., 25:490-494, 1984. [58] C E . Floyd, R.J. Jaszczak, K.L. Greer, and R.E. Coleman. Deconvolution of Comp-ton scatter in SPECT. J. Nucl. Med., 26(4):403-408, 1985. [59] P. Msaki, B. Axelsson, C M . Dahl, and S.A. Larsson. Generalized scatter correc-tion method in SPECT using point scatter distribution functions. J. Nucl. Med., 28(12):1861-1869, 1987. [60] T. Mukai, J.M. Links, K.H. Douglass, and H.N. Wagner Jr. Scatter correction in spect using non-uniform attenuation data. Phys. Med. Biol., 33:1129-1140, 1988. Bibliography 168 [61] S.R. Meikle, B.F. Hutton, D.L. Bailey, R.R. Fulton, and K. Schindhelm. SPECT scatter correction in non-homogeneous media. In A.C.F. Colchester and D.J. Hawkes, editors, Information Processing in Medical Imaging, XII IPMI Interna-tional Conference, pages 34-44, Kent, 1991. Springier-Verlag. [62] P. Msaki. Position-dependent scatter response functions: will they make a differ-ence in SPECT conducted with homogeneous cylindrical phantoms? Phys. Med. Biol, 39:2319-2329, 1994. [63] E. Frey and B. Tsui. A practical method for incorporating scatter in a projector-backprojector for accurate scatter compensation in SPECT. IEEE Trans. Nucl. Sci., 40:1107-1116, 1993. [64] R.J. Fleming. A technique for using ct images in attenuation correction and quan-tification in SPECT. Nucl. Med. Comm., 10:83-97, 1989. [65] E. Frey and B. Tsui. Parameterization of the scatter response function in SPECT imaging using Monte Carlo simulation. IEEE Trans. Nucl. Sci., 37:1308-1315, 1990. [66] E.C. Frey and B.M.W. Tsui. A practical projector-backprojector modeling attenu-ation, detector response, and scatter for accurate scatter compensation in SPECT. In IEEE Conference record of the 1991 Nuclear Science Symposium and Medical Imaging Conference, volume 3, pages 1777-1781, 1991. [67] E. Frey, Z.W. Ju, and B.M.W. Tsui. A fast projector-backprojector pair modeling the asymmetric, spatially varying scatter response function for scatter compensa-tion in spect imaging. IEEE Trans. Nucl. Sci., 40:1192-1197, 1993. [68] E. Frey, B. Tsui, and M. Ljungberg. A comparison of scatter compensation meth-ods in SPECT: subtraction-based techniques versus iterative reconstruction with accurate modeling of the scatter response. In IEEE Conference record of the 1992 Nuclear Science Symposium and Medical Imaging Conference, pages 1035-1037, November 1993. [69] F. Beekman, E. Eijkman, M. Viergever, G. Borm, and E. Slijpen. Object shape dependent PSF model for SPECT imaging. IEEE Trans. Nucl. Sci., 40:31-39, 1993. [70] F.J. Beekman, E.C. Frey, C. Kamphuis, B.M.W. Tsui, and M.A. Viergever. A new phantom for fast determination of the scatter response of a gamma camera. IEEE Trans. Nucl. Sci., 41:1481-1488, 1994. Bibliography 169 [71] C E . Frey, Z.W. Ju, and B.M.W. Tsui. Modeling the scatter response function in inhomogeneous scattering media for SPECT. IEEE Trans. Nucl. Sci., NS-4L1585-1595, 1994. [72] E. C. Frey and B. M. W. Tsui. A new method for modeling the spatially-variant, object-dependent scatter response function in SPECT. In 1996 Nuclear Science Symposium Conference Record, pages 1082-1086, November 1996. [73] D. Kadrmas, E . C Frey, and B.M.W. Tsui. Application of reconstruction-based scatter compensation to Tl-201 SPECT. In 1996 IEEE Nuclear Science Symposium Conference Record, pages 1663 - 1667, 1996. [74] D. J. Kadrmas, E. C. Frey, and B. M. W. Tsui. Application of reconstruction-based scatter compensation to thallium-201 SPECT: implementations for reduced reconstructed image noise. IEEE Trans. Med. Imag., 17:325-333, 1998. [75] D. J. Kadrmas, E. C. Frey, S.S. Karimi, and B.M.W. Tsui. Fast implementations of reconstruction-based scatter compensation in fully 3d SPECT image reconstruc-tion. Phys. Med. Biol, 43:857-873, 1998. [76] X. Song, E. C Frey, W.T. Wang, Y. Du, and B.M.W. Tsui. Validation and evalu-ation of model-based cross-talk compensation method in simultaneous 9 9 m tc stress and 2 0 1 t l rest myocardial perfusion SPECT. In 2002 IEEE Nuclear Science Sympo-sium and Medical Imaging Conference Record, volume 3, pages 1370-1374, Novem-ber 2002. [77] X. Song, E. C. Frey, X. He, W. P. Segars, and B. M. W. Tsui. A mathematical observer study for evaluation of a model-based compensation method for crosstalk in simultaneous dual isotope SPECT. In 2003 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, pages M9-4, October 2003. [78] C. Bai, G. L. Zheng, and G. Gullberg. A slice-by-slice blurring model and kernel evaluation using the klein-nishina formula for 3D scatter compensation in parallel and converging beam SPECT. Phys. Med. Biol., 45(5):1275-1307, 2000. [79] C. Bai, G.L. Zeng, and G.T. Gullberg. The modeling of multiple order compton scatter in SPECT. IEEE Trans. Nucl. Sci, 48:38-42, 2001. [80] J.W. Beck, R.J. Jaszczak, R.E. Coleman, C F . Starmer, and L.W. Nolte. Analy-sis of SPECT including scatter and attenuation using sophisticated Monte Carlo modeling methods. IEEE Trans. Nucl. Sci, 29(1):506-511, 1982. [81] M. Ljungberg and S-E. Strand. A Monte Carlo program for the simulation of scintillation camera characteristics. Comp. Meth. Prog. Biomed., 29:257-272, 1989. Bibliography 170 [82] D.R. Haynor, R.L. Harrison, T.K. Lewellen, A.N. Bice, CP. Anson, S.B. Gillispie, R.S. Miyaoka, K.R. Pollard, and J.B. Zhu. Improving the efficiency of emission tomography simulations using variance reduction techniques. IEEE Trans. Nucl. Sci, 37(2):749-753, 1990. [83] D.J. deVries, S.C Moore, R.E. Zimmerman, S.P. Mueller, B. Friedland, and R.C. Lanza. Development and validation of a Monte Carlo simulation of photon trans-port in an Anger camera. IEEE Trans. Med. Imag., 9(4):430-438, 1990. [84] H. Wang, R.J. Jaszczak, and R.E. Coleman. Solid geometry-based object model for Monte Carlo simulated emission and transmission tomographic imaging systems. IEEE Trans. Med. Imag., ll(3):361-372, 1992. [85] J . C Yanch, A.B. Dobrzeniecki, C. Ramanathan, and R. Behrman. Physically real-istic Monte Carlo simulation of source, collimator and tomographic data acquisition for emission computed tomography. Phys. Med. Biol, 37(4):853-870, 1992. [86] H. W. DeJong, M. Gieles, M. Ljungberg, and F. J. Beekman. Validation of a fast simulator for clinical SPECT simulation. J. Nucl. Med., 42:141P, 2001. Abstract. [87] H. W. DeJong and F. J. Beekman. Rapid SPECT simulation of downscatter in non-uniform media. Phys. Med. Biol, 46(3):621-635, 2001. [88] H. W. de Jong, E. T. P. Slijpen, and F. J. Beekman. Acceleration of monte carlo SPECT simulation using convolution-based forced detection. IEEE Trans. Nucl Sci, 48(l):58-64, 2001. [89] H.W.A.M. de Jong, W. Wen-Tung, E . C Frey, and F.J. Beekman. Efficient sim-ulation of spect down-scatter including photon interactions with crystal and lead. In 2001 IEEE Nuclear Science Symposium Conference Record, volume 3, pages 1517-1521, November 2002. [90] F.J. Beekman, H.W.A.M. de Jong, and S.van Geloven. Efficient fully 3-D iterative SPECT reconstruction with monte carlo-based scatter compensation. IEEE Trans. Med. Imag., 21:867-877, 2002. [91] H.W.A.M. de Jong, F.J. Beekman, and T.P. Slijpen. Acceleration of monte carlo spect simulation using convolution-based forced detection. In 1999 IEEE Nuclear Science Symposium Conference Record, volume 3, pages 1532-1536, October 1999. [92] H.W.A.M. de Jong and F.J. Beekman. Efficient photon cross-talk calculation in SPECT. In 2000 IEEE Nuclear Science Symposium Conference Record, volume 2, pages 13-30 - 13-34, October 2000. Bibliography 171 [93] F.J. Beekman, C. Kamphuis, and E. C. Frey. Scatter compensation methods in 3D iterative SPECT reconstruction: a simulation study. Phys. Med. Biol., 42:1619-1632, 1997. [94] T. Riauka and Z. Gortel. Photon propagation and detection in single-photon emis-sion computed tomography - an analytical approach. Med. Phys., 21:1311-1321, 1994. [95] R.G. Wells and A. Celler. Analytical calculation of scatter distributions in SPECT projections. In 1995 IEEE Nuclear Science Symposium Conference Record, pages 1087-1091, October 1995. [96] J. Sled, A. Celler, S. Barney, and M. Ivanovic. Monte Carlo simulation in SPECT: a comparison of two approaches. In Rodney Shaw, editor, Medical Imaging 1994: Physics of Medical Imaging. Proc. SPIE 2163, 1994. [97] C E . Floyd, R.J. Jaszczak, C C Harris, and R.E. Coleman. Revised scatter fraction results for SPECT. Phys. Med. Biol., 32(12):1663-1666, 1987. [98] J.A. Sorenson and M.E. Phelps. Physics in Nuclear Medicine. W.B. Sauders Company, 1987. [99] R. L. Harrison, D. R. Haynor, S. B. Gillispie, S. D. Vannoy, M. S. Kaplan, and T. K. Lewellen. A public-domain simulation system for emission tomography photon tracking through heterogeneous attenuation using importance sampling. J. Nucl. Med., 34(5):60P, 1993. [100] R. L. Harrison, S. D. Vannoy, D. R. Haynor, S. B. Gillispie, M. S. Kaplan, and T. K. Lewellen. Preliminary experience with the photon history generator: module of a public-domain simulation system for emission tomography. 1993 IEEE Nuclear Science Symposium Conference Record, pages 1154-1158, 1993. [101] H. Johns and J. Cunningham. The Physics of Radiology. Charles C. Thomas, Springfield, Illinois, fourth edition, 1983. [102] R.G. Wells, A. Celler, and R. Harrop. Analytical photon distribution calculations: a fast alternative to Monte Carlo simulations. Med. Phys., 23(5):810, May 1996. [103] R.G. Wells, A. Celler, and R. Harrop. Experimental validation of an analytical method of calculation of photon distributions. In 1996 IEEE Nuclear Science Symposium Conference Record, pages 1402-1406, November 1996. [104] R.G. Wells, A. Celler, and R. Harrop. Experimental validation of an analytical method of calculating SPECT projection data. IEEE Trans. Nucl. Sci, 44(3): 1283-1290, 1997. Bibliography 172 [105] P.H Pretorius, W. Xia, M.A. King, B.M.W. Tsui, T.S. Pan, and B.J. Villegas. Eval-uation of right and left ventricular volume and ejection fraction using a mathemati-cal cardiac torso phantom for gated pool SPECT. J. Nucl. Med., 38(10):1528-1534, 1997. [106] W. Cheney and D. Kincaid. Numerical Mathematics and Computing. Brooks/Cole, Pacific Grove, Calif., second edition, 1985. [107] S.Blinder, A.Celler, Wells G, D.Thomson, and R.Harrop. Experimental verification of 3D detector response compensation using the OSEM reconstruction method. In 2001 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, pages M13C-11, November 2001. [108] H.C. Gifford, M.A. King, R.G. Wells, W.G. Hawkins, M.V. Narayanan, and P.H. Pretorius. LROC analysis of detector-response compensation in SPECT. IEEE Trans. Med. Imag., 19:463-473, 2000. [109] Zhengrong Liang, Juihsi Cheng, and Jinghan Ye. Validation of the central-ray approximation for attenuated depth-dependent convolution in quantitative SPECT reconstruction. Phys. Med. Biol, 42(2):433-439, 1997. [110] J.C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright. Convergence prop-erties of the nelder-mead simplex method in low dimensions. SIAM Journal of Optimization, 9(1): 112-147, 1998. [Ill] K.L. Dixon, L.N. Baldwin, B. Coquinco, E.J. Vandervoort, A. Fung, and A. Celler A. A robust and versatile method for the quantification of myocardial infarct size. In IEEE Medical Imaging Conference, Portland, Oregon, November 2003. [112] W.S. Snyder, M.R. Ford, and G.G. Warner et al Estimates of absorbed frac-tions for monoenergetic photon sources uniformly distributed in various organs of a heterogeneous phantom. MIRD Pamphlet No.5. J. Nucl. Med., 10(Suppl.3):5-52, 1969.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Implementation of an analytically based scatter correction...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Implementation of an analytically based scatter correction in SPECT reconstructions Vandervoort, Eric 2004
pdf
Page Metadata
Item Metadata
Title | Implementation of an analytically based scatter correction in SPECT reconstructions |
Creator |
Vandervoort, Eric |
Date Issued | 2004 |
Description | Photon scatter is one of the main effects that degrade the quality and quantitative accuracy of SPECT images. The objective of the work performed in this thesis was to develop and validate a scatter correction technique that uses an accurate physics-based model of photon scatter as part of an iterative image reconstruction algorithm. This work was performed in two distinct stages. The first stage was to develop a scatter calculation based on a simplified version of an existing technique, the analytic photon distribution (APD) method. The scatter distributions generated using the approximate technique were compared to those obtained using the original APD method. On average the differences between the two methods were small relative to the statistical uncertainty in SPECT projection data (between 6 to 30 times smaller depending on the activity distribution) for a typical perfusion study. The second stage was to implement this scatter correction in an image reconstruction algorithm, and test its performance using computer simulated projection data, experimental data obtained from physical phantoms, and patient data. Images corrected for scatter, attenuation and collimator blurring were compared to images corrected only for attenuation and collimator blurring. In the simulation studies results were compared to an ideal situation in which only the primary (unscattered) photon data is reconstructed. In all cases, the scatter corrected images demonstrate substantially improved image contrast relative to no scatter correction. In simulations the scatter corrected and primaryonly images had similar contrast (to within 3%) and noise properties. Images with no scatter correction had contrast values 29% poorer on average than the ideal correction. Currently this scatter correction takes 3 to 4 hours for a typical clinical myocardial study. In the future additional code optimization and some improvements in computer processor speed could allow this scatter correction to be implemented clinically. |
Extent | 14272333 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-11-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0091460 |
URI | http://hdl.handle.net/2429/15176 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2004-0136.pdf [ 13.61MB ]
- Metadata
- JSON: 831-1.0091460.json
- JSON-LD: 831-1.0091460-ld.json
- RDF/XML (Pretty): 831-1.0091460-rdf.xml
- RDF/JSON: 831-1.0091460-rdf.json
- Turtle: 831-1.0091460-turtle.txt
- N-Triples: 831-1.0091460-rdf-ntriples.txt
- Original Record: 831-1.0091460-source.json
- Full Text
- 831-1.0091460-fulltext.txt
- Citation
- 831-1.0091460.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0091460/manifest