"Science, Faculty of"@en . "Physics and Astronomy, Department of"@en . "DSpace"@en . "UBCV"@en . "Krzywinski, Martin"@en . "2009-05-25T20:19:58Z"@en . "1998"@en . "Master of Science - MSc"@en . "University of British Columbia"@en . "An evaluation of new reconstruction algorithms for 3D positron emission tomography is presented. Using phantom and patient data, the standard 3D filtered backprojection algorithm is used to assess the merits of faster alternatives to full 3D reconstructions: FORE + 2D FBP, FORE + OSEM and FORE + SAGE. The phantom figures of merit that are used to evaluate these new algorithms are: resolution, noise, contrast recovery, signal-to-noise ratio, image\r\nuniformity and time activity curve variability. Using patient scans, the images are compared to 3D FBP reconstructions by examining the differences in radial profiles, region of interest activity, image root-mean-square, Logan slopes and computational demands. Both phantom\r\nand patient data analysis indicates that FORE + 2D FBP performs equally well as 3D FBP\r\nin terms of all figures of merit examined, with the exception that image uniformity and time activity curve variability are better than for 3D FBP. One clear benefit in using FORE .+ 2D FBP is that reconstruction is 10 times faster than 3D FBP. Based on the results obtained 3D FBP can be replaced by FORE + 2D FBP with no significant impact on subsequent analysis. Distribution volume ratios obtained with FORE + 2D FBP were different from 3D FBP results by 0.2 \u00C2\u00B1 0.4 % for cortical input function and \u00E2\u0080\u00941.5 \u00C2\u00B1 0.4 % for cerebellar input function. Both OSEM and SAGE performed well, with increased resolution and distribution volume differences being < 3% for either input function."@en . "https://circle.library.ubc.ca/rest/handle/2429/8150?expand=metadata"@en . "8436113 bytes"@en . "application/pdf"@en . "EVALUATION OF NEW ITERATIVE AND REBINNING RECONSTRUCTION ALGORITHMS IN FULLY THREE-DIMENSIONAL POSITRON EMISSION T O M O G R A P H Y By Martin Krzywinski B. Sc. (Physics-Mathematics) University of Ottawa, Ottawa, Ontario, 1996 A THESIS S U B M I T T E D IN 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 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 August 1998 \u00C2\u00A9 Martin Krzywinski, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this .thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: 11 A b s t r a c t A n evaluation of new reconstruction algorithms for 3D positron emission tomography is presented. Using phantom and patient data, the standard 3D filtered backprojection algorithm is used to assess the merits of faster alternatives to full 3D reconstructions: F O R E + 2D F B P , F O R E + O S E M and F O R E + S A G E . The phantom figures of merit that are used to evalu-ate these new algorithms are: resolution, noise, contrast recovery, signal-to-noise ratio, image uniformity and time activity curve variability. Using patient scans, the images are compared to 3D F B P reconstructions by examining the differences in radial profiles, region of interest activity, image root-mean-square, Logan slopes and computational demands. Both phantom and patient data analysis indicates that F O R E + 2D F B P performs equally well as 3D F B P in terms of all figures of merit examined, with the exception that image uniformity and time activity curve variability are better than for 3D F B P . One clear benefit in using F O R E .+ 2D F B P is that reconstruction is 10 times faster than 3D F B P . Based on the results obtained 3D F B P can be replaced by F O R E + 2D F B P with no significant impact on subsequent analysis. \" Distribution volume ratios obtained with F O R E + 2D F B P were different from 3D F B P results by 0.2 \u00C2\u00B1 0.4 % for cortical input function and \u00E2\u0080\u00941.5 \u00C2\u00B1 0.4 % for cerebellar input function. Both O S E M and S A G E performed well, with increased resolution and distribution volume differences being < 3% for either input function. i i i Contents A b s t r a c t i i List of Tables v i List of Figures v i i 1 Introduction 1 1.1 Goal of Thesis 1 1.2 Thesis Outline 2 1.3 Motivation 2 1.4 Previous Work 3 1.5 Positron Emission Tomography 4 1.5.1 Physics of P E T 4 1.5.2 Acquisition Mode: 2D and 3D P E T 8 1.6 Image Reconstruction Methods 14 1.6.1 Transform Methods 15 1.6.2 Algebraic Methods 16 1.6.3 Statistical Methods 18 1.6.4 Rebinning Methods 19 2 Iterative Algorithms 23 2.1 E M Algorithm 23 2.2 O S E M Algorithm 26 2.3 S A G E Algorithm 29 2.3.1 S A G E 1 29 2.3.2 S A G E 3 30 2.3.3 Regularized SAGE3 31 2.3.4 Notation 35 3 Figures of Mer i t for Algori thm Evaluation 36 3.1 Evaluation Philosophy 36 3.1.1 Phantom Studies 38 3.1.2 Digital Phantoms 39 3.1.3 human subject Brain Scans 40 3.1.4 Scanning Hardware 40 CONTENTS iv 3.2 Phantom Evaluation Criteria 40 3.2.1 Resolution 41 3.2.2 Noise 41 3.2.3 Contrast 42 3.2.4 Signal-to-Noise Ratio (SNR) 42 3.2.5 Axia l Uniformity 43 3.2.6 Radial Uniformity 43 3.2.7 In-plane Uniformity 43 3.2.8 Time Activity Curve Comparison 44 3.2.9 Calibration Factors 45 3.3 Human Subject Data Evaluation Criteria 45 3.3.1 Radial Profile Difference 45 3.3.2 Region-of-Interest Difference 46 3.3.3 Image Root-Mean-Square (RMS) 47 3.3.4 Kinetic Parameter Comparison 48 3.3.5 Computational Demands 52 4 Experimental Results 53 4.1 Choice of Algorithm Parameters 53 4.1.1 F O R E 53 4.1.2 2D F B P 53 4.1.3 O S E M 54 4.1.4 S A G E 3 54 4.2 Phantom Studies 54 4.2.1 Resolution 54 4.2.2 Noise 58 4.2.3 Contrast 63 4.2.4 Signal-to-Noise Ratio 64 4.2.5 Axia l Uniformity 64 4.2.6 Radial Uniformity 64 4.2.7 N E M A Uniformity 67 4.2.8 Time Activity Curve Comparison 67 4.2.9 Summary of Phantom Studies 67 4.3 human subject Studies 67 4.3.1 Radial Profile Difference 70 4.3.2 Region-of-Interest Difference 71 4.3.3 Image R M S 72 4.3.4 Kinetic Parameter Comparison 72 4.3.5 Computational Demands 73 5 Discussion and Conclusion 80 5.1 Algorithm Performance 80 5.1.1 F O R E and 2D F B P 81 5.1.2 F O R E and O S E M 81 5.1.3 F O R E and S A G E 83 CONTENTS 5.2 Algorithm Selection Bibliography vi List of Tables 4.1 Resolution for Various Reconstructions at the Center of the F O V . The parameter / represents the spatial filter and (3 is the regularization parameter for S A G E . . . 56 5.1 Comparison of Phantom FOMs for 3D F B P and F O R E + 2D F B P . Speed is given in minutes per frame 81 5.2 Comparison of Phantom FOMs for 3D F B P and F O R E + O S E M . A denotes the percent improvement in going from 3D F B P to F O R E and 2DFBP. Negative (positive) values of A signify a decrease (increase) in quality of the particular F O M . Speed is given in minutes per frame 82 5.3 Comparison of Difference in Phantom F O M s between 3D F B P and F O R E + O S E M . Negative (positive) values of A signify a decrease (increase) in quality of the par-ticular F O M 82 5.4 Comparison of Phantom FOMs for 3D F B P and F O R E + S A G E . Speed is given in minutes per frame 83 5.5 Comparison of Difference in Phantom F O M s between 3D F B P and F O R E + S A G E . Negative (positive) values of A signify a decrease (increase) in quality of the particular F O M 84 5.6 Comparison of difference in human subject study quantitation F O M s between 3D F B P and F O R E + 2D F B P / O S E M / S A G E 85 5.7 Comparison of difference in distribution volume ratios between 3D F B P and F O R E + 2D F B P / O S E M / S A G E 85 vii List of Figures 1.1 The event is registered as originating from along a line of response. The details of the decay and annihilation are shown in the figure inset 5 1.2 A point source gives rise to a sinusoidal sinogram 7 1.3 Scatter, randoms and attenuation all lead to a decrease in the number of true events detected 8 1.4 For each frame in a dynamic scan, the same ROI is placed on a given image plane. The activity within this region is found and plotted against time 9 1.5 Positioning terminology in P E T 10 1.6 The difference between 2D and 3D acquisition. The black star represents the actual decay location. The light star is the location as seen by the scanner after the event has been assigned to a plane 11 1.7 For sinograms beyond a certain azimuthal angle, 6 > 0 (see text), some projec-tions are missing, or incomplete. Here (j) is the standard polar angle 12 1.8 (1) Only the complete sinograms are used to reconstruct the image. (2) The missing projections are calculated based on the image. (3) A l l the sinograms are used to reconstruct the final image 13 1.9 The sinogram is a set of projections of the activity within the F O V taken at different angles 15 2.1 The effect of iteration number (here 3, 25 and 100 iterations are shown) on reconstructed image quality with the E M algorithm. After a certain number of iterations, the algorithm's performance is degraded by the presence of noise in the projections 26 2.2 The effect of number of O S E M iterations and subsets on image reconstruction. . 28 2.3 The effect of various values of f3 on the reconstruction of a uniform cylinder with S A G E : 32 2.4 The effect of number of S A G E iterations and regularization on final image quality. 34 3.1 The ROIs used for Logan analysis and ROI difference analysis 47 4.1 R M S of a single brain image plane between S A G E and 3D F B P 55 4.2 Variation in radial resolution with radial distance in the F O V 58 4.3 Variation in radial resolution with axial distance in the F O V 59 LIST OF FIGURES vi i i 4.4 Variation of axial resolution with radial distance in the F O V 60 4.5 Variation of axial resolution with axial distance in the F O V 61 4.6 Image noise for various reconstructions at two count rates. Note the scale differ-ence in the two plots 62 4.7 Hot spot contrast for various reconstruction algorithms 63 4.8 Axia l uniformity for various reconstructions at high count rates 65 4.9 Radial uniformity for various algorithms 66 4.10 N E M A uniformity for various reconstructions in two count regimes 68 4.11 Variability in the T A C for various reconstructions 69 4.12 Radial profiles taken vertically through a raclopride brain image. Several recon-structions are shown. The solid line represents 3D F B P and the scatter points are the algorithm indicated above the plot 70 4.13 Average radial profile difference for various reconstruction methods 71 4.14 Average normalized ROI differences for three ROI groups 75 4.15 Image R M S and signed difference for various reconstructions averaged over 10 human subject scans 76 4.16 Average difference in Logan slopes (DVR) between 3D F B P Hann 0.4 and various reconstruction methods. The error bars represent the standard deviation of the D V R among the 10 human subjects 77 4.17 Average Logan slope dispersion of the difference between 3D F B P Hann 0.4 and various reconstruction methods 78 4.18 Reconstruction times for standard 3D F B P and F O R E reconstructions. The time is significantly reduced because 2D reconstruction methods can be utilized after rebinning with F O R E 79 1 Chapter 1 Introduction 1.1 Goal of Thesis This thesis will present an evaluation and comparison of image reconstruction methods for 3D positron emission tomography (PET) data. The reconstruction schemes under study are based on rebinning a 3D data set to a 2D data set with Fourier Rebinning (FORE) , followed by 2D versions of one of the following algorithms: Ordered-Subsets Expectation Maximization (OSEM), Space-Alternating Generalized Expectation Maximization (SAGE) and 2D Filtered Backprojection (2D F B P ) . Based on the performance of these methods, as tested with phantoms and human subject studies, it will be possible to decide whether the current 3D reconstruction methodology, 3D Filtered Backprojection (3D F B P ) can be replaced. 3D F B P is disadvan-taged because it is computationally very costly. To motivate the replacement of the current reconstruction standard, any new method should yield images of at least equivalent quality, where quality is defined by a series of figures of merit. Factors such as speed, image quality and implementation should all be considered. CHAPTER 1. INTRODUCTION 2 1.2 Thesis Outline To understand the new algorithms and the significance of the tests applied to them, some background in P E T is required. A short summary of this imaging modality will be presented in this chapter, followed by a brief exposition of the general types of reconstruction algorithms available. Subsequent to the background, the candidate algorithms will be presented. A detailed understanding of them is not required to understand the results and therefore the presentation will be kept short. Again, the reader is referred to detailed references. The evaluation methodology is described next. The phantom and human subject scans figures of merit used to gauge performance are listed and described. The use of human subject studies is presented. The following results section presents the performance of each method. The discussion section will examine the results and propose which of the tested meth-ods provides optimal reconstructed images. 1.3 Motivation The work described in this thesis is motivated by the fact that (a) the current recon-struction standard in P E T is very computationally demanding and (b) new and still clinically untested algorithms exist which are both faster and have better noise characteristics. Recon-struction speed is of great importance and a fast algorithm is extremely desirable. Even with significant computer resources, the current reconstruction speed makes the processing of sever studies per day at a typical P E T center difficult to accomodate. For example, currently our center can process 48 frames per day, assuming constant computer uptime. It takes 10 minutes CHAPTER 1. INTRODUCTION 3 per frame to process the data for reconstruction and 20 minutes per frame to reconstruct the data to obtain an image. Scans are typically 16 frames long and it is not unusual to perform 3 such scans in a single day. This taxes the computers maximally and if a reconstruction needs to be redone, backlog can occur. Some scans have 25 frames and if two of these are done per day the computational requirements would exceed available resources. In addition to speed, decreased noise or better resolution with the same noise level would be highly valued, particularly for brain scans where the region of interest is anatomically small and contains relatively few counts. The algorithms which will be examined are both faster and offer lower noise images at equivalent resolution compared to the current standard. In order to replace the standard reconstruction algorithm, 3D filtered-backprojection, with confidence, any new method must be fully tested to ensure that the images can be used. Any such analysis of algorithms necessarily requires bringing to bear a full battery of phantom tests to establish the absence of any systematic artifacts, followed by human subject studies to ensure that the algorithms yield consistent results from biologically significant data. From the point of view of image analysis the switch in reconstruction, or for that matter any data analysis or methods, should be as transparent as possible, with no systematic impact on the data. 1.4 Previous Work During the period of 1997-1998, new algorithms have have been developed and refined for use in P E T reconstruction. Their acceptance into clinical and research centers as viable reconstruction methods has been slow largely due to the overwhelming use of a single recon-CHAPTER 1. INTRODUCTION 4 struction standard and lack of comparison in realistic scanning situations. Some studies exist [1] [2] [3] [4] [5] in which these algorithms have been used and compared but they are sparse and none include both phantom and clinical data. This thesis will provide the information required to answer the question whether these methods can be immediately used. 1.5 Positron Emission Tomography 1.5.1 Physics of P E T P E T , unlike some imaging modalities like M R I and X-rays, collects functional, rather than anatomical, information which can be used to probe into the biochemical function of specific cells or to extract information about kinetic parameters of metabolic pathways. Several excellent review articles exist [6] [7] [8] and the reader is encouraged to use them to supplement the introduction presented in this thesis. . The physical basis for P E T is fundamentally simple [9] [10] ([11], chapter 6). The human subject is injected with a chemical called a tracer which has been labeled with a positron-emitting radioactive isotope. The tracer is either a close analogue of the endogenous species present in the body, and can be used to study a metabolic pathway, or has a high specificity for a particular binding site whose population is of interest. For example, glucose metabolism in the body is studied with 1 8F-deoxygluxose. This compound is similar to glucose except that one of the protons has been replaced with the 1 8 F isotope and one of the - O H groups is missing. The 1 8 F allows the detection of the tracer in the body through the use of P E T and the missing - O H traps the tracer in a specific metabolic step. Once the tracer has had time to diffuse to the organ of interest and trapping has reached equilibrium, the human subject is placed in the CHAPTER 1, INTRODUCTION 5 Figure 1.1: The event is registered as originating from along a line of response. The details of the decay and annihilation are shown in the figure inset. detector and the tracer activity distribution is monitored with the scanner. Another common radioisotope that is used is n C , which has a half-life of 20 minutes, compared to 110 minutes for 1 8 F . A comprehensive list of labels and tracers can be found in chapter 9 of [11]. Figure 1.1 illustrates how the radioactivity, X(x,y), within the field of view (FOV) is detected. A P E T scanner is constructed from stacked rings of detectors. Figure 1.1 shows only one such ring and for illustration we will consider only those events coming from the plane of the ring. How events spanning more than a single ring are handled depends on the acquisition mode (vide infra). While in the body, the decay of the radiotracer nucleus is described by p ^ n + e+ + ise (1.1) where a proton decays to emit a neutron, n, and a positron, e + (anti-matter equivalent of an electron). A n electron neutrino, ve, is also emitted in the process[12]. The emitted positron CHAPTER 1. INTRODUCTION 6 subsequently annihilates with a nearby electron producing two photons. Because the positron is nearly at rest immediately before the annihilation, the total linear momentum of the e+e~ pair is almost zero and the two photons are emitted nearly at 180\u00C2\u00B0 to each other. There is about an average 2.5\u00C2\u00B0 off-set in this angle, shown as 9 in Figure 1.1. In addition, the decay and annihilation are not exactly spatially coincident, because positron travels about 2 mm, in the case of 1 8 F , from the decay site before annihilating 1. The two photons travel through the body and, if unhindered, are detected. The scan-ner logs an annihilation when two of its detectors register photons within a narrow coincidence time window. The position of the two detectors is the only information we have about the decay and therefore we only know the line along which the original decay must lie. This is called the line of reponse (LOR) and is parametrized as shown in Figure 1.1, where r is the distance of the L O R from the origin and is the angle that the L O R makes with the y-axis (the orientation of the definition of is arbitrary). Because there is a finite number of detectors in the ring, the number of L O R s is also finite. The number of events registered for a given L O R is histogrammed into a structure called the sinogram. For example, if we only had a point source in our field of view at (x, y), the only L O R s with events in them are described by the relationship r = x cos (f> + y svt\ (f> (1.2) The sinogram would therefore contain binned events along a single sinusoidal path shown in Figure 1.2. Compton scatter within the F O V is primarily responsible for scattering the photons 1This distance depends on the positron energy. The higher the energy the further the positron travels away from the decay before annihilation. CHAPTER 1. INTRODUCTION 7 Figure 1.2: A point source gives rise to a sinusoidal sinogram. away from their original LORs. Either the scattered photon is detected and the L O R is misreg-istered or multiple scattering, or possibly single large angle scatter, reduces the photon energy to below the set detector limit and the event is not registered at all. Typically as much as 30% of the detected events are scatters in a typical brain scan acquired without septa (vide infra). Unregistered events are broadly categorized as being due to attenuation. Figure 1.3 illustrates the principle of scatter. In addition to missing some events, the scanner can register events which did not actually take place. This can happen when two photons from two separate events in the scanner are detected within the coincidence time window. This is termed as randoms and pictured in Figure 1.3. Processing the sinogram involves correcting for randoms, deadtime, scatter, detector non-uniformities as well as attenuation. Once the sinogram is corrected so that it best reflects only the actual decays that took place in the scanner it is ready for reconstruction, which is the next step in data processing. Reconstructed data contain information about the tracer distribution in the object CHAPTER 1. INTRODUCTION 8 not detected Figure 1.3: Scatter, randoms and attenuation all lead to a decrease in the number of true events detected. of interest. When kinetics of metabolic processes are of interest, the fundamental object used for this analysis is the time activity curve (TAC) shown in Figure 1.4. In order to construct a T A C one needs to perform a dynamic scan. Dynamic scans are composed of sequential frames, with each frame being a single scan. Frames vary in number and duration depending on the tracer used and kinetics of the biological process under investigation. For example, the first frames might be short to sample rapid changes in the T A C and the last few frames might be longer, as is often the case when tracer distribution changes slowly at later times. A typical 16 frame study acquired in 3D mode (vide infra) contains about 16N2 sinograms, where N is the number of detector rings in the scanner, and the reconstruction time of such a scan is very long, typically longer than the scan duration. 1.5.2 Acquisition Mode: 2D and 3D P E T The example using Figure 1.1 was simplified in that it considered a single detector ring and decays occurring only in the ring plane. Indeed, one could perform a scan in this way but CHAPTER 1. INTRODUCTION 9 Figure 1.4: For each frame in a dynamic scan, the same ROI is placed on a given image plane. The activity within this region is found and plotted against time. many events would not be registered. At the decay site, photons are emitted isotropically and the chances that they will be detected by the detectors in the same plane is small. Sacrificing events like this yields low statistics images which are noisy and carry very little interpretable information for anatomical structures which are small and contain few counts, as is the case in many brain scans. Therefore, scanners are typically built as a stack of detector rings. Originally, P E T data were acquired in so-called septa-in or 2D mode. Between each detector ring were positioned heavy metal (high Z) septa which effectively absorbed any photons along a L O R which spanned more than several detector planes. For a scanner with N detector rings, the detector LORs were assigned to one of 2JV \u00E2\u0080\u0094 1 unique planes. N of these planes corresponded to direct planes which were exactly the detector ring planes. There were N \u00E2\u0080\u0094 1 indirect planes which were taken to lie between the direct planes. In this way the F O V was CHAPTER 1. INTRODUCTION 10 [ Positioning Terminology ~| tr&nsaxiai/direct sinogram Figure 1.5: Positioning terminology in P E T . segmented axially into slices. Events were binned into the direct planes if they were registered by detectors with a ring difference of either A = 0 or A = 2. Events sequestered into indirect planes were those from detectors with a ring difference of A = 1 or A = 3. The binned events are collected into sinograms, one for each of the 2N \u00E2\u0080\u0094 1 image planes. Some of this terminology is shown in Figure 1.5. The 2D acquisition method is shown in Figure 1.6. It can be seen that in 2D mode the scanner moves the axial position of any events associated with A > 0. The advantage of this is that all the data are segregated by the scanner hardware into 2N \u00E2\u0080\u0094 1 independent slices which can then be reconstructed with any 2D reconstruction method. More importantly, however, the number of scattered events is greatly reduced, because a large fraction of scattered events will be stopped by the septa as it is likely to span more than 3 detector rings. However, there is one main disadvantages with 2D P E T : a significant number of events are rejected (or even not detected) and the sensitivity is therefore sacrificed. Not long after 2D P E T was being routinely used, the motivation to increase the number CHAPTER 1. INTRODUCTION 11 Indirect plane ^ / direct pla: imnppoD rejected event direct plane events indirect plane events [ 3D Acquisition sinogram B \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 sinogram C sinogram D 1DDDDDDDDDDDDDDD all events accepted Figure 1.6: The difference between 2D and 3D acquisition. The black star represents the actual decay location. The light star is the location as seen by the scanner after the event has been assigned to a plane. of acquired counts by increasing sensitivity urged the development of 3D or septa-less P E T . In terms of scanner hardware, all that was required was to remove the septa and accept events with ring differences A > 4. In terms of reconstruction, an entirely new method had to be devised. The relatively slow transition from 2D P E T to 3D P E T was limited by the reconstruction issues which followed the removal of septa. The difficulty stemed from the presence of not 2N \u00E2\u0080\u0094 1 but N2 individual image planes. Normally the scanner may be configured not to not accept some extremely oblique image planes, such as those with A = N \u00E2\u0080\u0094 1. What compounded the problem of 3D reconstruction most of all is that some image planes had incomplete, also termed missing, projections. This is caused by the fact that the detector does not afford full coverage of the F O V , but rather has an axial acceptance angle of about a = 10\u00C2\u00B0, depending on the machine. This means that there are no sinograms for angles which are closer than | - a to the normal vector of the detector ring plane. However, even CHAPTER 1. INTRODUCTION 12 sinograms at smaller angles may have some of their projection bins unfilled due to insufficient coverage of the F O V . Figure 1.7 shows this phenomenon, where sinograms for angles 9 < 9 < a have missing projections. Figure 1.7: For sinograms beyond a certain azimuthal angle, 9 > O (see text), some projections are missing, or incomplete. Here tj> is the standard polar angle. The reason why missing projections are problematic is that it is not possible to recon-struct data which have a spatially variant point spread function (PSF) using Fourier methods. The extension of 2D F B P to 3D F B P was made difficult by this fact and until a reasonable way to solve this problem could be found, fully 3D reconstruction would have to wait. To solve the spatial variance problem one could potentially disregard all of the sino-grams for which 9 > \u00C2\u00A9 and accept only those without missing projections for reconstruction. Doing this, however, sacrifices sensitivity since valuable events are dismissed and nullifies the very benefit of 3D P E T . The methodology which was developed [13] centered around estimating detector rings Missing Projections CHAPTER 1. INTRODUCTION 13 Figure 1.8: (1) Only the complete sinograms are used to reconstruct the image. (2) The missing projections are calculated based on the image. (3) A l l the sinograms are used to reconstruct the final image. the missing data using the complete data available for 9 < O. One does this by reconstructing the low statistics image with only complete sinograms and then using this reconstructed image forward projects to fill in the missing data for sinograms 9 > O. Then, once all the sinograms are complete, the image is reconstructed for the second and final time. This is illustrated in Figure 1.8. Even without considering the implementation, this is clearly a tedious process in which nearly two complete reconstructions have to be performed for a single image. The com-putational effort is made entirely in order that sensitivity is not sacrificed by making use of as many detected events as possible. 3D P E T images, because of the increase in sensitivity, have less noise for the same amount of activity in the F O V when compared to 2D P E T . The standard reconstruction method for 3D P E T is 3D F B P . A good overall description CHAPTER 1. INTRODUCTION 14 of 3D P E T is given in [14]. Because this algorithm has heavy computational demands, a faster method to process a 3D data set has long been sought after. In the next section, a description of various algorithms is presented to contrast different philosophies in reconstruction. Immediately following that a more detailed examination of the algorithms that will be studied in this thesis is given. 1.6 Image Reconstruction Methods Image reconstruction describes the process of converting the projection data in the sinogram into a viewable image which represents tracer activity distribution, X(x,y,z), in the F O V . Since P E T does not collect direct information about \(x,y,z) but only the values of line integrals of this quantity 2, some method must be devised to invert the projection process (generally called backprojection). Ordinarily, recovering the image from a given projection is not possible because an entire dimension has been collapsed during the projection step This is shown in Figure 1.9. However, because projections from many different angles exist, each projection is correlated with all other ones and it is possible to find a way to invert P E T scan data into an image. methods, and statistical methods. Depending on the implementation, an algorithm can be either deterministic or iterative, in cases where it is not possible to formulate a purely deterministic approach. In addition, data may be manipulated with rebinning methods before undergoing 2This is the so-called Radon transform. (1.3) Reconstruction methods can be broadly categorized into transform methods, algebraic CHAPTER 1. INTRODUCTION 15 | Sinogram as Projections] Figure 1.9: The sinogram is a set of projections of the activity within the F O V taken at different angles. reconstruction. Each of these will be now briefly described, before the focus is put on the algorithms tested in this thesis. 1.6.1 Transform Methods 2D Filtered Backprojection Using the Central Section Theorem, which states that the Fourier transform of the projection p, z, A ) where (r, (f>) are the transaxial parameters shown in Figure 1.1 and z and A are the quantities ^ (axial position) and i \u00E2\u0080\u0094 j (ring dif-ference), respectively, given that the event was registered by the rings i and j. F O R E is built upon a different parametrization, namely /oo A(r cos0 \u00E2\u0080\u0094 t sin 0, r sin 4> \u00E2\u0080\u0094 tcosfi, z + t6)dt (1.9) -oo 4The striatum is an anatomical structure in the central brain region which contains neurons that use dopamine as a neutrotransmitter. These particular neurons are of interest in neurodegenerative diseases such as Parkinson's. C H A P T E R 1. I N T R O D U C T I O N 21 where 6 = tan 6, with 8 being the angle between the sinogram and the transaxial plane and X(x,y,z) is the activity distribution. The development of F O R E starts with considering the 2D Fourier transform of the sinogram, P(w,k,z,6) = p(p,(f),z,6)e-ik't,e-iajsd^ds (1.10) J-Rn JO where R& is the effective radius of the F O V and k and u are the Fourier index and radial frequency, respectively. Using the symmetry p(r, (f) + Ti,z,8) = p{\u00E2\u0080\u0094r, (j), z, \u00E2\u0080\u00948) one can write eq. (1.10) with the use of eq. (1.9) in cylindrical coordinates as P(cv,k,z,8) = f ^ j ^ J 0 R \u00C2\u00B0 \ ( p c o s ( 3 , p s m ( 3 , z - p 8 S m ( f > ) x e -**0 e -*\u00C2\u00BBp cos e-Wpdpdpd in eq. (1.12) can be taken separated out which yields the expression P(u;,k,wz,6) = 2n(-i)k f 0 R \u00C2\u00AB F(p,p,uz)e-We-*k*rctan^ r \u00E2\u0080\u0094 ^ r - ^ ( L 1 3 ) x Jfc (^ptuy/l + S'2u'j/uj'2^j pdpd/3 where Jfc is a Bessel function and F(p,8,uz)= \{p cos 13, p sin (3, z)e~lul*zdz (1-14) J-L/2 with L being the axial length of the F O V . F O R E is based on the exact rebinning formula obtained from eq. (1.13). We can set 8 = 0 to find the relationship between the Fourier transformed sinograms and their direct counterparts. This is given by P ( W ^ , W z , 0 ) = e - \" \u00C2\u00A3 a r c t a n ^ ? {ojJl + ^-,k,uz,8\ (1.15) CHAPTER 1. INTRODUCTION 22 The expression in eq. (1.15) is an analytical relationship, expressed in Fourier space, between 2D and 3D data sets. Its implementation is complicated and time consuming and can be approximated by expanding in a Taylor series in the variable a = ^f-. If we take the zeroth order approximation and set a = 0 we find 7>(w, k,toz, 6) = V (w, k,u>z, 0) + 0(a) (1.16) W which is exactly the SSRB algorithm. In the first order approximation one takes fcarctana = ka + 0(a3) (1-17) which leads to V(to,k,Loz,6) = e-*kaV(u,k,uz,0)+O(a2) (1.18) The inverse transform along z is easy to take since the exponential factor e~lka represents a shift in the z direction. In this way inverting eq. (1.18) yields the main F O R E result P(CJ, k, z, 6) = P(u, k,z- kS/tJ, 0) (1.19) Implementing eq. (1.19) is particularly straightforward and fast since no interpolation is required in the to and k directions and only a single interpolation along z is needed. In addition, unlike the exact rebinning result in eq. (1.15), the Fourier transform along z is not taken in eq. (1.19) and hence we do not need to know the data for all values of z. This obviates the need for the time consuming task of finding the missing projections. 23 Chapter 2 Iterative Algorithms In this chapter the algorithms that will be tested are described in more detail. The Expecatation-Maximization (EM) algorithm will be described, since it provides a foundation for all of the iterative methods under study. Subsequently, O S E M and S A G E will be presented. 2.1 E M Algorithm The theory of the E M algorithm has been presented in the seminal work of Dempster et al. in 1977 [20]. The essence of the method is as follows [18]: denote the measured projections by Pi (j \u00E2\u0080\u0094 1...P). Let the number of decays detected in projection tube i emitted from pixel j (i = 1...F) be given by riij. Strictly, n = n^- is not an observed quantity but is related to the observed projections by F 3=1 However, this quantity itself is not vital, as it can be related to the actual image through the expectation operator. Let the actual activity in pixel i be given by / j . From eq. (1.6) it can be CHAPTER 2. ITERATIVE ALGORITHMS 24 seen that E[mj} = fij = Lijfj (2.2) where E[ ] is the expectation operator and Lij represents the system matrix in eq. (1.6). Notice that n is the realization of a random vector whose mean is the image itself and it is not this vector, n, but rather the expectation of this vector which is sought. Because n represents Poisson variables, which are independent over both of the two indices, the probability that a particular set of values corresponds to the measured projection data is F P fnij j=li=l w where the external sum is over all possible combinations of riij such that eq. (2.1) holds. This sum is required because there is potentially more than one way to redistribute counts in all the bins in such a way that the measured projections are still the same. Although eq. (2.3) looks elaborate, particularly the vague definition of the outer sum, its explicit evaluation is not required. The E M method maximizes the function Z(f|n) = logL(f |n) (2.4) and the reader is referred to the original paper [18],[19] for the full development of the algorithm which gives the image on the k iteration in terms of the image from the previous iteration by f(fc) = f ( f c - i ) 0 jLT^p0Lf(fc-i)jj k = 1:2,... (2.5) where (g> and 0 signify component-wise multiplication and division. Eq. (2.5) can also be CHAPTER 2. ITERATIVE ALGORITHMS 25 presented as \u00E2\u0080\u00A2(0) n = 1...P fo r k = 1,2,... F Pn \u00E2\u0080\u0094 Lnjfj 3=1 n = 1...P (2.6) fo r i = l...F The E M method depicted in eq. (2.6) has certain desirable constraints already built in. Mainly, non-negativity of the image and perservation of counts from one image guess to another. The former is particularly attractive considering that F B P methods, because of a finite frequency bandwidth in the filters, produce negative pixel values in the image. On the other hand, the E M method is disadvantaged because, although it converges to a unique solution in the case of noiseless projections, it fails to do so for real, noisy data. After a certain number of iterations, the E M algorithm begins to fit the image exclusively to the noise in the data and the quality of the image deteriorates, as shown in Figure 2.1. For this reason, the number of iterations for this method must be carefully controlled and the process must be stopped before this deterioration occurs. tively slow. Many attempts at increasing E M convergence have been made [30],[31],[32]. The E M converges slowly because during the entire execution of eq. (2.5) projections calculated from the old image are used. Even when the last pixel of the new image is calculated, it is with the projections of the old image. This can be rephrased by saying that the maximization data space is which is exceedingly large [33]. Second, r e g u l a t i o n is not possible E M suffers from two other pratical disadvantages. First, the convergence rate is rela-CHAPTER 2. ITERATIVE ALGORITHMS 26 Figure 2.1: The effect of iteration number (here 3, 25 and 100 iterations are shown) on re-constructed image quality with the E M algorithm. After a certain number of iterations, the algorithm's performance is degraded by the presence of noise in the projections. to include analytically into the E M framework [34]. Regularization is a method to penalize certain types of images in an attempt to prevent the algorithm guesses from deteriorating after a certain number of iterations. Regularization will be discussed in conjuction with the SAGE algorithm. 2.2 O S E M Algorithm The slow convergence properties of E M urged modifications to the method. One such alteration is provided by the the ordered-subsets methodology, embodied in the Ordered-Subset Expectation Maximization (OSEM) algorithm. Instead of updating the entire image with the old projections, as done in E M , the projections are divided into subsets and the entire image is updated once for each subset with the projections used for the subsequent subset having been updated to include changes in the previous image subset [35]. If the projection subsets are given by Sm (m = 1...M) then a single OSEM step is given by CHAPTER 2. ITERATIVE ALGORITHMS 27 \u00C2\u00A3 ( f c m o d i V , m ) f(kmodN,m\u00E2\u0080\u00941) J2 L n i ( P 0Lf( f c m o d J V > m - 1 ) ) Jfc = 1,2, m = 1.. .M (2.7) The notation f ( . k m o d N > m ) [s meant to represent the image at iteration (kmodN) after m number of updates during that iteration with the provision that f k ' \u00C2\u00B0 is equivalent to J f c - 1 > M . In terms of pseudo-code, (0,0) Pn \u00E2\u0080\u0094 \u00C2\u00A3 Lnjfj 3=1 fo r k = 0,l,... fo r m = 1...M F n = l...P - r f(k,rn\u00E2\u0080\u0094l) Q Pn \u00E2\u0080\u0094 Z^i ^njjj fl fc \u00C2\u00A3>n 3=1 (2-f o r \u00C2\u00AB = l . . . F J-\"ni -n\u00C2\u00A3Sm Pn-fi = * f i One O S E M iteration, which is defined as the traversal of all projection subsets, is equivalent in speed to a single E M iteration. However, O S E M updates the image once for each subset so that in a single iteration the image may be updated as often as 16 times, or more, although practically the algorithm provides optimum results for some subset number, usually smaller than 16. O S E M provides an increase of the rate of convergence of about N/p where N is the number of projections.and p is the size of each subset. How the subsets are chosen reflects on the convergence rate. Commonly, orthogonal subsets are used so that each subsequent subset contains as much new information as possible. The free parameters in O S E M are the number of subsets, S, and the number of iterations, / . Their effect on the reconstructed image is shown in Figure 2.2. Because O S E M , and E M in C H A P T E R 2. I T E R A T I V E A L G O R I T H M S s u b s e t s 12 JttJlMk, ggg|\u00C2\u00A3 .-aflMtts-. jfi^fe. 16 Figure 2.2: The effect of number of O S E M iterations and subsets on image reconstruction. CHAPTER 2. ITERATIVE ALGORITHMS 29 general, are not linear, filtering or regularization is required. In the case of O S E M , regularization cannot be analytically included, so a posteriori (post reconstruction) filtering is done with an isotropic Gaussian filter. The width of the filter is another O S E M parameter. 2.3 S A G E Algorithm Space-Alternating Generalized Expectation-Maximization (SAGE) has several vari-ants, such as SAGE1 , S A G E 2 and SAGE3 ([33]), coined in the manner by the original authors. SAGE1 is described because it is a relatively straightforward evolution of the E M algorithm, although not necessarily the one that offers the most benefit. SAGE2 offers intermediate bene-fit, but since SAGE3 , the most optimized S A G E method, is tested in this thesis SAGE2 is not described. 2.3.1 SAGE1 Statistical models can be used to calculate not only the true decays but also scatter, randoms and attenuation, since statistical models can be developed for each of these. Prior studies have demonstrated, however, that the more variables E M was asked to calculate, the slower the convergence rate was. The general statement is that small and less informative data spaces result in faster convergence [36] [37] [38]. In this case, the data space is the set of parameters over which the maximization step is performed. O S E M has a smaller data space than E M , because the maximization uses only a subset of projections rather than the full set. S A G E attempts to decrease the size of this data space further by updating the image one pixel at a time. If we index the pixels by i = 1...F, the simplest S A G E method [33] is therefore given by CHAPTER 2. ITERATIVE ALGORITHMS 30 k = 1 2 j(kmodF,i) _ y.(fcmodF,i-l) \" j ^ T ^ Q ^ ( f c m o d F . i - l ) ^ ' ' ' ' ' ' ^ Q) 1 = 1...F One iteration would consist of updating each pixel individually and using this changed pixel in the calculation of p 0 L f for the next pixel. SAGE1 can also be written as Pn=ELnJri0) n = l...P 3=1 for k = 1,2,... fo r j = 1...F P r> r- - V T \u00E2\u0080\u0094 C0 - ^n3 ^ n=L fn J 3 - C3h (2.10) Pr = Pn + ( / j f c ) - f f ~ l ) ) L n j Vn : L n j \u00C2\u00B1 0 2.3.2 SAGE3 Increased rate of convergence is obtained by choosing small data spaces. For SAGE1 the data space over which the maximization is performed is {nij}f=1 where j corresponds to the current pixel that is being updated. Other data spaces can be created. S A G E 3 [33] which is the variant that is being studied uses y ] L i n ft (fc) (fc) \u00E2\u0080\u00A2 n^j z\ ' = mm 3 '-wo Li zV-PoissoniLijifj + Zj)) ( 2 ' U ) B\f ~ PoisSOn Yl L i n f n - LijZj J WJ J The data space is {z\^\B^}n=1 and depends on the iteration. The quantity Pi = Zij+Bii (2.12) CHAPTER 2. ITERATIVE ALGORITHMS 31 has the correct distribution and the maximization can be done over the new data space. This new space attempts to merge the information from the image and its possible background (due to randoms, for example) to minimize the information content. Furthermore, eq. (2.11) shows that the background equivalent for pixel j, Zj, is taken to contain the values of all the other pixels. This is not intuitive but results in a data space responsible for faster convergence than by just including, for examples, randoms as the background for pixel j. S A G E 3 is nearly idendical to S A G E 1 in eq. (2.10) except that the step / j f c ) = c J / j f c _ 1 ) is replaced by / f = [c, ( / f \" X ) + Zf-V) - zf~\ ( 2 . 1 3 ) x x > 0 0 x < 0 where [x], denotes the operator [x], = < 2.3.3 Regularized SAGE3 As mentioned, in an effort to further increase convergence, other data spaces have been found. In addition, in order to prevent the classic deterioration of the image after many iterations as observed in E M , regularization is included. Because S A G E data spaces are small, adding regularization is possible analytically, without resorting to expensive numerical proce-dures during reconstruction [34]. The principle of regularization is to maximize the function /'(f) = /(f) \u00E2\u0080\u0094 (3R(f) rather than /(f), where R(f) is some measure of image roughness and (3 is a regularization parameter which weighs how much regularization is to be included. The range of (3 varies with the units used and may range from 2 - 2 0 to 2 2 0 . For our units we find that 2 - 9 < P < 2 - 3 is a useful range which we sample in differing powers of 2. Image roughness can be measured in many ways [39] [40] and a standard way is to do so by using differences of CHAPTER 2. ITERATIVE ALGORITHMS 32 Figure 2.3: The effect of various values of (3 on the reconstruction of a uniform cylinder with S A G E . neighbouring pixels, a w ^ E E f ^ - / / (2-14) where wzj are weight factors, usually inversely proportional to the distance between pixels i and j , and Ni denotes the appropriate indices for pixels in the neighbourhood of i. Any image which is locally rough, has a higher value of R(f) and thus low value of /'(f). Since the method maximizes /'(f) such noisy images will be penalized and the algorithm will tend to converge to images with a small value of R(i). This is illustrated in Figure 2.3, where a uniform cylinder is reconstructed with various amounts of regularization. In this figure, images with j3 < 0 were done with 10 iterations and all other images with 25 iterations. From Figure 2.4 it is evident that for extensive regularization, when j3 > 0 for example, convergence towards a final image is slow. The inclusion of regularization obviates the need for filtering. The present S A G E method that is tested is P M L - S A G E 3 (Penalized Maximum Likelihood) CHAPTER 2. ITERATIVE ALGORITHMS 33 [34]. The algorithm is given by Pn= n = l...P j = l F Wj- = E Wji i = l f o r k = 1,2,... f o r j = 1...F E t Fn (2.15) n = l Pn v 7 P (fc) _ , 'r i = i Pi Ak) zi = \u00E2\u0080\u00A2 m i 3 n -F~-fj ff = (-Bj + yJBj+Cjij^+Zj) /Pwr - zd Pn=Pn+ ( / j f c ) - f f i n j Vn : L n i ^ 0 Notice that the parameters Bj and are dependent on the iteration. The main attraction of S A G E is that regularization can be included analytically. This is due to the fact that the data space is small and equations that used to be coupled in regularized E M are now independent. The apparent increase in complexity in going from eq. (2.10) to eq. (2.15) is not accompanied by an increased computational demand. In fact, efficient implementations of these two variants of S A G E have execution times which are indistinguishable [33]. Wi th any S A G E algorithm, the adjustable parameters are the number of iterations and, if regularization is included (vide infra), one also chooses the value of the regularization parameter, j3. The impact of these two parameters on the final image is shown in Figure 2.4. CHAPTER 2. ITERATIVE ALGORITHMS 34 iterations 10 25 50 100 -9 -3 if t* the slope yields the distribution volume, D V . In a similar manner, one can write eq. (3.20) for the non-receptor region JpCB(t)dt _ J?Cp(t)dt _ 1_ CB(T) ~ A CB(T) k2 [ 6 Uj where CB represents the activity measured by the scanner in receptor-free region. For raclo-pride, for example, the receptor-free region can be either the cortex or the cerebellum (vide infra). Here X = k\/k2 for the receptor-free region. If one isolates JQ Cp(t)dt in eq. (3.20) and substitutes this into eq. (3.19) the result is /\"J A(t)dt J n T CB(t)dt + ^P-A(T) A(T) + { 6 - Z i ) where DVR = DV/X = B'max/K'd + 1 where B'max is the free receptor concentration (receptors CHAPTER 3. FIGURES OF MERIT FOR ALGORITHM EVALUATION 51 not occupied by endogenous neurotransmitter) and K'd is the effective equilibrium dissociation constant (K'd = K&I'/NS where f^s is the free fraction of tracer). [48]. Logan Analysis Implementation The Logan-type bloodless analysis described above was performed on human subject scans using cerebellar and tissue input functions. 8 regions of interest were placed in the striatum with 6 background tissue regions and 3 cerebellar regions. The D V R values were calculated for each method and directly compared to the results from 3DFBP. 10 human subjects were used for this comparison. Each of these scans used one of four tracers. In total there were 7 raclopride scans and one scan each with Sch 23390, tetrabenazine and methylphenidate. The Logan slope was calculated using both cerebellar and cortical input functions for each of the 8 ROIs in the human subject scans. The slope values were then averaged over the ROIs for a given human subject to yield the average value Asp = ^ s 3 D F B P \u00E2\u0080\u0094 S P ) R O I -The dispersion in this value was taken to be the standard deviation. Once this value was found, it was averaged over all human subjects, A s = ( A s p ) p and the error was taken to be the standard deviation. The cortical input is taken from the same planes as the striatial activity and when A s c o r is used as a F O M we are testing whether the algorithm reconstructs various regions of the same plane properly. When A s c e r is used, where the cerebellar input is taken from different planes than the striatal activity, we are checking whether the algorithm preforms correctly across many planes. Figure (3.1) shows how the different ROIs where placed. We seek A s to be as close to zero as possible for both cortical and cerebellar inputs with a small standard deviation. We report A s and its error, which measures patient-to-patient CHAPTER 3. FIGURES OF MERIT FOR ALGORITHM EVALUATION 52 difference, as well as the average dispersion in A s p , which measures ROI-to-ROI differences for each human subject. 3.3.5 Computational Demands The speed of the algorithm was measured by the time taken to reconstruct a scan. Image reconstruction is generally composed of two separate steps: correction and reconstruction. In the correction step, the data is corrected for attenuation, scatter and normalization. In the reconstruction step, the corrected sinograms are processed towards an image. A l l of the processing is done at our center on an Ultrasparc I (128 M R A M ) and a CSPI Supercard 1. The Ultrasparc performs all calculations, except for (i) attenuation correction (ii) 2D F B P and (iii) 3D F B P , which are carried out on the Supercard. 53 Chapter 4 Experimental Results 4.1 Choice of Algorithm Parameters 4.1.1 F O R E F O R E requires several parameters to be specified. These are mainly cutoff values which divide the Fourier frequency space into three distinct regions in which a different rebinning step is being performed. The presence of these parameters is closely related to the fact that the implementation of F O R E is based on the first-order approximation, rather than the exact rebinning. It was found that varying these parameters did not influence the image in any appreciable way and therefore the default parameters, as optimized by the original authors of the F O R E code, were used. 4.1.2 2D F B P The current reconstruction method is 3D F B P with a Hann filter with 0.4 cutoff. In order to mimic this as closely as possible with F O R E and 2D F B P , the same filter was used in CHAPTER 4. EXPERIMENTAL RESULTS the 2D F B P step. 54 4.1.3 O S E M The number of iterations and subsets were varied in order to find the optimum choice. A single iteration was seen to be insufficient - in this case not enough passes were made through the projections to fully reconstruct an image. In particular, the radial uniformity for a single iteration showed a marked decrease towards the edge of the cylinder. This effect was not seen for two or more iterations. In some respect, Figure 2.2 can be used to visually estimate the parameter space which is of use. 4.1.4 SAGE3 Because regularized S A G E is used, the number of iterations has a natural minimum that should ensure convergence and the maximum number of iterations is simply defined from a practical sense of computing time. Examining image R M S for various iterations and values of (3, shown in Figure 4.1, shows that about 5-6 iterations is sufficient to reach a stable S A G E image. 4 . 2 Phantom Studies 4.2.1 Resolution The point source phantom gives an indication of the reconstructed resolution, as men-tioned in Section 3.2.1. We are interested in two main trends here: (a) resolution at the center of the F O V , which gives a general impression of how well the system 1 will resolve small ob-l rrhe system here refers to both hardware and software, since the camera itself has an intrinsic resolution. CHAPTER 4. EXPERIMENTAL RESULTS 55 Image RMS fer Various SAGE Iterations 0 2 4 6 8 10 iterations Figure 4.1: R M S of a single brain image plane between S A G E and 3D F B P . jects and (b) variation of the resolution across the F O V . Table 4.1 shows the resolution of the algorithms at the center of the F O V as measured using the (0,0) cm point source. It was noticed that F O resolution did not depend on the number of iterations and subsets if there were more than 8 subsets. The table shows F O resolution for S = 12 and 1 \u00E2\u0080\u0094 2. The / parameter is the filter and unless otherwise given it is isotropic and Gaussian. FS parameters are log 2 (3 and 10 iterations were used for each reconstruction. Two of the S A G E reconstructions (log 2 (3 = \u00E2\u0080\u0094 6 and \u00E2\u0080\u00945) have been additionally filtered, to examine the effect of a posteriori filtering added to regularization. Center of F O V volume resolution for 3D F B P and F2D were nearly identical. Since F O R E is a rebinning algorithm, some compromise in axial resolution is expected but only minimally at the center of the F O V . Rebinning algorithms generally do very well at the center of the F O V because the axial mispositioning of events is smallest in this region. The volume resolution is reported in Table 4.1 to facilitate comparisons. Resolution is an important figure of CHAPTER 4. EXPERIMENTAL RESULTS 56 F W H M Parameters Radial (mm) Tangential (mm) Axia l (mm) Volume (mm 3) F O R E / O S E M S = 12,1 = 2 / = 0 mm 5.4 4.3 5.2 120 / = 2 mm 5.8 4.8 5.4 150 / = 3 mm 6.5 5.6 5.9 220 / = 4 mm 7.2 6.3 6.8 310 / = 5 mm 7.9 7.0 7.7 430 / = 6 mm 8.5 7.7 8.8 580 / = 7,9,2 mm 9.5 10.1 5.4 520 3D F B P Ramp 0.5 5.7 6.7 5.4 210 Hann 0.4 9.5 10.4 5.6 550 Hann 0.5 8.4 9.2 5.5 420 F 0 R E / 2 D F B P Ramp 0.5 6.5 6.0 5.3 210 Hann 0.4 10.1 9.9 5.5 550 Hann 0.5 8.8 8.7 5.4 420 F O R E / S A G E I = 10 log2(3 = - 9 5.5 4.9 5.2 140 log2(3 = - 6 6.2 5.4 5.4 180 log2f3 = - 5 6.4 5.7 5.4 195 = - 4 6.5 5.7 5.4 205 log2P = - 3 6.6 5.8 5.5 210 = 0 6.74 6.0 5.4 220 log2(3 = \u00E2\u0080\u00945, / = 3 mm 7.6 7.3 6.0 330 Zog2/3 = \u00E2\u0080\u00946, / = 4 mm 8.0 7.6 6.8 410 Table 4.1: Resolution for Various Reconstructions at the Center of the F O V . The parameter / represents the spatial filter and (3 is the regularization parameter for S A G E . CHAPTER 4. EXPERIMENTAL RESULTS 57 merit, but cannot be used by itself to evaluate an algorithm. Usually extremely high resolution reconstructions are noisy and are of little use clinically. Filtering or regularization is usually applied to decrease the noise. The resolution results will serve together with noise and human subject data to choose the best algorithm. For FO, the filter fx \u00E2\u0080\u0094 7 mm fy = 9 mm and fz = 2 mm was used to obtain a resolution-matched reconstruction compared to 3D F B P with Hann 0.4. This filter can be used to compare the algorithms on an equal footing at the resolution currently used with 3D F B P . How the resolution varies across the F O V is an important consideration, particularly because a rebinning algorithm is being applied. Figures (4.2) and (4.3) show how the radial resolution varies with radial and axial distance in the F O V , respectively. Both radial and axial dependent trends are the same for 3D F B P as for algorithm which uses F O R E . Figures (4.4) and (4.5) show how axial resolution varies within the F O V . It can be seen that axial resolution worsens as one moves radially out from the center of the F O V by about 1 mm in going from x = 0 cm to x = 10 cm. This trend is the same for 3D F B P as it is for F O R E , showing that the rebinning done by F O R E does not alter the variability in axial resolution. Resolution in water, a scatter medium, was also measured in the same way but was found to be identical to that measured in air. Furthermore, the profiles obtained from the point sources, in both air and water, were smooth and well suited for Gaussian fits.The resolution results show that the new algorithms reconstruct with the same consistency as 3D F B P within the F O V and show which parameters can be used to obtained a desired resolution. CHAPTER 4. EXPERIMENTAL RESULTS 58 Radial Dependence ef Radial Resolution \u00E2\u0080\u00A2a =9 \"i r 0 5 10 Radial distance from center of F0V (cm) O 3DFBP Ramp 0.5 \u00E2\u0080\u00A2 3DFBP Hann 0.4 \u00E2\u0080\u00A2 F0RE/2DFBP Ramp 0.5 \u00E2\u0080\u00A2 F0RE/2DFBP Hann 0.4 \u00E2\u0080\u00A2 F0RE/0SEM 1=2 S=12 + F0RE/SAGE 1=10 log2p=-6 Figure 4.2: Variation in radial resolution with radial distance in the F O V . 4.2.2 Noise Figure 4.6 shows the amount of noise in each reconstruction for both low (L) and high (H) counts using the experimental conditions outlined in Section 3.2.2. The current 3D F B P Hann 0.4 method has a noise value of 0.12 (H) and 0.58 (L), where (H) and (L) indicates high and low counts, respectively. Notice that the ramp filter, although offering a better resolution, produces images with 210% (H) and 320% (L) more noise. F2D images have the same amount of noise as those obtained with 3D F B P . This fact, together with the identical resolution shown in Table 4.1 suggests that these two reconstruction CHAPTER 4. EXPERIMENTAL RESULTS 59 Axial Dependence of Radial Resolution 10 -4 Cu 7 \u00E2\u0080\u00941 2 -3 5 5 -4 i r 0 2 4 Axial distance from center of P0V (cm) O 3DFBP Ramp 0.5 \u00E2\u0080\u00A2 3DFBP Hann 0.4 \u00E2\u0080\u00A2 F0RE/2DFBP Ramp 0.5 \u00E2\u0080\u00A2 FGRE/2DFBP Hann 0.4 A F0RE/0SEM 1=2 S=12 + FORE/SAGE 1=10 log2p=-6 Figure 4.3: Variation in radial resolution with axial distance in the F O V . methods will perform similarly, if not identically, with further phantom tests. Of course, F2D must still be tested with biologically meaningful human subject data. O S E M noise is seen to vary with both subsets and iterations. Commonly with E M algorithms, the fewer passes one makes through the data the less sharp the reconstructed image. In the extreme case of S = I = 1, the noise properties are excellent but, of course, the reconstruction is far from complete. Based on radial uniformity studies, we have found that S = 8 with one iteration provides incomplete reconstruction, whereas S = 12 with one iteration gives better results. For this reason, we do not test 5 = 8 but only S = 12 and 16 with 2 and CHAPTER 4. EXPERIMENTAL RESULTS 60 Radial Dependence ef Axial Resolution fa \u00E2\u0080\u00A2a ~i 1 r 0 5 10 Radial Distance from center of F0V (cm) O 3DFBP Ramp 0.5 \u00E2\u0080\u00A2 3DFBP Hann 0.4 A F0RE/2DFBP Ramp 0.5 \u00E2\u0080\u00A2 F0RE/2DFBP Hann 0.4 A F0RE/0SEM 1=2 S=12 + FORE/SAGE 1=10 log2P=-6 Figure 4.4: Variation of axial resolution with radial distance in the F O V . 4 iterations each. Figure 4.6 illustrates that noise is proportional to the product S \u00E2\u0080\u00A2 I so that (5,7) = (12,4) has less noise than (16,4) but more than (16,2). Because F O should be filtered, the noise of filtered images is also included in the figure. As expected, a posteriori filtering reduces the image noise substantially, although at the cost of resolution. To make the comparison between 3D F B P and F O one can compare either 3D F B P with Hann 0.4 to resolution matched F O . Figure 4.6 shows that high count F O image noise is 40% lower for (5,1) = (12,2) and 14% lower for (12,4) in comparison to 3D F B P with Hann 0.4. For low count rates, the noise is 38% lower for (12,2) and 26% lower for (12,4). Alternatively, CHAPTER 4. EXPERIMENTAL RESULTS 61 Axial Dependence of Axial Resolution 0 2 Axial distance from center of FOV (cm) O 3DFBP Ramp 0.5 \u00E2\u0080\u00A2 3DFBP Hann 0.4 \u00E2\u0080\u00A2 F9RE/2DFBP Ramp 0.5 \u00E2\u0080\u00A2 FGRE/2DFBP Harm 0.4 A F9RE/0SEM 1=2 S=12 + FORE/SAGE 1=10 log2fJ=-6 Figure 4.5: Variation of axial resolution with axial distance in the F O V . one can see that even unfiltered O S E M reconstructions have lower noise than 3D F B P with ramp 0.5, illustrating the improved noise handling of this algorithm. S A G E noise largely depends on regularization. Extremely regularized images, such as those with log 2 (3 > 0, have very little noise but also have very little detail as a result of the large regularization parameter. It was found that for log 2 (3 > \u00E2\u0080\u00943 the images are too smoothed to be used, even though the point source reconstruction yielded an acceptable resolution 2. In the high count regime, FS(\u00E2\u0080\u00944,10) (volume resolution V = 200 mm 3 ) has image noise at par with Because regularized iterative algorithms are not linear, one cannot directly predict the smoothness of the image only from the resolution found using a point source. CHAPTER 4. EXPERIMENTAL RESULTS 62 Noise Measurements High Counts Low Counts I I I I I I I I I I I I I I I I I 3DFBP 2DFBP OSEM SAGE \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 n O \u00E2\u0080\u00A2 g A A O \u00E2\u0080\u00A2 \u00E2\u0080\u009E \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 I I I I I I I I I I I I I 3DFBP 2DFBP OSEM T T T SAGE Filtering (mm) \u00E2\u0080\u00A2 Nene a 7,9,2 o 3,3,3 A 4,4,4 Figure 4.6: Image noise for various reconstructions at two count rates. Note the scale difference in the two plots. FO(12,2,3) (V = 216 mm 3 ) and FO(12,4,4) (V = 312 mm 3 ) . FS(-5 ,10) (V = 195 mm 3 ) has image noise of 0.13 which compares to unfiltered FO(12,2,0) (V = 122 mm 3 ) . Once FS(-5,10) is filtered with a 3 mm filter (V = 330 mm 3 ) the noise drops to the same value as FS(\u00E2\u0080\u00944,10). Less regularized FS(\u00E2\u0080\u00946,10) has relatively high image noise and needs to be filtered with a 4 mm filter to bring it down to the same level as FS(\u00E2\u0080\u00944,10) or FO(12, 2,3). Results in the low count regime, have the same relative pattern except that the absolute noise is increased about 10 fold. The noise results show that both S A G E and O S E M can be tuned to have similar noise properties with the same volume resolution. Furthermore, both iterative schemes offer lower noise images at equal or better resolution than 3D F B P . CHAPTER 4. EXPERIMENTAL RESULTS 63 4.2.3 Contrast Constrast results, measured as outlined in Section 3.2.3 are shown in Figure 4.7. The Contrast Ratio Measurements 5 4 g 3 2 -4 o O 1 High Counts Low Counts \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00C2\u00B0 \u00E2\u0080\u00A2 I I I I I I I I I I I I I I I I %? v \u00E2\u0080\u00A2 3DFBP 2DFBP 0SEM I I I I I I I I I I I I I I I I SAGE %? V J^ '7 3DFBP 2DFBP 0SEM SAGE Filtering (mm) \u00E2\u0080\u00A2 None \u00E2\u0080\u00A2 7,9,2 o 3,3,3 A 4,4,4 Figure 4.7: Hot spot contrast for various reconstruction algorithms. C R is quite constant at high counts, except for the excessively regularized case of log 2 /3 = 0. At low counts, S A G E overestimates the C R significantly. A l l the other methods have C R values very similar to 3D F B P . The theoretical value for C R for a 5:1 hotspotbackground ratio is 4 and because of partial volume effects the measured values are about 3.2. Also notice in Figure 4.7 that filtering reduces the CR. This is expected for small hotspots, since the activity is smeared outside of the hotspot during filtering. CHAPTER 4. EXPERIMENTAL RESULTS 64 4.2.4 Signal-to-Noise Ratio Using noise and contrast data, the SNR values for each reconstruction can be calcu-lated. However, because the C R is fairly constant for all the methods, SNR will be inversely proportional to noise and the SNR F O M does not present new information. Highly regularized S A G E images which have very small noise, such as those with log 2/? = 0, will have very high SNR, even though the C R is lower. However, because the image is very blurred they cannot be used. 4.2.5 Axial Uniformity Axia l uniformity is an important measure of consistency across image planes. From Figure 4.8 it can be seen that all algorithms outperform 3D F B P with respect to axial unifor-mity by about 20%. Unfiltered F O and O S E M have nearly the same uniformity, regardless of parameters, and filtering images produced with either method yields the same result. Axia l uniformity at low count rates is significantly lower but relative trends between algorithms are unchanged. 4.2.6 Radial Uniformity Radial uniformity at high counts is shown in Figure 4.9. When O S E M is resolution-matched to 3D F B P it is found to have a better radial uniformity. S A G E with log 2 (3 = \u00E2\u0080\u0094 5 has the same radial uniformity as 3D F B P with Hann 0.4 and when regularization is increased to log 2/3 = \u00E2\u0080\u00944 radial uniformity exceeds that of 3D F B P and is at par with FO(12,2,4). CHAPTER 4. EXPERIMENTAL RESULTS 2.4-2.2-1.80-g ^ 1.75-| .2 1.70 q 3 1.65 \u00E2\u0080\u0094 1.60-1.55 Axial Uniformity High Counts 3DFBP F0RE/0SEM 2DFBP FORE/SAGE: o o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 : A A o o A A I I I I I I I I I I I I I I I I Filtering (mm) \u00E2\u0080\u00A2 None \u00E2\u0080\u00A2 7,9,2 O 3,3,3 A 4,4,4 Figure 4.8: Axia l uniformity for various reconstructions at high count rates. CHAPTER 4. EXPERIMENTAL RESULTS Radial Uniformity High Ceunts l l l l l l l l l l l l l l ^^^^^^% * v V v \u00E2\u0080\u0094 4,x) possibilities because of inferior phantom study results. CHAPTER 4. EXPERIMENTAL RESULTS N E M A Uniformity High Ceunts Lew Ceunts I I I I I I I I I I I I I I I I 3DFBP 2DFBP 0SEM SAGE I I I I II I I I I I I I I I I 3DFBP 2DFBP GSEM SAGE Filtering (mm) \u00E2\u0080\u00A2 none \u00E2\u0080\u00A2 7,9,2 O 3,3,3 Figure 4.10: N E M A uniformity for various reconstructions in two count regimes. CHAPTER 4. EXPERIMENTAL RESULTS Average T A C Variability 1.75 1.70 1.65 -=| 1.60 1.55 q 1.50 1.45 -= 1.40 1 I I FORE/SAGE log2f3=-6 10 iter FORE/SAGE log2P=-5 10 iter 0.9 -s \u00E2\u0080\u00A2 r-t > 0.7 - \ 0.6 -J 0.5 16,4 ...^ F0RE/SAGE log2p=-4 10 iter 3DFBP Hann 0.5 16,1 F0RE/0SEM (subsets,iterations Figure ef merit avg p(stdev fT f t ) /avg fTj where R 0 I activity f frame 1...12 p plane 3...29 F0RE/2DFBP Ramp 0.5 F0RE/2DFBP Hann 0.5 F0RE/2DFBP Hann\"OA i 1 r 2 3 4 xyz filter (mm) Figure 4.11: Variability in the T A C for various reconstructions. CHAPTER 4. EXPERIMENTAL RESULTS 70 4.3.1 Radial Profile Difference One of the raclopride human subject scans was used to visually examine the radial profiles. The vertical profile running down the center of the image, parametrized by SR = {(*>.?) I * = 64} (recal that the image is 128 x 128), was extracted for several reconstruction methods. Figure 4.12 shows these profiles. To estimate the average difference in the radial Horizontal Profile of a Brain Scan F2D Hann 0.4 FO(12,2,0) F0(12,2,3) FO(12,4,4) FS(-6,10) FS(-4,10) 20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 Figure 4.12: Radial profiles taken vertically through a raclopride brain image. Several recon-structions are shown. The solid line represents 3D F B P and the scatter points are the algorithm indicated above the plot. CHAPTER 4. EXPERIMENTAL RESULTS 71 profiles, all 10 human subject scans were used and the same vertical profile as in Figure 4.12 was calculated and averaged over the central planes 10-20 for all tested reconstruction methods. The subset SR was further restricted to SR = | i = 64,25 < j < 100} so that values near zero would not be compared, since their relative difference would artificially increase A R . Figure 4.13 shows the average A R value for each reconstruction method. The correspondence between 3D F B P and F2D is excellent and most methods fall within 1% of 3D F B P , showing that there is no systematic quantitation offset. Radial Profile Difference i r i i i i r **V. <*Vv \"Av \"A* o, \u00C2\u00A9, a. a. \"a. \"s, ~<8-o t7-> 0-> Q-> < A < A r v . 'vs\ 'v . \u00E2\u0080\u00A2O) 'J/ '^i -O) 'ty '^i -O/ -Of -J} -Oj Figure 4.13: Average radial profile difference for various reconstruction methods. 4.3.2 Region-of-Interest Difference Figure 4.14 shows average, normalized ROI differences for 3 ROI groups. The differ-ences are slighly larger than found using radial profiles. Although both methods test quantita-tion, the radial profiles were averaged over the central planes (10-20) whereas the ROI differences were averaged over planes 3,5,10,15,20,25 and 28. Outer planes typically contain fewer counts CHAPTER 4. EXPERIMENTAL RESULTS and low statistics images are expected to yield larger differences between algorithms. 72 4.3.3 Image R M S Figure 4.15 shows image R M S and signed difference averaged over 10 human subject scans. The R M S is used to gauge the overall difference regardless of sign and the signed difference is used to test any systematic offset. F2D continues to prove itself as an excellent replacement for 3D F B P , with the R M S and signed difference being very small. S A G E shows a significantly larger signed difference than all the other algorithms showing that it reconstructs absolute activities smaller than 3D F B P . This large difference was not seen in the radial profiles but was seen in the ROI differences for the outer ROI group (but not central). This suggests that there is a larger difference between 3D F B P and FS as one moves radially outward. 4.3.4 Kinetic Parameter Comparison The kinetic parameters were calculated and compared as outlined in Section 3.3.4. Figures 4.16 and 4.17 show the average difference in Logan slopes and average dispersion in these slopes, respectively, for a group of 10 human subjects. The standard deviation of the D V R shown in Figure 4.16 represents variability between human subjects. The dispersion shown in Figure 4.17 represents the average variability for different ROIs for a single human subject. F2D yields the best results. For both cortical and cerebellar inputs, A s is consistent with zero and has a very small dispersion of < 0.5%. Most of the other algorithms give slope values greater than those found with 3D F B P , as would be expected for methods with a higher resolution. Both F O and FS perform consistently, but do show a systematic offset in the Logan slopes. CHAPTER 4. EXPERIMENTAL RESULTS 73 FO(12,4,3), for example, shows a difference A s c o r t = -2 .5 \u00C2\u00B1 1.8% and A s c e r e b = \u00E2\u0080\u00943.9 \u00C2\u00B1 1.1. Cortical input results for F O overlap, except for the (12,4,4) case which yields a lower slope than expected. Cerebellar input results are more scattered but show the same relative trends as cortical inputs. 4.3.5 Computational Demands A n important F O M is the amount of time required for reconstruction. Currently, each frame of a scan requires about 8 minutes for transmission processing and 20 minutes for 3D F B P . Rebinning with F O R E requires about 0.6 minutes per frame. F2D requires approximately the same amount of time per frame so the entire F2D reconstruction requires 1.2 minutes per frame. Both F O and FS are slighly slower because of the nature of the iterative reconstruction. For 4 iterations O S E M takes about 2.4 minutes per frame (per iteration it takes the same amount of time as 2D F B P ) . S A G E requires about 0.8 minutes per frame per iteration, so 5 iterations require 4 minutes. Figure 4.18 compares frame reconstruction times graphically assuming 8 minutes for transmission processing, 0.6 minutes for F O R E and 4 iterations for F O and 5 iterations for FS. The total reconstruction time for a 16 frame study, without counting transmission processing, is 320 minutes (3D F B P ) , 16 minutes (F2D), 48 minutes (FO(5,2,/)) and 64 min-utes (FS(log2(3,5)). If complete processing is included (normalization, scatter correction and attenuation) then 160 minutes should be added to each of these times. Another way to present the computational demands is to give the maximum number of frames that can be reconstructed in a 24 hour period using these algorithms, assuming 100% computer uptime. Wi th 3D F B P one can reconstruct 50 frames per day. F O R E radically CHAPTER 4. EXPERIMENTAL RESULTS 74 increases the potential number of frames that can be reconstructed: 155 with 2D F B P , 130 with O S E M and 115 with S A G E . \u00E2\u0080\u00A2S CHAPTER 4. EXPERIMENTAL RESULTS 75 Regional R 9 I Differences 4 -4 o -A -2 -4 Striatial ROIs \u00C2\u00A9 0 -2 10 Central R01s - Outer ROIs 6 4 - | 2 0 -2 i 1 1 1 1 ' i i1 1 1 1 1 1 A ^ ^ ^ A> A> Av A JuJ/ J<9/ -\"c?/ ^ ^ ^ ^ ^ Figure 4.14: Average normalized ROI differences for three ROI groups. CHAPTER 4. EXPERIMENTAL RESULTS 76 Image R M S and Signed Difference RMS 0.5 0.4 - j erf 0.2 0.1 0.0 4.5 o d \u00C2\u00AB St) Signed Difference i i i i i i i i r % \ \ \ % \ W W W <-> O Q fsp Q ? 0 \ A a \ <-e. Tjr 'Ig 'SQ 'SQ >S0 >S0 'OJ 'Jj Y'Oj *>y \">o, \"oj \"jj \"qj Figure 4.15: Image RMS and signed difference for various reconstructions averaged over 10 human subject scans. CHAPTER 4. EXPERIMENTAL RESULTS 77 CD O CI <*k 4 A> 4 4. A <9, S\u00E2\u0080\u009E. \u00C2\u00A9\u00E2\u0080\u009E. \u00C2\u00A9\u00E2\u0080\u009E. Q, t\u00C2\u00A3- ^ * ^ ^ iTp <7p r v *V Oj '\u00C2\u00A5j 'Oj 'Oj '<}j J Figure 4.16: Average difference in Logan slopes (DVR) between 3D F B P Hann 0.4 and various reconstruction methods. The error bars represent the standard deviation of the D V R among the 10 human subjects. CHAPTER 4. EXPERIMENTAL RESULTS 78 Legan Slepe Difference Dispersien Certical Input Function 6 -q Cerebellar Input Function i i i i i i i r A -^ r. ^ \"^ o <*!r> ^ ^ A> A> A> A> A S , S , S , S , 9 , S , < f < f < f * \u00C2\u00AB -s >Oj *&) J \u00C2\u00A5 j >Oj J & ) J \u00C2\u00A5 j >0) J ? j >0/ - / Figure 4.17: Average Logan slope dispersion of the difference between 3D F B P Hann 0.4 and various reconstruction methods. CHAPTER 4. EXPERIMENTAL RESULTS 79 I a a B Cemputatienal Demands ef Tested Algorithms 30 20 10 28 I | Reconstruction BB Processing (scatter, normalization, attenuation) 13 l u l l 3D FBP F2D F9 FS Recenstructien method Figure 4.18: Reconstruction times for standard 3D F B P and F O R E reconstructions. The time is significantly reduced because 2D reconstruction methods can be utilized after rebinning with F O R E . 80 Chapter 5 Discussion and Conclusion 5.1 Algorithm Performance When evaluating the performance of any algorithm one should identify the F O M s that are to be considered first, followed by secondary ones. For example, if speed is of the essence, a faster algorithm which has slighly lower resolution would be prefered to one which has perhaps higher resolution but is slower. On the other hand, if computing resources are not an issue, one must still weight the noise/resolution factors and choose what noise level or resolution is acceptable. If the decisive factor is noise then the parameters should be chosen such that the noise is the highest allowable to provide best resolution. A l l phantom F O M and F O M differences results are summarized the following tables: F O R E + 2D F B P in Table 5.1, F O R E + O S E M in Tables 5.2 and 5.3 and F O R E + S A G E in Tables 5.4 and 5.5. human subject scan F O M differences are given for all algorithms in Table 5.6 (quantitation) and Table 5.7 (DVR). CHAPTER 5. DISCUSSION AND CONCLUSION 81 Method F O M VR (mm 3) Noise C R N U M A X CV Uax Urad (%) Speed 3 D F B P 548 0.12 3.23 0.25 0.09 2.30 5.80 0.94 20 F2D 548 0.13 3.24 0.27 0.10 1.80 5.60 0.62 1 A ( % ) 0.0 -8 -0.3 -8 -11 22 3 34 95 Table 5.1: Comparison of Phantom F O M s for 3D F B P and F O R E + 2D F B P . Speed is given in minutes per frame. 5.1.1 F O R E and 2D F B P A l l of the phantom F O M s show that F2D either equals in performance to 3D F B P or supersedes it, as in the case of axial uniformity and T A C variability. Practically, it is the 95% speed increase (i.e. 20 times faster) over 3D F B P , with all other FOMs relatively unchanged, that makes the F O R E + 2D F B P choice so appealing. Tables 5.6 and 5.7 show that F2D quantitation is excellent, within 1% of 3D F B P and D V R differences are also very small. Importantly, the dispersion of D V R is < 0.5%, signifying that on average for a single human subject all of the ROI values used to calculate the average D V R fall very close to the 3D F B P values. 5.1.2 F O R E and O S E M Table 5.3 illustrates the need for a posteriori filtering for O S E M . Certainly the FO(12,2,0) and (12,4,0) options have very poor noise characteristics compared to 3D F B P . The N E M A uni-formity values suggest that these images are rougher than 3D F B P . Looking at the noise and radial uniformity FOMs, the only cases which are as good or better than 3D F B P are (12,2,3), (12,2,4) and (12,4,4). When using the wider 4 mm filter, the resolution is not not as good CHAPTER 5. DISCUSSION AND CONCLUSION 82 Method F O M VR (mm 3) Noise C R NUmax CV Uax UTad T ACVar (%) Speed 3 D F B P 548 0.12 3.23 0.25 0.09 2.30 5.80 0.94 20 FO(12,2,0) 120 0.14 3.40 0.27 0.093 1.80 5.2 0.80 2 FO(12,2,3) 220 0.10 3.34 0.20 0.077 1.71 4.9 0.76 2 FO(12,2,4) 310 0.087 3.29 0.18 0.068 1.67 4.7 0.68 2 FO(12,4,0) 120 0.22 3.36 0.32 0.12 1.81 6.9 0.81 3 FO(12,4,3) 220 0.14 3.31 0.25 0.094 1.72 6.3 0.68 3 FO(12,4,4) 310 0.11 3.27 0.21 0.083 1.67 6.0 0.61 3 Table 5.2: Comparison of Phantom FOMs for 3D F B P and F O R E + O S E M . A denotes the percent improvement in going from 3D F B P to F O R E and 2DFBP. Negative (positive) values of A signify a decrease (increase) in quality of the particular F O M . Speed is given in minutes per frame. Method A(%) VR Noise C R NUmax CV Uax Urad Speed FO(12,2,0) 78 -20 5 -6 -3 22 10 15 90 FO(12,2,3) 61 14 3 18 14 26 16 19 90 FO(12,2,4) 43 27 2 29 24 27 19 28 90 FO(12,4,0) 78 -83 4 -24 -29 21 -19 14 85 FO(12,4,3) 61 -17 2 0 -5 25 -10 27 85 FO(12,4,4) 43 0 1 28 8 27 -5 36 85 Table 5.3: Comparison of Difference in Phantom F O M s between 3D F B P and F O R E + O S E M . Negative (positive) values of A signify a decrease (increase) in quality of the particular F O M . but noise is reduced. For 4 iterations, O S E M image noise is the same as 3D F B P but volume resolution is improved by 43%. Radial uniformity, however, is reduced by 5% for the (12,4,4) case. Based on phantom F O M s it appears that 2 iterations is sufficient, with 4 iterations resulting in increased image noise and decreased uniformity. Performing only 2 iterations re-quires 2 minutes per frame whereas 4 iterations require 3 minutes per frame. This is still 10 CHAPTER 5. DISCUSSION AND CONCLUSION 83 Method F O M VR (mm 3) Noise C R N U m a x CV Uax Urad (%) Speed 3 D F B P 548 0.12 3.23 0.25 0.09 2.30 5.80 0.94 20 FS(-6,10,0) 180 0.16 3.33 0.30 0.11 1.77 6.70 1.68 7 FS(-6,10,4) 410 0.11 3.23 0.21 0.08 1.60 6.10 1.67 7 FS(-5,10,0) 195 0.13 3.24 0.27 0.10 1.78 5.90 1.52 7 FS(-5,10,3) 330 0.10 3.17 0.22 0.08 1.69 5.60 1.47 7 FS(-4,10,0) 205 0.11 3.07 0.24 0.09 1.77 4.90 1.50 6.50 Table 5.4: Comparison of Phantom F O M s for 3D F B P and F O R E + S A G E . Speed is given in minutes per frame. and 7 times faster than 3D F B P . Quantitation tests using human subject scans do not reveal any benefit in performing 4 iterations instead of 2. In fact, quantitation is a little better with only 2 iterations as seen in Table 5.6. D V R difference between 3D F B P and F O R E + O S E M is smaller using 4 iterations as seen in Table 5.7, but only for cerebellar input function. In fact, with FO(12,4,4) cerebellar input D V R difference is -0.6% whereas FO(12,2,4) yields a -3.7% difference. Generally, O S E M appears to overestimate the cerebellar input D V R yielding a negative difference. 5.1.3 F O R E and S A G E The results for S A G E are slighly less favourable than for O S E M . Based on noise characteristics, Table 5.5 shows that the log 2 /3 = \u00E2\u0080\u00944 case is the only one with noise lower than 3D F B P , if no filtering is used. In general, the filter can be seen to (a) decrease resolution (b) decrease noise and (c) increase uniformity. If the inclusion of a filter is acceptable, then log 2 j3 = \u00E2\u0080\u00945 with a 3 mm filter is nearly equivalent to log 2 /? = \u00E2\u0080\u00944 from the point of view of phantom FOMs. In filtering the log 2 0 = \u00E2\u0080\u0094 6 case with a 4 mm filter the resolution is increased CHAPTER 5. DISCUSSION AND CONCLUSION 84 Method A(%) VR Noise C R N U m a x CV Uax Urad T A C v a r Speed FS(-6,10,0) 67 -33 3 -20 -26 23 -16 -79 65 FS(-6,10,4) 25 13 0 17 9 30 -5 -78 65 FS(-5,10,0) 64 -9 0 -8 -13 23 -2 -62 65 FS(-5,10,3) 40 13 -2 12 6 27 3 -56 65 FS(-4,10,0) 63 11 -5 5 1 23 16 -60 65 Table 5.5: Comparison of Difference in Phantom FOMs between 3D F B P and F O R E + S A G E . Negative (positive) values of A signify a decrease (increase) in quality of the particular F O M . significantly and is only 25% higher than 3D F B P . Wi th none of the phantom F O M s showing dramatic improvement over O S E M , the fact that S A G E is significantly slower (requiring more iterations) makes this algorithm a less favourable replacement. human subject scan quantitation is good, but as already mentioned, outer ROI dif-ferences are large. Additionally, the signed image difference is large, caused most likely by differences in the images at larger radial distances. The D V R values calculated with FS are best for the log 2 (3 \u00E2\u0080\u0094 \u00E2\u0080\u00945 case with a 3 mm filter. As for O S E M , the D V R is systematically overestimated, except for log 2/3 = \u00E2\u0080\u00944 where it is 2.4% and 1.4% too small for cortical and cerebellar inputs. Table 5.4 shows that F O R E and 10 iterations of S A G E require require about 7 minutes. As mentioned, 5 iterations is sufficient to reach a stable image so this time can be expected to be closer to 4 minutes in any practical reconstructions. CHAPTER 5. DISCUSSION AND CONCLUSION 85 Method A(%) AR AROIs AROIc AROI0 R M S Ax F2D 0.04 1.1 1.3 1.2 0.05 -0.13 FO(12,2,0) 0.4 1.5 1.5 1.0 0.23 0.02 FO(12,2,3) -0.7 1.3 1.7 0.7 0.13 0.18 FO(12,2,4) -0.9 1.1 1.8 0.5 0.12 0.28 FO(12,4,0) 1.3 1.3 0.8 1.1 0.39 0.10 FO(12,4,3) 1.0 1.0 1.3 0.5 0.19 0.26 FO(12,4,4) 0.7 0.5 1.4 0.3 0.14 0.36 FS(-6,10,0) 1.3 1.7 1.9 4.3 0.19 3.6 FS(-6,10,4) 0.5 1.2 1.8 3.7 0.12 3.7 FS(-5,10,0) 0.7 2.8 2.4 4.8 0.13 4.1 FS(-5,10,3) 0.3 2.5 2.5 4.6 0.11 4.1 FS(-4,10,0) -0.9 1.0 0.8 3.4 0.12 3.7 Table 5.6: Comparison of difference in human subject study quantitation F O M s between 3D F B P and F O R E + 2D F B P / O S E M / S A G E . Method A(%) Ascor Dispersion A s c e r Dispersion F2D 0.2 0.4 -1.5 0.4 FO(12,2,0) -5.4 5.0 -8.9 5.6 FO(12,2,3) -2.5 4.2 -5.7 4.4 FO(12,2,4) -0.7 3.7 -3.7 3.9 FO(12,4,0) -2.5 4.4 -5.7 4.8 FO(12,4,3) 0.2 3.6 -2.3 3.9 FO(12,4,4) 1.8 3.2 -0.6 3.4 FS(-6,10,0) -6.9 2.8 -9.4 3.4 FS(-6,10,4) -2.3 2.3 -3.9 2.4 FS(-5,10,0) -2.5 2.2 -4.5 1.6 FS(-5,10,3) -0.1 2.1 -1.5 2.1 FS(-4,10,0) 2.4 2.0 1.4 1.9 Table 5.7: Comparison of difference in distribution volume ratios between 3D F B P and F O R E + 2D F B P / O S E M / S A G E . CHAPTER 5. DISCUSSION AND CONCLUSION 86 5.2 Algorithm Selection Many F O M s have been calculated for each algorithm to understand how it behaves under various artificial as well as practical conditions. A l l of the algorithms tested are faster than 3D F B P so if any one of these is chosen to replace 3D F B P , reconstruction time will be decreased significantly. A l l of the tests suggest that replacing 3D F B P with F O R E + 2D F B P does not change any of the image characteristics, nor any of the Logan analysis. The sole major difference between these two reconstruction methods is speed: F O R E + 2D F B P is 80% faster than 3D F B P , allowing our center to process 120 frames per day instread of only 48. Accuracy and reproducibility are important in clinical studies. Increased resolution may improve accuracy, but reproducibility may be compromised. On the other hand, if the image is greatly smoothed, accuracy can suffer while reproducibility is increased1. Both O S E M and S A G E can be tuned to yield equivalent phantom F O M s to those of 3D F B P with associated improvements in resolution. In fact, with O S E M volume resolution can be increased by as much as 61% (FO(12,2,3)) and 40% for S A G E (FS(-5,10,3) where D V R agreement with 3D F B P is best). O S E M quantitation is excellent, whereas S A G E quantitation is poorer. In addition, both algorithms appear to overestimate D V R values by 2-3% compared to 3D F B P . For studies where high resolution and minimum noise is desired F O R E + O S E M is a strong candidate, although the biological constants will be slighly different from the F2D results. Certainly, F O R E + 2D F B P can immediately replace 3D F B P , with an increase in overall reconstruction speed as well as the quality of some figures of merit, such as T A C variability 1If one smooths all every human subject scan with an extremely wide Gaussian, all the scans will look the same every time (excellent reproducibility) but none of the analysis will be correct (very poor accuracy). CHAPTER 5. DISCUSSION AND CONCLUSION 87 and image uniformity. On the other hand, it is important to realize that practically a 2% difference in D V R values will never be statistically significant in any P E T study of a biological system, owing to (a) large inherent variability in anatomy and physiology between human subjects (b) detector drift and normalization and (c) fluctuations in the data due to the statistical nature of radioactivy decay as well as the presence of scatters and randoms. Indeed, typically differences of <10% in D V R are not interpretable without very large populations of >50 subjects. 88 Bibliography [1] Matej S, Herman G T , Narayan T K , Furuie SS, Lewitt R M , Kinahan P E (1994) Evaluation of task-oriented performance of several fully 3D P E T reconstruction algorithms. Phys Med Biol 39: 355-367. [2] Kinahan P E , Matej S, Karp JS, Herman G T , Lewitt R M (1995) A comparison of transform and iterative reconstruction techniques for a volume-imaging P E T scanner with a large axial acceptance angle. I E E E Trans Nuc Sci 42: 2281-2287. [3] Defrise M , Geissbuhler A , Townsend D W (1994) A performance study of 3D reconstruction algorithms for positron emission tomography. Phys Med Biol 39: 305-320. [4] Matej S, Karp JS, Lewitt R M , Becher A J (1998) Performance of the Fourier rebinning algorithm for P E T with large acceptance angles. Phys med Biol 43: 787-795. [5] Reader A J , Visvikis D , Erlandsson K , Ott R J , Flower M A (1998) Intercomparison of four reconstruction techniques for positron volume imaging with rotating planar detectors. Phys Med Biol 43: 823-834. [6] Ter-Pogossian M M , Raichle M E , Sobel B E (1980) Positron-emission tomography. Scientific American 243: 169-181. [7] Phelps M E , Mazziotta J (1985) Positron emission tomography: human brain function and biochemistry. Science 228: 799-809. [8] Phelps M E (1991) P E T : A biological imaging technique. Neurochem Res 16: 929-940. [9] Phelps M E , Hoffman E J , Mullani N A , Ter-Pogossian M M (1975) Application of annihila-tion coincidence detection to transaxial reconstruction tomography. J Nuc Med 16:210-233. [10] Phelps M E , Hoffman E J , Huang S-C, Kuhl D (1978) E C A T : A new computerized to-mographic imaging system for positron-emitting radiopharmaceuticals. J Nucl Med 19: 635-647. [11] Edited by Phelps M E , Mazziotta J C , Schelbert HR. Positron Emission Tomography and Autoradiography: Principles and applications for the brain and heart. Raven Press, New York, 1986. [12] Perkins, D H . Introduction to High Energy Physics. 3rd edition. Addison-Wesley Publishing. 1987 [13] Rogers J G , Harrop R, Kinahan P E (1987) The theory of three-dimensional image recon-struction for P E T . I E E E Trans med Imag MI-6: 239-243. B I B L I O G R A P H Y 89 Bailey D L , Miller M P , Spinks T J , Bloomfield P M , Livieratos L , Young H E , Jones T (1998) Experience with fully 3D P E T and implications for future high-resolution 3D tomographs. Phys Med Biol 43: 777-786. Kaczmarz S (1937) Angenaherte auflosung von systemen linearer gleichungen. Bull Acad P o l o n a i s e Sci Lett Classe Sci Math Natur Series A , pp. 355-7. Herman G T , Meyer L B (1993) Algebraic reconstruction techniques can be made compu-tationally efficient. I E E E Trans Med Imag 12:600-9. Herman G T , Lent A , Lutz P H (1978) Relaxation methods for image reconstruction. Coram of A C M 21:152-8. Shepp L A , Vardi Y (1982) Maximum likelihood reconstruction for emission tomography. I E E E Trans Med Imag MI-2:113-122. Lange K , Carson R (1984) E M reconstruction algorithms for emission and transmission tomography. J Comp Assist Tomography 8:306-16. Dempster A P , Laird N M , Rubin D B (1977) Maximum likelihood from incomplete data via the E M algorithm. J Roy Stat Soc B 39:1-38. Liow J-S, Strother SC, Rehm K , Rottenberg D A (1997) Improved resolution for P E T volume imaging through three-dimensional iterative reconstruction. J Nuc Med 38:1623-31. Kaufman L (1987) Implementing and accelerating the E M algorithm for P E T . I E E E Trans Med Imag MI-6: 37-51. Llacer J , Andreae S, Veklerov E , Hoffman E J (1986) Towards a practical implementation of the M L E algorithm for positron emission tomography. I E E E Trans Nuc Sci 33: 468-477. Reader A J , Erlandsson K , Flower M A , Ott R J (1998) Fast accurate iterative reconstruction for low-statistics positron volume imaging. Phys Med Biol 43: 835-846. Daube-Witherspoon M E , Muehllehner G (1987) Treatment of axial data in three-dimensional P E T . J Nucl Med 28: 1717-1724. Sossi V , Stazyk M W , Kinahan P E , Ruth T J (1994) The performance of the single-slice rebinning technique for imaging thehuman striatum as evaluated by phantom studies. Phys Med Biol 39: 369-380. Erlandsson K , Esser P D , strand S-E, van Heertum R L (1994) 3D reconstruction for a multi-ring P E T scanner by single-slice rebinning and axial deconvolution. Phys Med Biol 39: 619-629. Lewitt R M , Muehllehner G, Karp JS (1994) Three-dimensional image reconstruction for P E T by multi-slice rebinning and axial image filtering. Phys Med Biol 39: 321-339. Defrise M , Kinahan P E , Townsend D W (1995) A new rebinning algorithm for 3D positron emission tomography: principle, implementation and performance. Proceedings of the 1995 International Meeting on Fully 3D Image R e c o n s t r u c t i o n in Radiology and Nuclear Medicine pp 235-240. Lange K , Bahn M , Little R (1987) A theoretical study of some maximum likelihood algo-rithms for emission and transmission tomography. I E E E Trans Med Imag MI-2 pp. 106-14. Lewitt R M , Muehllehner G (1986) Accelerated iterative reconstruction for positron emis-sion tomography based on the E M algorithm for maximum likelihood estimation. I E E E Trans Med Imag MI-5 pp. 16-22. BIBLIOGRAPHY 90 [32] Nuyts J , Suetens P, Mortelmans L (1993) Acceleration of maximum likelihood reconstruc-tion, using frequency amplification and attenuation compensation. I E E E Trans Med Imag 12: 643-652. [33] Fessler J A , Hero A O (1994) Space-alternating generalized expectation-maximization algo-rithm. I E E E Trans Sig Proc 42:2664-77. [34] Fesser J A , Hero A O (1995) Penalized maximum-likelihood image reconstruction using space-alternating generalized E M algorithms. I E E E Trans Image Proc 4:1417-29. [35] Hudson H M , Larkin RS (1994) Accelerated image reconstruction using ordered subsets of projection data. I E E E Trans Med Imag 13:601-9. [36] J A Fessler (1992) Hidden-data spaces for maximum-likelihood P E T reconstruction. 1992 I E E E Nuclear Science Symposium and Medical Imaging Conference. Conference Record. [37] Fessler J A , Clinthorne N H , Rogers W L (1993) On complete-data spaces for P E T recon-struction algorithms. I E E E Trans Nuc Sci 40:1055-61. [38] Fessler J A , Hero A O (1993) New complete-data spaces and faster algorithms for penalized-likelihood emission tomography. 1993 I E E E Nuclear Science Symposium and Medical Imag-ing Conference. Conference Record. [39] Fessler J A et al. (1996) Spatial resolution properties of penalized maximum-likelihood image reconstruction methods. I E E E Trans Image Proc 5:1346-58. [40] Fessler J A , Rogers W L (1994) Uniform quadratic penalities cause nonuniform image res-olution (and sometimes vice versa). In Conf Rec of the I E E E Nuc Sci Symp Med Image Conf vol. 4, pp 1915-1919. [41] Furuie SS, Herman G T , Narayan T K , Kinahan P E , Karp JS, Lewitt R M , Matej S (1994) A methodology for testing for statistically significant differences between fully 3D P E T reconstruction algorithms. Phys Med Biol 39: 341-354. [42] Brix G et al. (1997) Performance evaluation of a whole-body P E T scanner using the N E M A protocol. J Nucl Med 38:1614-1623. [43] Ehrin E et al (1985) Preparation of 1 1 C-labelled raclopride, a new potent dopamine receptor antagonist: preliminary P E T studies of cerebral dopamine receptors in the monkey. Int J Appl Radiat Isot 36: 269-273. [44] Farde L , Halldin C, Stone-Elander S, Sedvall G (1987) P E T analysis of human dopamine receptor subtypes using n C - S C H 23390 and 1 ^-raclopride. Psychopharm 92: 278-284. [45] Logan J , Fowler JS, Volkow N D , Wang G-J , Ding Y-S , Alexoff D L (1996) Distribution volume ratios without blood sampling from graphical analysis of P E T data. J Cereb Blood Flow Metab 16: 834-840. [46] Volkow N D , Fowler JS, Gatley SJ, Logan J , Want G-J , Ding Y-S , Dewey S (1996) P E T evaluation of the dopamine system of the human brain. J Nuc Med 37: 1242-1256. [47] Logan J et al (1990) Graphical analysis of reversible radiologand binding from time-activity measurements applied to [N-1 1C-methyl]-(-)-cocaine P E T studies in human subjects. J Cereb Blood Flow Metab 10:740-747. [48] Logan J et al (1994) Effects of blood flow on [ n C ] raclopride binding in the brain: model simulations and kinetic analysis of P E T data. J Cereb Blood Flow Metab 14: 995-1010. BIBLIOGRAPHY 91 [49] Ollinger J M , Goggin A S (1996) Maximum likelihood reconstruction in fully 3D P E T via the S A G E algorithm. 1996 I E E E Nuclear Science Symposium. Conference Record. [50] Kaufman L (1993) Maximum likelihood, least squares and penalized least squares for P E T . IEEE Trans Med Imag 12:200-214. [51] Fessler J A (1994) Penalized weighted least-squares image reconstruction for positron emis-sion tomography. IEEE Trans Med Imag 13:290-300. [52] Kinahan P E , Karp JS (1994) Figures of merit for comparing reconstruction algorithms with a volume-imaging P E T scanner. Phys Med Biol 39: 631-642. "@en . "Thesis/Dissertation"@en . "1998-11"@en . "10.14288/1.0085160"@en . "eng"@en . "Physics"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "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."@en . "Graduate"@en . "Evaluation of new iterative and rebinning reconstruction algorithms in fully three-dimensional positron emission tomography"@en . "Text"@en . "http://hdl.handle.net/2429/8150"@en .