UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Fourier regularization applied to spect reconstruction Togane, Dylan Dalmar 2000

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2000-0595.pdf [ 20.69MB ]
Metadata
JSON: 831-1.0085173.json
JSON-LD: 831-1.0085173-ld.json
RDF/XML (Pretty): 831-1.0085173-rdf.xml
RDF/JSON: 831-1.0085173-rdf.json
Turtle: 831-1.0085173-turtle.txt
N-Triples: 831-1.0085173-rdf-ntriples.txt
Original Record: 831-1.0085173-source.json
Full Text
831-1.0085173-fulltext.txt
Citation
831-1.0085173.ris

Full Text

FOURIER REGULARIZATION APPLIED TO SPECT RECONSTRUCTION By Dylan Dalmar Togane B. Sc. Concordia University, 1997 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF PHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA June 2000 © Dylan Dalmar Togane, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics and Astronomy The University of British Columbia 6224 Agricultural Road Vancouver, B.C., Canada V6T 1Z1 Date: _ i Abstract Fourier Regularized Emission Computed Tomography (FRECT), a reconstruction method for Emission Computed Tomography (ECT), is presented and compared with other re-construction methods. The compromise which characterizes ECT reconstruction is between image resolution and stability with respect to noise. FRECT reconstruction consists of minimizing an objective function that has a "goodness of fit" term which fits the data to a specified resolution and a "Fourier regularization" term which stabilizes the reconstruction process. Goodness of fit terms that can be coupled with the Fourier regularizer include the Least Squares term (FRLS) and the Maximum Likelihood term (FRML). Measures are taken to ensure that the fit and regularization terms act on separate areas of the Fourier domain. In this work FRECT is applied to Single Photon Emission Computed Tomography (SPECT) data. Analytically simulated data, Monte-Carlo simulations, phantom experi-ments and clinically acquired SPECT data are reconstructed. The compromise between resolution and stability is measured for Filtered Back-Projection (FBP), Ordered Subset Expectation Maximization (OSEM), and two variations of FRECT. The sensitivity with respect to experimental noise in the data is found empirically for all reconstruction methods, and analytically for those methods that allow such analysis. ii Thanks to Pierre Marechal, Dr. Anna Celler, Dr. Ronald Harrop, Dr. Dave Axen, Dr. Phil Gregory, Troy Farncombe, Stephan Blinder, Paul-Jean Letourneau, Shira Avni, Dane Doleman, Sulya Fenichel, Manda Harmon, Sara Davis, Julie Kaisla, Rita and Nick, Nicole and Michelle Pollard, Dax Xenis, Katrin Bowen, Johnathan Rider, Suzanne Klocke, Jean Maeght, Niklas Roeber, Lorenz Drig Parinias, Mike Wyszinski, Sebastien Gagnon, Tara Fortier, Jason Burke, Jan Seuntjens, Blake Walters, Lei Xing, Sanyu Kiruluta, Steve Garber, Fran Deitz, Carolyn Garber, Rachel Garber, Guy Trepanier, Mohamud Togane and Madeline Lefebvre. iii TABLE OF CONTENTS Abstract i i Table of Contents iv List of Tables v i i List of Figures v i i i 1 Introduction 1 1.1 Aim of This Work 1 1.2 Structure of the Thesis 2 1.3 Nuclear Medicine Imaging 3 1.3.1 Applications 4 1.4 Radioactive Decay and Interactions of Photons with Matter 4 1.5 The Anger Camera 8 1.5.1 History 8 1.5.2 The Collimator 10 1.5.3 The Scintillating Crystal 10 1.5.4 PMTs, Electronics and Computer : . 11 1.5.5 Detector Efficiency 12 1.6 Emission Computed Tomography Image Reconstruction 12 2 The Fourier Slice Theorem and Computed Tomography 15 2.1 The Fourier and Radon Transforms . . 15 2.2 The Fourier Slice Theorem 19 2.3 Data and Noise in the Fourier Domain 20 iv 2.4 The Modulation Transfer Function 22 2.5 The Radon Transform and the System Matrix 24 2.6 Collimator Blurring 25 2.7 Attenuation 26 3 S P E C T Reconstruction Methods 28 3.1 Filtered Back-Projection (FBP) 28 3.2 Reconstruction as an Inverse Problem 31 3.2.1 Least Squares and Weighted Least Squares 32 3.2.2 Maximum Likelihood 33 3.3 Iterative methods 34 3.3.1 ML-EM and OS-EM 34 4 Regularized S P E C T Reconstruction 37 4.1 Classical Regularizers 38 4.2 Fourier Regularization 38 5 Sensitivity to Noise 41 5.1 FBP 42 5.2 Iterative Methods 43 6 Reconstructions 45 6.1 Simulation of SPECT data 45 6.1.1 Analytic Simulations 45 6.1.2 Monte-Carlo Methods 46 6.2 Digital Phantoms 46 6.2.1 Analytic Circles 46 v J 6.2.2 Monte-Carlo MCAT Phantom 47 6.3 Reconstruction Evaluation Criteria 48 6.3.1 Difference from Original Object and Re-projection Misfit to Data 48 6.3.2 Noise Response 50 6.4 Reconstruction Parameters 51 6.4.1 Selecting Cutoff Frequency and Filter 51 6.4.2 Regularization Parameters and Stopping Criteria 51 6.4.3 Parameter Ranges 52 7 Results and Discussion 54 7.1 Fit to Original Object and Fit to Data . . . 56 7.2 Response to Noise 71 7.3 Summary 89 8 Conclusion , 91 8.1 General Comments 91 8.1.1 Future Work 92 Bibliography 94 9 Appendix A 97 vi LIST OF T A B L E S 1.1 Nuclear Medicine Radio-pharmaceuticals and Applications - part 1 . . . 13 1.2 Nuclear Medicine Radio-pharmaceuticals and Applications - part 2 . . . 14 1.3 Types of Radioactive Decay 14 5.1 Gradients and Hessians of £ and p 44 6.1 Reconstruction Parameter Settings Used in Tests of Reconstruction Methods 53 7.1 FBP Results 54 7.2 Two-Parameter Results 55 7.3 Characteristics of "Norm(V) vs. / c " Plots 72 7.4 Summary of Results - Part 1 89 7.5 Summary of Results - Part 2 90 vii LIST OF FIGURES 1.1 Photo-electric Effect 6 1.2 Compton Scattering 6 1.3 Pair Production 7 1.4 Schematic Drawing of the Anger Camera's Basic Components 9 2.1 Geometry of the Radon Transform . 17 2.2 Radon Transformed Square and Profile 18 2.3 Frequency Domain Sampling in Emission Tomography 20 2.4 Digital Circle Phantom 20 2.5 Radial Fourier Transform of Sinogram and Noise . 21 2.6 Modulation Transfer Function Test , . 23 2.7 High Resolution Low Energy Modulation Transfer Function 23 2.8 Sparsity Structure of the Radon Matrix 25 2.9 Collimator Blurring Geometry 26 2.10 Attenuated Radon Transform Geometry 27 3.1 Back-Projection Operator 29 3.2 Sine Function . 31 3.3 ML-EM and OS-EM 35 3.4 Increasing Number of Iterations of OS-EM 36 4.1 Partitioning the Fourier Domain 39 6.1 Digital Circle Phantom 47 6.2 One Slice of the MCAT Phantom's Activity Map 48 viii 7.1 SelectedFBP Reconstructions (top), RMS Difference Between FBP Re-constructions and Analytic Circle Phantom vs fc (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and Analytic Circle Sinogram vs fc (bottom, right) 57 7.2 Selected FBP Reconstructions (top), RMS Difference Between FBP Re-constructions and MCAT Digital Phantom vs / c (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and MCAT Monte-Carlo Sinogram vs fc (bottom, right) 58 7.3 Least RMS Difference Between 4 Subset OSEM Reconstructions and Dig-ital Circle Phantom vs fc (left) and Corresponding Iteration No. vs fc (right) 59 7.4 Least RMS Difference Between FR-LS Reconstructions and Digital Circle Phantom vs fc (left) and Corresponding a vs fc (right) 60 7.5 Least RMS Difference Between FR-ML Reconstructions and Digital Circle Phantom vs fc (left) and Corresponding a vs fc (right) 61 7.6 Least RMS Difference Between 4 Subset OSEM Reconstructions and Dig-ital MCAT Phantom vs fc (left) and Corresponding Iteration No. vs fc (right) . 62 7.7 Least RMS Difference Between FR-LS Reconstructions and Digital MCAT Phantom vs fc (left) and Corresponding a vs fc (right) 63 7.8 Least RMS Difference Between FR-ML Reconstructions and Digital MCAT Phantom vs fc (left) and Corresponding a vs fc (right) 64 7.9 Least RMS Difference Between Re-projection of 4 Subset OSEM Recon-structions and Sinogram Data vs fc (left) and Corresponding Iteration No. vs fc (right) 65 ix 7.10 Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) . . . . 66 7.11 Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) . . . . 67 7.12 Least RMS Difference Between Re-projection of 4 Subset OSEM Recon-structions and Sinogram Data -vs fc (left) and Corresponding Iteration No. vs fe (right) 68 7.13 Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) . . . . 69 7.14 Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) . . . . 70 7.15 Noise Response of FBP Reconstructions vs. fc for Analytic Circle Data . 75 7.16 Noise Response of FBP Reconstructions vs. fc for MCAT Data 76 7.17 Noise Response of 4-Subset OSEM Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object 77 7.18 Noise Response of 4-Subset OSEM Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data . . . 78 7.19 Noise Response of FR-LS Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object . . . 79 7.20 Noise Response of RC-LS Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data 80 7.21 Noise Response of FR-ML Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object . . . 81 x 7.22 Noise Response of RC-ML Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data 82 7.23 Noise Response of 4-Subset OSEM Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object . . . 83 7.24 Noise Response of 4-Subset OSEM Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Expected Misfit to Data 84 7.25 Noise Response of FR-LS Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object 85 7.26 Noise Response of RC-LS Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Expected Misfit to Data 86 7.27 Noise Response of FR-ML Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object 87 7.28 Noise Response of RC-ML Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Expected Misfit to Data 88 9.1 Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital Circle Phantom vs fc and Iteration No 98 9.2 Surface Plot of RMS Difference Between FR-LS Reconstructions and Dig-ital Circle Phantom vs fc and a 99 9.3 Surface Plot of RMS Difference Between FR-ML Reconstructions and Dig-ital Circle Phantom vs fc and a 100 9.4 Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs fc and Iteration No 101 9.5 Surface Plot of RMS Difference Between FR-LS Reconstructions and Dig-ital MCAT Phantom vs fc and a 102 xi 9.6 Surface Plot of RMS Difference Between FR-ML Reconstructions and Dig-ital MCAT Phantom vs fc and a 103 9.7 Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs fc and a 104 9.8 Surface Plot of RMS Difference Between Reprojection of FR-LS Recon-structions and Sinogram Data vs fc and a 105 9.9 Surface Plot of RMS Difference Between Reprojection of FR-ML Recon-structions and Sinogram Data vs fc and a 106 9.10 Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs fc and a 107 9.11 Surface Plot of RMS Difference Between Reprojection of FR-LS Recon-structions and Sinogram Data vs fc and a . . : 108 9.12 Surface Plot of RMS Difference Between Reprojection of FR-ML Recon-structions and Sinogram Data vs fc and a 109 xii CHAPTER 1 I N T R O D U C T I O N 1.1 A I M OF THIS W O R K Single Photon Emission Computed Tomography (SPECT) is a diagnostic medical imag-ing technique based on the detection of gamma radiation emitted by tracers which are introduced into the body of the patient. 3D images' of the activity distribution within the patient's body are reconstructed from sets of 2D projections of the detected gamma radiation. An important challenge in the reconstruction process is the separation of the meaningful information from the unavoidable statistical noise in the data. In this work a new reconstruction method (Fourier Regularized Emission Computed Tomography) for SPECT is presented and compared with reconstruction methods in current use. The Fourier Slice Theorem, which relates the Radon and Fourier transforms, plays a central role in the theory behind Fourier regularization. For this reason, FRECT is well suited to tomographic imaging systems whose action can be described in terms of the Radon transform. This work can be described as an attempt to define a pseudo-inverse of the operator representing the SPECT acquisition system. This inversion should be done in a way that will not excessively amplify inaccuracies or inconsistancies due to noise in the measured 1 Chapter 1. Introduction 2 data. The objectives of the work are (i) to present some of the theory behind FRECT and (ii) to present results characterizing the performance of FRECT and other SPECT reconstruction methods on SPECT data, allowing comparison. 1.2 STRUCTURE OF THE THESIS This thesis has nine chapters. The first chapter contains background information about Nuclear Medicine imaging physics and detectors, and a description of how Fourier regu-larization is intended to complement the existing techniques. Chapter 2 presents the Fourier Slice Theorem and its relevance to ECT imaging data. Characteristics of experimental noise and data in the Fourier domain are also discussed. The variations on the Radon transform which implement attenuation correction, colli-mator blurring correction, and scatter correction are outlined with varying degrees of detail. A variety of formulations of the image reconstruction problem are presented in Chap-ter 3. The principles behind some existing reconstruction methods (FBP, ML-EM, OS-EM) are described. The general formulation of image reconstruction as an inverse prob-lem is also presented here. In Chapter 4, Fourier regularization for ECT reconstruction is presented, as well as three classical regularization terms. The Least Squares and Maximum Likelihood fit func-tions, as well as classical regularizers are discussed and their application demonstrated. Methods used to calculate each reconstruction method's sensitivity to experimental noise in the ECT data are presented in Chapter 5. Linear reconstruction methods are identified and special features of these methods are discussed. Reconstruction procedures for simulated and experimental data sets are detailed in Chapter 1. Introduction 3 Chapter 6. Also presented are the methods used to select values for the reconstruction parameters. Finally, the post-reconstruction calculations that will be used to characterize the reconstructions are described. Results of the calculations described in Chapter 6 are presented and discussed in Chapter 7. In conclusion, Chapter 8 summarizes the work and results. Issues related to clinical application and future research a^re also discussed. Additional results are presented without discussion in Chapter 9. 1.3 NUCLEAR MEDICINE IMAGING Diagnostic medical imaging began in 1895 when Wilhelm Roentgen produced X-Ray im-ages of the bones in his hand. Medical imaging now includes a variety of non-invasive imaging techniques including X-ray computed tomography (CT, pinoeered by Hounsfield [1]), planar nuclear medicine imaging (NM), positron emission tomography (PET), sin-gle photon emission computed tomography (SPECT), ultrasound (US), magneto- and electro-encephalography (MEG and EEG), and magnetic resonance imaging (MRI). These methods are used to provide images of the anatomy and physiology of a patient allowing for distinction between normal and diseased states. Techniques for non-invasive imaging have proven to be extremely useful tools in medicine, and hospitals world-wide are now equipped with medical imaging systems. NM, PET and SPECT are all based on the detection of radiation emitted from within the patient's body, and all fall into the category of Nuclear Medicine (unlike CT which uses an external radiation source). These techniques use radio-tracers, drugs labelled with a radioactive isotope which enter the patient's body by injection, inhalation or ingestion. Tracers are taken up by different systems or organs in the body and their Chapter 1. Introduction 4 detected distribution provides information about the target organ's function. They are used to study and detect some of the most common diseases in our society such as heart diseases, cancer, cerebro-vascular disorders and many others. 1.3.1 APPLICATIONS Stress-rest scans of the heart are a prime example of the application of SPECT imaging. A 2 0 1T1 or 9 9 m T c labelled drug, sestamibi is used. Sestamibi is a lipophilic, positively charged complex that is taken up by the cells of the heart during their normal function (necrotic or dead cells do not attract the drug). The quantity of drug taken in by the cells is proportional to the amount of blood flow through them. As the radio-isotope attached to the drug decays, the emitted radiation is detected making it possible to see where and in what quantity the drug is present. This type of study can indicate areas of the heart which are no longer functioning or no longer adequately supplied with blood. Besides the detection of areas of ischemia and infarction in the heart muscle, other areas where SPECT imaging is useful include the detection of areas of epilepsy and stroke in the brain, and the detection of tumor locations in the liver, kidney, thyroid and throughout the body. Table 1.1 shows some common isotopes, and the organs and medical applications for which they are used.[2] 1.4 RADIOACTIVE DECAY AND INTERACTIONS OF PHOTONS WITH M A T T E R Radioactive decay is a process in which an unstable [parent] nucleus trans-forms into a more stable [daughter] nucleus by emitting particles and/or pho-tons and releasing nuclear energy. Atomic electrons may become involved in Chapter 1. Introduction 5 some types of radioactive decay, but it is basically a nuclear process caused by nuclear instability. [3] Various processes by which a nucleus can decay into a more stable state are summa-rized in Table 1.3. In addition to the processes listed in Table 1.3, Internal Conversion, (IC) and Nuclear Fission are described here for completeness. Internal conversion occurs when the daughter nucleus is in a meta-stable state, and the photon emitted when the meta-stable state decays is absorbed by an orbital electron (called a conversion electron). A heavy nucleus undergoes fission when it spontaneously splits into two lighter nuclei and two or three fission neutrons. The decay processes that are most used in Nuclear Medicine are 7 emission and positron (/3+) emission. The gamma-ray in 7 emission is the photon referred to in "Single Photon Emission Computed Tomography", and positron emission is the type of radioac-tive decay that PET imaging is based on. In the case of PET, the emitted positron travels a very short distance (usually a few mm) before meeting an electron and annihilating, which produces two 511 KeV photons that are detected by the PET imaging system. The operation of the SPECT detection system is described in the next section. Once produced by nuclear decay, the emitted photons must travel some distance from the decay site within the patient to the detection system, outside the patient. During this journey the photon may interact with the medium through which it travels. The four types of interactions that can occur between atoms and photons having energies of slightly greater than 1 MeV or below are (i) the photo-electric effect, (ii) Compton scattering, (iii) pair production, and (iv) coherent (Rayleigh) scattering. In the photo-electric effect, the photon is absorbed by an electron orbiting the atom. The photon's energy is converted into kinetic energy causing the electron to be ejected from the atom. The electron's kinetic energy is equal to that of the incident photon Chapter 1. Introduction 6 Figure 1.1: Photo-electric Effect (illustration from [3]) minus the binding energy of the electron in the atom. Figure 1.2: Compton Scattering (illustration from [3]) Compton scattering occurs between photons an electrons in the loosely bound outer shells. The loosely bound electron absorbs the photon and is ejected from the atom. A photon of lower energy than the original is emitted in some new direction. The energy of the scattered photon depends on the direction in which it is emitted as well as the energy of the incident photon. When a photon having 1.022 MeV or greater energy interacts with the electric field of a charged particle, there is a chance that it will disappear and be replaced by a positron-electron pair. This process is called pair production and can be thought of as Chapter 1. Introduction 7 0 . 5 l l - M e V ANNIHILATION PHOTONS Figure 1.3: Pair Production (illustration from [3]) the reverse of positron-electron annihilation. The positron will typically travel a short distance before meeting an electron and annihilating. Coherent or Rayleigh scattering occurs when a low energy photon (less than 50 KeV) interacts with the electric field of an atom as a whole. In this process, very little recoil energy is absorbed by the atom because of it's great mass. For this reason, the photon is deflected with almost no loss of energy. When travelling through a medium, the chances that the photon will undergo any of the interactions described above depends on the energy of the photon and the composition of the medium. These interactions are collectively described by a material's attenuation coefficient (//). The differential change in intensity (/) of a photon beam passing through a thickness dx of matter is [3] dl/I = -pdx. (1.1) In order to find the radiation intensity after having passed through an inhomogeneous Chapter 1. Introduction 8 medium (p, = /i(x)), Equation E l is changed into an integral formula I = IQe-* "M*1', (1.2) where I0 is the initial beam intensity, and the integration is carried out over the path that the radiation travels. [3] 1.5 T H E ANGER CAMERA The Anger camera is the detection system used in SPECT. Because the operation of the Anger camera determines the nature of SPECT data, all reconstruction methods incorporate knowledge of Anger camera's operation in some way. Further, the history of the Anger camera is closely related to the history of SPECT. Figure 1.4 shows the basic components of the Anger camera. 1.5.1 HISTORY In the 1940's a single detector was positioned at various locations around the head of a patient. Only crude spatial information, such as left-right symmetry, was obtained. The Rectilinear scanner was introduced by Ben Classen in the 1950's. A mechanically controlled detector scanned a raster-like pattern over the area of interest. This technique suffered from long imaging times due to the sequential nature of the scanning. Hal Anger used a pin-hole in a lead block to project a gamma-radiation image of the source distribution in 1953. The image was projected onto a scintillating screen with photographic film behind it. Extremely long exposure times were necessary because of inefficiencies in the system (most losses occurred in the film). Therapeutic doses were required to produce interpretable images. The screen and film were later replaced with a single Nal crystal and Photo-Multiplier Chapter 1. Introduction 9 X - P O S I T I O N S I G N A L Y - P O S I T I O N S I G N A L C A T H O D E R A Y T U B E G A T I N G C I R C U I T R Y P U L S E - H E I G H T A N A L Y Z E R Z - P U L S E P O S I T I O N L O G I C C I R C U I T S Q P M - T U B E A R R A Y N o I ( T J 0 C R Y S T A L C O L L I M A T O R P A T I E N T Figure 1.4: Schematic Drawing of the Anger Camera's Basic Components (illustration from [3]) Chapter 1. Introduction 10 Tube (PMT) array. This formed the basis for the "Anger Camera" which is now the standard clinical nuclear imaging device. In 1963, Kuhl and Edwards presented the first tomographic images using the Anger camera. 1.5.2 T H E COLLIMATOR The lead collimator is necessary for source localization. Its function is to absorb photons that arrive outside a defined acceptance angle with respect to the camera face. Various collimator designs exist, including parallel hole, fan beam, diverging, pin-hole, and others. A typical parallel hole collimator is basically a thick sheet of lead with a large number of hexagonal holes. A typical collimator is 2.4-3.2 cm thick, has 1.4-2.2 mm wide holes, and 0.2-0.4 mm thick septa. The acceptance angle is usually on the order of three degrees. The collimator greatly reduces the number of photons that reach the scintillating crystal. The collimator also scatters some of the photons, but not enough to appreciably increase the background noise from which the useful data must be extracted. Another source of experimental error associated with the collimator is septal penetration, which occurs when a photon arriving outside the acceptance angle is not absorbed by the collimator's septa. 1.5.3 T H E SCINTILLATING CRYSTAL The scintillating crystal is used to convert gamma rays into scintillation photons that are detected by the PMT array. It is a single Nal crystal with a small amount of Thallium (Tl) added to improve performance at room temperature. One of the reasons Nal is used is its low production cost. Chapter 1. Introduction 11 A typical Anger camera uses a 40x50x0.9 cm crystal. The crystal thickness is chosen according to the energy of the radiation to be detected. Thinner crystals offer better spatial resolution, but are inefficient detectors of high energy radiation due to the higher probability that photons will pass through them without scintillating. The thickness given above is appropriate for low energy (140 KeV) radiation. 1.5.4 PMTs, ELECTRONICS AND C O M P U T E R The PMT array converts scintillation photons into electric pulses and amplifies them for further processing. A typical Anger camera will have 50 to 100 PMTs. A single scintillation event may be detected by many PMTs, and the detected intensities will vary depending on the PMT positions relative to that of the scintillation. The position logic electronics calculate the location of scintillation's centroid based on knowledge of the PMT locations and the relative signal strengths. The intensity of the scintillation is related to the energy of the detected gamma photon, allowing the electronic pulse height analyzer to filter out events that do not fall within a specified energy range. The principal use of this feature is to reject photons that have Compton scattered before reaching the detector. The position logic circuits will not function properly if more than one scintillation event is detected simultaneously. The detected pulse height will be the sum of the pulses produced by the two scintillations (pile up effect), resulting in erroneous energy calculations. The over-all result is that the Anger camera, like all radiation detectors, has a period of dead time after each detected event during which electronic circuits prevent new events from being accepted. This time is on the order of milliseconds. The role of the computer is to control the mechanics that move the camera head (an extremely massive device), store the acquired data, and reconstruct images of the source Chapter 1. Introduction 12 distribution. 1.5.5 D E T E C T O R EFFICIENCY The number of photons detected by an Anger camera is typically four orders of magnitude lower than the number of emitted photons. A number of factors contribute to this loss of sensitivity, but the dominant cause is the lead collimator. [3] In the standard experiment for measuring the Anger camera's efficiency as a radiation detector, a petri dish ten centimeters in diameter is placed ten centimeters from the detector face. Counts are collected until a measurement of the desired precision has been acquired, and the ratio of the detected counts to the average counts emitted in petri dish is calculated. Some typical results are 170°'^' (where c.p.m. is counts per minute), which is approximately equal to 5 cf'^'. [4] 1.6 EMISSION COMPUTED TOMOGRAPHY IMAGE RECONSTRUCTION The basic problem in ECT image reconstruction is to deduce what source distribution within the patient produced the measured distribution of radiation outside the patient. The main approach is to think of the acquisition system as a mathematical transformation from the source distribution to the acquired data, and then to attempt to invert the transformation. The first step is to find a mathematical representation of the acquisition system. Having defined this mathematical operator, there are a number of ways to invert it, each with different characteristics. Because the acquired data is always corrupted by experimental noise, it is important to quantify the effect that this noise will have on the reconstructed images. Chapter 1. Introduction 13 Organ Radio-pharmaceutical Diagnosis Brain and Spinal Cord 9 9 mTc-Sodium Pertecnetate ( 9 9 mTc04 ) 9 9 m T c . D T p A blood flow, tumor, brain death 123I-Isopropyl-p-iodoamphetamine (IMP) cerebral perfusion and metabolism 99mTc-Hexamethylpropylene Amine Oxime (HMPAO) cerebral perfusion (stroke, infarction, tumor) 9 9 mTc-Ethyl Cyteinate Dimer (ECD) cerebral perfusion lsF-Flouordeoxyglucose (FDG) glucose metabolism oxygen consumption 1 5 0-H 2 0 , 1 3 &Xe blood flow n i In-DTPA leakage of cerebro-spinal fluid (CSF), hydrocephalus Thyroid 131I-Sodium Iodide (and other iodine isotopes) thyroid function and structure 9 9 mTc-Sodium Pertechnetate morphology Lung 99mTc-Macroaggregated Albumin ( 9 9 m Tc-MAA) perfusion defects i 3 3 Xe and 1 2 YXe gas airway obstructions 99mTc-Aerosol abnormal distribution in bronchial diseases Liver 1311-Rose Bengal liver disease, biliary obstruction 9 9 m Tc-IDA derivatives (-HIDA, -PIPIDA, -BIDA, -DISIDA) gallbladder and biliary tree imaging 9 9 mTc-Sulfur Colloid and 9 9 mTc-Albumin Colloid tumor detection e7Ga-gallium abscesses and tumors Spleen 9 9 mTc-Sulfur Colloid and 9 9 mTc-Albumin Colloid morphology 99mTc-labeled Red Blood Cells morphology Pancreas 75Se-Selenomethionine digestive enzymes Kidney 1 2 3I- or 131I-Orthoiodohippurate (Hippuran) filtration rate and effective renal plasma flow (EPRF) 99mTc-Mercaptoacetylglycylglycine (MAG3) 9 9 m Tc-DTPA 99mTc-Gluceptate 9 9 m Tc-DMSA Table 1.1: Nuclear Medicine Radio-pharmaceuticals and Applications - part 1 [2] Chapter 1. Introduction 14 Organ (cont'd) Radio-pharmaceutical (cont'd) Diagnosis (cont'd) Skeleton 99mTc-Phosphate and Phosphonate complexes ( 9 9 m Tc-MDP and 9 9 m Tc-HDP) bone blood flow, and bone formation rate Heart 201T1-Thallous Chloride perfusion studies yymTc-Sestamibi yymTc-Teboroxime S2Rb-Rubidium Chloride PET perfusion studies 13N-Ammonia 18F-Fluorodeoxyglucose (FDG) PET metabolic studies 99mTc-Pyrophosphate myocardial infarct U 1 l n - and.131I-labeled Antimyosin Antibody y y m T c - R B C gated equilibrium cardiac blood pool Table 1.2: Nuclear Medicine Radio-pharmaceuticals and Applications - part 2 [2] Name Nuclear Notation Description /3~ Emission and (/3~,7) Emission Ay !Ll A V 2 A ->• z + l r Ay X v* X A v The e~ is called a beta particle (/?"). The daughter nucleus is in an excited state and then emits a photon. Positron (/?+) and (/3+,7) decay A Y ?1 A v z A Z-lY AZX £ i^Y* 4 Z_XY The daughter nucleus is sometimes in an excited state [(P+,j) decay]. Electron Capture (EC) and (EC,7) decay ix E4 i_xv iX E4 i^Y* 4 i_xY An orbital e~~ is captured by the nucleus, transforming a proton into a neutron. The daughter nucleus is sometimes in an excited state [(EC/y) decay]. Decay by a emission The a-particle consists of two neutrons and two protons (a ^He nucleus). Table 1.3: Types of Radioactive Decay [3] C H A P T E R 2 T H E F O U R I E R SLICE T H E O R E M A N D C O M P U T E D T O M O G R A P H Y This chapter will provide the mathematical context for emission computed tomography data and reconstruction. The Fourier transform is fundamental to the idea of spatial resolution, and the Radon transform provides a mathematical model of the way the Anger camera transforms the observed source distribution into the acquired data. 2.1 T H E FOURIER AND RADON TRANSFORMS The Fourier transform of a function /(r) is defined as /inf /(r)e-i27r<k'r>dr, (2.1) •inf where (•, •) denotes a standard scalar product, / returns a scalar value as a function of the spatial variable r, F (the Fourier transform of /) returns a scalar quantity as a function of the wave number k which has units of and will sometimes be referred to as spatial frequency. Spatial frequency can be loosely thought of as the inverse of spatial 15 Chapter 2. The Fourier Slice Theorem and Computed Tomography 16 resolution in the sense that if F contains high spatial frequencies, then / describes an object with fine resolution. The inverse Fourier transform is /inf F(k)e i 2 7 r < k ' r )dk. (2.2) -inf These two equations can be interpreted as a mapping from the spatial domain to the frequency domain and back. They apply to one, two or three-dimensional distributions. The Radon transform is ubiquitous - fields where it has been applied include Astron-omy, Molecular Biology, Geophysics, Material Science, Optics, and many others. Until the early 1970s, many investigators outside the field of Mathematics were not aware of Radon's original 1917 paper ([5]). From 1972 to 1973, it was pointed out that Radon's work was fundamental to the problem of reconstruction from projections ([6, 7, 8, 9]). Budinger, Gullberg and Huesman, and others did some of the earliest work in applying the Radon transform to Emission Computed Tomography [10, 11, 12, 13]. [14] The Radon transform of the two-dimensional (£ = (x, y)) function / ( £ ) is defined as /inf /(0*(p-<*,0)de, (2.3) •inf where 0 is a unit vector in the direction of 9, p (sometimes referred to as the radial variable) is the smallest distance between line of integration (L in Figure 2.1) and the origin and (•, •) denotes a standard scalar product. The function f(9,p) returns the line integral through / ( £ ) along a line at the angle 9s- and passing at a distance p from the origin (L in Figure 2.1). The Anger camera is designed to image a radiation distribution. Because the radiation being imaged passes through it's medium, each point on the detector face will detect the sum of the activity directly in front of it. The camera is rotated around the source distribution, so the description of activity being imaged is transformed from a function of Chapter 2. The Fourier Slice Theorem and Computed Tomography 17 Figure 2.1: Geometry of the Radon Transform Euclidean co-ordinates to a function of the camera's angle and the position on the camera face. The Radon transform can be used as an approximate mathematical description of the Anger camera. Figure 2.2 shows an example of the two dimensional Radon transform. In this case, the object being transformed is a square - it could be an enlarged view of a single image pixel. The form shown on the right is the Radon transform of the square. Below this is a profile taken at 9 — OATT. The profile is clearly divided into four regions, two of which are constantly increasing and decreasing slopes, corresponding to the values of p which result in the line of integration passing through two sides of the square that are at a right angle to each other. In the uniform, non-zero region the line of integration passes through two parallel sides of the square. The three dimensional detector geometry that will be studied in this work does not correspond to the geometry of the three dimensional Radon transform. Instead, a stack of two-dimensional Radon transforms is appropriate. What could be referred to as the cylindrical Radon transform is obtained by adding dependence on the z variable to both Chapter 2. The Fourier Slice Theorem and Computed Tomography 18 Figure 2.2: Radon Transformed Square and Profile / and / . In this case, f(9,p,z) is the set of two-dimensional projections of the function / at all angles 9. In SPECT, the two-dimensional Radon transform of an object is sometimes referred to as a sinogram. A general mathematical definition of the Radon transform is: Briefly, if / is a function defined on n-dimensional Euclidean space R", the Radon transform of / is determined by integrating / over all hyper-planes of dimension n — 1. Thus, if / is defined on the plane K 2 , the / [(the Radon transform of /)] is determined by the line integrals of / . And if / is defined on R 3 then / is determined by the surface integrals of / over two-dimensional planes. [14] Chapter 2. The Fourier Slice Theorem and Computed Tomography 19 2.2 T H E FOURIER SLICE T H E O R E M The Fourier Slice Theorem relates the Radon and Fourier transforms. It allows us to see the acquired data in the Fourier domain. With the Fourier Slice Theorem, the measurement process can be thought of as sampling the Fourier transform of the imaged object. This view of the problem is the basis for an understanding of Filtered Back-Projection and Fourier Regularization. In two dimensional space, the Fourier Slice Theorem states that the two dimensional Fourier transform of / is equal to the Radon transform of / followed by the Fourier transform of f(0,p) with respect to the radial variable p /inf mP)e-i2^dp, (2.4) •inf where kx = kpcos(9) and ky — kpsin(9). In case of the cylindrical Radon transform, the Fourier transform is taken with respect to both p and z, and everything else is unchanged. In tomography, f{9,p) is known (measured) for a finite number of 0j's. Equivalently, the Fourier transform of the acquired data F(9i,kp) is known on a star-shaped domain such as represented in Fig. 2.3. Further, the two-dimensional Radon transform is only sampled for discrete values of p - in the Fourier domain the consequence is that only a finite number of points on each spoke are actually measured. The measurement domain represented in Fig. 2.3 can be regarded as the experimental frequency coverage, as defined in [15]. It is known (cf. [15], again) that interpolating and extrapolating the Fourier transform of a (compactly supported) function is an ill-posed problem, meaning that there is simply not enough information present for the problem to be solved directly - some form of a priori information or constraints on the solution are necessary. Chapter 2. The Fourier Slice Theorem and Computed Tomography 20 63 Figure 2.3: Frequency Domain Sampling in Emission Tomography 2.3 DATA AND NOISE IN THE FOURIER DOMAIN Having seen the measurement process in the Fourier domain, the next step is to look at how experimental noise entries into this picture. In order to demonstrate the distributions of noise and data in the Fourier domain, the sampled Radon transform was applied to a digital phantom, as shown in Figure 2.4. e Figure 2.4: Digital Circle Phantom - Left: light background circle has an intensity of 1/4, light gray circles have an intensity of 1, dark gray are of intensity 2, and black circles have an intensity of 3. Right: corresponding noiseless sinogram. Chapter 2. The Fourier Slice Theorem and Computed Tomography 21 A discrete Fourier transform was applied for each of the projection angles. The transformed projection angles were then added together. Only the positive frequencies are shown in Figure 2.5. The left hand side of Figure 2.5 shows a semi-logi0 plot of f(kp) when the Radon transform of the phantom is corrupted by a random noise vector scaled so as to produce 16% error, on the mean. The right hand side of Figure 2.5 shows a semi-logi0 plot of f{kp) when no noise is present, and the Fourier transformed random noise vector. Frequency in Nyquist units (cycles / pixel) Frequency in Nyquist units (cycles / pixel) Figure 2.5: Radial Fourier Transform of Sinogram and Noise - Left: Absolute value of the discrete Fourier transform of a noisy sinogram. Right: Absolute value of the discrete Fourier transforms of noiseless data (solid) and random noise (dashed). The important features of the data and noise in Figure 2.5 are (i) the noise is uni-formly distributed over all frequencies (sometimes called white noise), and (ii) there is an obvious break-point where the data containing useful information rises above the constant-amplitude experimental noise. At low frequencies, the noise is orders of magni-tude less intense than the useful information. The data is more densely sampled at low frequencies, resulting in higher signal strength. Any data acquired by a system which can be approximated by the sampled Chapter 2. The Fourier Slice Theorem and Computed Tomography 22 Radon transform, and whose measurements are corrupted by uncorrelated Poisson noise will have these properties. 2.4 T H E MODULATION TRANSFER FUNCTION The most complete specification of an imaging device's spatial resolution is expressed by the modulation transfer function (MTF) [3]. It describes the fidelity with which a pattern will be imaged, as a function of the spatial frequencies in the pattern. Expressed another way, it describes the spatial resolution of the imaging device in the Fourier domain. Sorenson ([3]) compares the MTF to a frequency response curve, used for evaluating audio equipment. In the case of the frequency response curve, pure tones are sent into the audio component, and the relative amplitude of the output signal is measured. The frequency response curve is the ratio of output and input signal amplitudes over a range of frequencies. A system with a constant or "flat" curve over the range of frequencies provides the most faithful sound reproduction. In the case of an imaging system, the contrast of the imaged test pattern is measured: C = TJ—^7~T' ( 5 ) where Imax and Imin are the maximum and minimum radiation intensities emitted by the test pattern, as shown in Figure 2.6. The Modulation Transfer Function of an imaging system is defined as: MTF(k) = Cout(k)/Cin(k), (2.6) where the test pattern shown in Figure 2.6 is used to determine the MTF at a specific value of k. Each component of the imaging system has it's own MTF, and the system's total MTF is given by their product. [3] Chapter 2. The Fourier Slice Theorem and Computed Tomography 23 Figure 2.6: Modulation Transfer Function Test Spatial Frequency (cycles / mm) Figure 2.7: High Resolution Low Energy Modulation Transfer Function - line source on detector face (A), 20 cm from detector (B) and 55 cm from detector (C) Chapter 2. The Fourier Slice Theorem and Computed Tomography 24 Another method for measuring a system's MTF is to acquire it's Line Spread function (LSF) and take the Fourier transform of the LSF to obtain the MTF. This method was used to produce the MTF of an Anger camera equipped with a Low Energy High Resolution collimator, shown in Figure 2.7. 2.5 T H E R A D O N T R A N S F O R M AND T H E S Y S T E M M A T R I X In practice, the acquired data and the image being reconstructed are represented by lists of numbers which can be arranged as vectors. The transformation from the image vector to the corresponding data vector can be represented by a matrix operator.. This matrix models the acquisition system, and is referred to as the system matrix. The pure Radon transform can be thought of as modeling an Anger camera with a perfect collimator, no attenuating medium, and no scatter. In this case, an individual element of the system matrix would correspond to the Radon transform of image pixel j at detector bin i. Each detector bin corresponds to a particular (9,p) co-ordinate. A reconstruction of a 128 x 128 image from a 128 bin x 32 projection sinogram might require a system matrix consisting of over sixty-seven million elements. However, because an individual detector bin only "sees" a small number of the image pixels most of the system matrix's elements are zero. The system matrix based on the pure Radon transform is sparse, allowing it to be stored in a form which requires much less computer memory than if all the matrix elements were non-zero. Figure 2.8 shows the sparsity structure of a Radon matrix. Chapter 2. The Fourier Slice Theorem and Computed Tomography 25 0 50 100 150 200 250 nz = 4904 Figure 2.8: Sparsity Structure of the Radon Matrix 2.6 C O L L I M A T O R BLURRING With the pure Radon transform, a point source will produce a point projection, while the actual Anger camera would image a point source as a Gaussian distribution on the camera face (Figure 2.9). This effect is due to the collimator's acceptance angle, and the resolution of the scintillation crystal and PMT array. The distribution that a point source produced on the camera face is called the Point Spread Function (PSF), and the LSF is distribution a line source produces on the camera face. The width of the LSF and PSF has two components, one due to the collimator's acceptance angle, which depends on the distance from the source to the collimator, and one due to the resolution of the scintillation crystal and supporting electronics, which depends on the nature of the crystal and the photon energy and is independent of the Chapter 2. The Fourier Slice Theorem and Computed Tomography 26 source distance. The PSF can be described by the equation PSF(x, z) = exp x2 + z2 2(v2 + *2Ay))\ ( 2 ' 7 ) where crt is the intrinsic response of the scintillation crystal, PMTs and electronics, and <Tc is the depth dependent collimator response. [16] detector source Figure 2.9: Collimator Blurring Geometry 2.7 A T T E N U A T I O N In order to model the attenuation effect described by Equation 1.2, the Radon transform can be replaced by the attenuated Radon transform A(z,P,e) = exp [ -J r ^^M0*(p-<* ,O)#] U(e,p) = [nMe,P) = fA(tP,o)f{£)6(jp-(o,z))dt. where 0L is a unit vector perpendicular to 0, pointing along L in Figure 2.10, and £' • 01 > £ • 6L denotes the segment of L in between the point £ and the physical camera face (Denoted L' is Figure 2.10). [14] Chapter 2. The Fourier Slice Theorem and Computed Tomography 27 The properties of this transformation have been extensively studied (see [17], for example), and it's representation as a matrix operator is straightforward and has the same sparsity structure as the matrix representation of the pure Radon transform. C H A P T E R 3 S P E C T R E C O N S T R U C T I O N M E T H O D S 3 .1 F I L T E R E D B A C K - P R O J E C T I O N ( F B P ) Among the essential features of the FBP method are its linearity and the fact that a certain level of resolution in the reconstructed image is specified by means of the filtering operation. One of the operations that makes up the FBP method is the back-projection, which can be represented by the adjoint of the Radon transform, TI*: [n*f](t) = / f(e,(o,s))M, (3.i) Jo [14] where f(9,p) is an object in the Radon domain. Because the Radon transform is related to line integral of transformed function with-out the incorporation of any new imaginary term, / will only be complex if / was complex to begin with. For this reason, the discrete matrix representation of TZ has no imaginary component. Because the adjoint of R corresponds to the Hermitian conjugate of R, and R has no imaginary component, R* -» R} [18]. The action of the Radon transform and it's adjoint, the back-projection operator are illustrated in Figure 3.1. 28 Chapter 3. SPECT Reconstruction Methods 29 (A) (B) fC> 20 40 60 80 100 120 20 40 60 80 100 120 0.02 0.015 0.01 0.005 0 Figure 3.1: Back-Projection Operator - (A) the original object /(£) (B) its sampled Radon transform [/](#i..4,p) (C) the back-projection [R*/](£) In a continuous setting, the FBP reconstruction is defined by i=K*(g*f), (3.2) where / represents the Radon transform of the original object, f is the reconstructed image, g{0,p) is a convolution kernel in the Radon domain (described below), and * is a convolution with respect to the radial variable, p. The FBP method therefore essentially consists in applying the adjoint of the Radon transform to convolved projection data. The convolution of two functions u(x) and v(x) is defined u(x)-kv(x) = J u(x)v(x — y)dy. (3.3) Therefore, the convolution of two functions in the Radon domain with respect to p is u{e,p)*v(6,P) = ju(9,p)v(d,p-q)dq. (3.4) It can be shown ([19]) that a convolution in the spatial domain corresponds to a multiplication in the Fourier domain: /inf U(k)V(k)ei2*{Kr)dk, (3.5) •inf Chapter 3. SPECT Reconstruction Methods 30 where $ _ 1 represents the inverse Fourier transform and U and V are the Fourier trans-forms of u and v. Using Equation 3.5, a discrete FBP reconstruction can be expressed x = ff^^GY]], (3.6) where x is the vector representing the reconstructed image, G and Y are the Fourier trans-forms of the convolution kernel and acquired data with respect to the radial variable p, and <3>~1 represents the inverse Fourier transform with respect to the radial variable p. As is apparent from Figure 2.3, the discrete version of / has gaps between the spokes whose size are proportional to the radial distance from the origin, ||k||. This effect is due to the polar sampling of the Radon transform. One of the important functions of the convolution kernel g is to weight high frequencies more heavily than low, the equivalent of replacing dxdy with rdrdd for integration in polar co-ordinates. The component of g which performs this function is referred to as the ramp filter. The spatial equivalent of the ramp filter is the sine function, shown in Figure 3.2. The convolution kernel g can be seen as being composed of two parts: g = rs, where r is the ramp filter and s is the noise dampening part of the filter, described below. Having understood the function of the ramp filter, reconstructing the image / can be thought of as equivalent to interpolating the Fourier transform of / , i.e. to synthesizing F in the angular gaps of the experimental frequency coverage (see Figure 2.3). The second function of the convolution kernel g is to dampen the effects of exper-imental noise on the reconstruction. As shown in Section 2.3, the experimental error dominates the data's high spatial frequencies. For this reason, the noise dampening part of the convolution kernel in the Fourier domain s is simply a low-pass filter. Because part of g is a smoothing kernel in the spatial domain, the expected result of the FBP is not to reconstruct / directly but a resolution-limited (smoothed) version, Chapter 3. SPECT Reconstruction Methods 31 0.25 0.2 0.15 0.1 0.05 -Q.05[ -0.1 50 100 150 200 Figure 3.2: Sine Function s*f. It is interesting to note that, as pointed out in [20] and [21], the FBP algorithm can be derived from the following fundamental relationship: where * is a normal 2-dimensional convolution (as opposed to *, the convolution with respect to p in the Radon domain). 3.2 RECONSTRUCTION AS AN INVERSE P R O B L E M Another approach to the image reconstruction process is to formulate it as an inverse problem. In this case, we define an objective function, o(x), whose minimum corresponds to the "best" possible trial image (x), given the data acquired (y) and a system matrix (R). The reconstruction problem then becomes a minimization problem where various trial images are tested until the one that minimizes o is found. •(K*g)*f = n*(g*Hf) (3.7) x = argmin{o(x) = e(i?x,y) x G R n , x G (3.8) Chapter 3. SPECT Reconstruction Methods 32 where argmin means the argument minimizing o(x) subject to the constraints x G R n and x G V (V is a sub-set of W1 which may contain positivity constraints), and e is a function which measures the difference between the data vector corresponding to x (Rx) and the data that was actually acquired (y). Stated in this way, the reconstructed image (x) is the x that minimizes o. The reconstruction process then becomes a search for the image x which minimizes the objective function. The field of mathematics which addresses this type of problem is called optimization theory. There exist a number of optimization algorithms which are designed to find the value of a variable which minimizes an objective function. The optimization algorithms which use the gradient of the function being minimized are said to be quasi-Newtonian, in reference to the Newton-Raphson method for finding the roots of an equation. The Expectation Maximization algorithm from mathematical statistics is another optimization method used to solve inverse problems. 3.2.1 L E A S T SQUARES A N D W E I G H T E D L E A S T SQUARES In order to fully define o, the fit function e must be selected. A common choice is the squared L2 norm of the difference between Rx and y: els(Rx, y ) = || Jfcc - y ||2 = ]T([ifcc]i - yi)2, (3.9) i also known as the squared Euclidean norm or the least squares function. The gradient (needed by quasi-Newtonian minimization algorithms) of eis with respect to x is: Vxels{Rx, y) = 2R* \\ [Rx] - y ||. (3.10) The least squares fit term is a good starting point, as it makes no assumptions about the quantities being compared. Chapter 3. SPECT Reconstruction Methods 33 A variation on the least squares fit term is the weighted least squares fit term proposed in [22]. This fit term weights the fit to each sinogram bin according to the inverse of the measurement's variance: ^ ( y - A O ' E - ^ y - i f c ) (3.11) where E is a diagonal matrix with the ith entry being of, an estimate of the variance of the ith measurement y{. This fit term is tailored to PET imaging, where the statistics of the measurements can be described by Gaussian statistics due to the nature of the radiation counting system used and the corrections applied to the data after acquisition [22]. The maximum likelihood fit term which is based on Poisson statistics will be used instead of weighted least squares, as the data reconstructed in this work obeys Poisson statistics. 3.2.2 M A X I M U M L I K E L I H O O D In the case of tomographic imaging, the quantities dealt with are most often observations of radioactive decay, a stochastic process which can be described by Poisson statistics. The number of photons emitted from the jth voxel in the source distribution (XJ) can be regarded as an independently distributed Poisson variable. The number of photons observed by a single detector element yi is a sum of XjS weighted by the system matrix: yt = [Rx\i = V \ Ri,jxj a n d will therefore also obey Poisson statistics [23]. In this case, the re-projection of the reconstructed image (i?x) can be seen as the set of Poisson distribution means that is most likely to generate the observed count distribution y. The probability of observing yi counts from a Poisson distribution with mean pi = [Rx]i, Chapter 3. SPECT Reconstruction Methods 34 for a collection of observations (y) is [24]: (3.12) where / / j is the sum of the Poisson means contributing to data element Instead of attempting to maximize this probability directly, it is more practicalto deal with the log of this function [24]: Also, the sign of the function is changed so that the maximization becomes a minimization and the ln(yi\) term is dropped because it is a constant that does not depend on / i ; : Taking the gradient of this function with respect to x (for use with quasi-Newtonian algorithms) gives us [23]: 3.3 ITERATIVE METHODS Once the objective function is defined, the next step is to select an algorithm with which to find the image x which minimizes it. This is an important feature of the reconstruction process which influences computational efficiency. 3.3.1 M L - E M A N D OS-EM The Maximum Likelihood Expectation Maximization (ML-EM) reconstruction method was first formulated by L. Shepp and Y. Vardi in 1982 [23]. It consists of using the lnP(y,p) = Y^iVilnPi] - ^ P i ~ J^ Ms/iO- (3.13) emi(Rx,y) = -lnP(y,Rx) - J^fn^!) = ]P[[-Rx]; - yiln[Rx\i . (3.14) (3.15) Chapter 3. SPECT Reconstruction Methods 35 Figure 3.3: ML-EM and OS-EM - (A) 1 iteration of the M L - E M algorithm (B) 32 iterations of ML-EM (C) 4 iterations of OS-EM with 8 sub-sets Expectation Maximization algorithm to optimize the maximum likelihood fit function presented in Section 3.2.2. Ordered Subset Expectation Maximization (OS-EM) is a variation of ML-EM first presented in [25]. It consists in the application of ML-EM to ordered sub-sets of the data being reconstructed. In one iteration of OS-EM, all the sub-sets of the data are used, but the algorithm's convergence rate is increased by an order of magnitude as shown in Figure 3.3. Because of noise in the reconstructed data (y), the optimization problem that these algorithms attempt to solve is ill-posed, meaning that there is no source distribution x that will exactly satisfy the equation y = Rx. As described in Section 2.3, noise domi-nates the high frequencies of the experimental data. The consequence of the reconstruction problem's ill-posedness is that as the ML-EM and OS-EM algorithms progress, the image is degraded by high frequency artifacts. This can be seen in Figure 3.4. The effects of noise in the reconstructed data are obvious Chapter 3. SPECT Reconstruction Methods 36 O S - E M 8 subsets, 1 iteration O S - E M 8 subsets, 4 iterations O S - E M 8 subsets, 8 iterations Figure 3.4: Increasing Iterations of OS-EM (with 8 sub-sets) when the middle image in Figure 3.4 is compared to the noiseless reconstruction shown in Figure 3.3 C. In practice, the ML-EM and OS-EM algorithms are stopped before the reconstructed image begins to show these effects. A post-reconstruction filter is sometimes applied in order to further reduce these effects. C H A P T E R 4 R E G U L A R I Z E D S P E C T R E C O N S T R U C T I O N Up to this point, we have described objective functions consisting only of a fit term. A slightly more complex formulation of the problem uses an objective function consisting of two terms: x € W1 >, (4.1) x = argmin |o(x) = ^e(Rx,y) + ^ap(x) where e is the fit function, p is the regularizer (whose negative is often called an entropy), and a is a scalar parameter which controls the balance between the two terms. The role of the regularizer is to compensate for the ill-posedness of the reconstruction problem when formulated using the fit term alone. In the case of SPECT reconstruction the two principal sources of instability are (i) statistical noise in the acquired data and (ii) the incompatible sampling patterns in the Fourier domain (polar data co-ordinates vs. Cartesian image co-ordinates). 37 Chapter 4. Regularized SPECT Reconstruction 4.1 CLASSICAL REGULARIZERS 38 When selecting or defining a regularizer, care must be exercised to create a regularized objective function which is concave. In practical terms this means that the function should be designed to have a unique, well defined,.minimum. Examples of classical regularizers are the Boltzmann-Shannon, the Kullback-Leibler, the Tikhonov and the generalized quadratic Tikhonov neg-entropies, respectively [26]: n pbs(x) = ^Xjlnxj, pkl(x) = y^Xjln^-, j=i j=i x°i • (4.2) Pii(x) = | |x|| 2 , pq(x) = (x,Qx). In the definition of p^i, the XQJ'S are the components of some reference model of the object. In the definition of pq, Q is a symmetric positive semi-definite matrix, which should be chosen so that the resulting objective function is concave, as stated above. The gradients of the regularizers listed above (for use with quasi-Newtonian mini-mization algorithms) are [26]: V xp b 5(x) = lnx, Vxp fci(x) = lnx- lnxo, (4.3) V x/3tt(x) = 2x, V xp,(x) = 2Qx, with the In function understood to be component-wise (lnx = [Inxi, lnx2, lnx lnx j]). 4.2 FOURIER REGULARIZATION One of the principal features of Fourier regularization is the partitioning of the Fourier domain into two areas: one where the fit term extracts meaningful information, and another where the regularizer dampens frequencies that were under-sampled or not sam-pled at all. The Fourier domain is partitioned (see Figure 4.1) in order to ensure that Chapter 4. Regularized SPECT Reconstruction -iq •r^— 1 fii K Figure 4.1: Partitioning the Fourier Domain - The white area shows the part of the Fourier domain constrained by fit to the experimental measurements, and the shaded area covers the frequencies constrained by the regularization term. the fit and regularization terms are not in direct competition with each other during the minimization process. Like FBP, the objective of the reconstruction process is to find a resolution-limited version of the source distribution s * f (where s is a convolution kernel) instead of / itself. The corresponding Fourier domain filter will be denoted S, and the corresponding convolution kernel in the Radon domain will be denoted s. The Fourier domain is partitioned by applying a low-pass filter to the acquired data before the minimization (s * <f>), and defining a regularizer that returns a measure of the high frequencies in / , ignoring the resolution range where the experimental information is adequate: (4.4) Chapter 4. Regularized SPECT Reconstruction 40 The reconstructed image is therefore defined as the solution to the following opti-mization problem: / ss argmin j o(f) = l-e{Tlf, s * /) + 11| (1 - S)F \ f = f(v)eV,TERs }. (4.5) In a discrete setting, this becomes x = argmin |o(x) = ^e(Rx, Csy) + | | |B s x| | 2 X € K n , x e 7> j , (4.6) where Cs is the matrix representation of the s * • operation and Bs is the matrix rep-resentation of the Fourier transform and the application of the 1 — 5 filter, and V is a sub-set of E" which may contain positivity constraints. The boundary dividing the areas contributing to the regularization and fit terms is radially symmetric in the Fourier domain. This feature is important with regard to the visual interpretation of the reconstructed image. A boundary that is non-circular in the Fourier domain will result in an image with varying spatial resolution at different angles. II l|2 The Fourier regularizer p/ = Bsx || is a special case of the quadratic regularizer pq with Q = B*BS. The gradient of pj is therefore: V x P / = 2B*Bsx. (4.7) C H A P T E R 5 S E N S I T I V I T Y T O N O I S E In this chapter the effect, on the reconstructed image, of experimental noise in the ac-quired data will be analyzed for a selection of reconstruction methods. For some methods the sensitivity to noise can be studied through mathematical analy-sis, but other reconstruction methods can not be described rigorously enough for analytic techniques to be applied. Where applicable, the analytic approach offers greater insight into a method's sensitivity to the reconstructed data. ML-EM and OS-EM are not suitable for the sensitivity analysis described here because they are stopped before a solution to the optimization problems described in Equation 3.8 is found. As will be seen, the sensitivity analysis requires that iterative methods converge (V x/(x) = 0). Using sensitivity analysis, a quantitative measure of a reconstruction method's sensi-tivity to noise is only available for small scale two dimensional reconstructions. Empirical methods are more suitable when a quantitative measure of a reconstruction method's re-sponse to noise (such as the covariance matrix) is desired. All methods can have their response to noise studied empirically by repeatedly recon-structing data that are corrupted by artificially generated random noise vectors. This 41 Chapter 5. Sensitivity to Noise 42 type of study will described in Chapter 6 and the corresponding results will be presented in Chapter 7. The image co-variance matrix Cx calculated in Chapters 6 and 7 is related to the sensitivity matricies described below by: CX = SCyS* , (5-1) where Cy is the sinogram co-variance matrix. [26] 5.1 F B P The discrete FBP reconstruction can be described as x = R*&-1g*py (5.2) where R*, ^ p 1 , g, and <frp are matrix representations of the elements making up the FBP algorithm as described in Equation 3.6. If y is replaced by y + 6y where 6y represents experimental noise, and x replaced by the perturbed solution x + Sx, Equation 5.2 becomes x + dx = R+Q-igGpty + Sy). (5.3) Subtracting Equation 5.2 from 5.3, it becomes apparent that the sensitivity matrix, Sfbp, relating experimental noise to perturbations in the reconstructed image is the same as the matrix representation of the FBP algorithm: Sx = Sfbp6y, with Sfbp = i i i * * - 1 ^ . (5.4) The sensitivity analysis of FBP reconstruction is extremely simple because FBP is a linear operation and can be represented by a matrix operator. Chapter 5. Sensitivity to Noise 43 5.2 ITERATIVE M E T H O D S The starting point for the sensitivity analysis of iterative methods is the condition that the problem expressed in Equation 4.1 be solved. In this case, V x / (x) = ^-Vxe(R5t, y) + | Vxp(x) = 0. (5.5) Replacing y with y + Sy, the same condition must hold with x replaced by the perturbed solution x + Sx: Y Vxe(i?[x + Sx],y + Sy) + | V xp(x + 8x) = 0. (5.6) Keeping only the first order terms in the development of Equation (5.6) yields (using Equation (5.5)) [R*H£R + aHe] Sx = -R*H£ Sy, (5.7) i in which we have denoted H£ and He the Hessian (matrix of second partial derivatives) of £ at (Rx,y) and the Hessian of g at x, respectively. We therefore obtain Sx = SSy with S = -[R+HcR + aHg]'1 R*H£. (5.8) The sensitivity matrix S depends on the Hessians of e and p, listed in the table below: Chapter 5. Sensitivity to Noise 44 Vxels(Rx,y) 2 i ? * | i ? x - y | | = 2R*R Vx;fm/(Z?x,y) = £ ^ [ 1 - ( * / [ « * ] , ) ] Hem, = Vi/x] V xp b s(x) = lnx H P b s = 1/x Vxpfci(x) = ln(x/x0) H p k l = x 0 /x VxPti(x) = ' 2x H P U 2 VxPg(x) 2Qx 2Q V X P S ( X ) 2B*B 8x = 2/3;^ Table 5.1: Gradients and Hessians of e and p C H A P T E R 6 R E C O N S T R U C T I O N S 6.1 SIMULATION OF S P E C T DATA During the course of research activities related to image reconstruction, data is often generated by using a computer to mathematically simulate the source distribution and detection system. This approach allows the physical situation to be modeled with varying degrees of detail. Physical effects such as scatter, attenuation, collimator blurring and statistical noise can be modeled in different ways, or not modeled at all. In this work, data is simulated in two ways, described below. 6.1.1 A N A L Y T I C SIMULATIONS In this case the camera's action is approximated by the system matrix which may be a matrix representation of the Radon transform or one of its variations described in Sections 2.5, 2.6, and 2.7. The source distribution is represented by a two or three dimensional matrix, as is the acquired data. The data is produced by carrying out the matrix multiplication y = Rx, where y is the acquired data, R is the system matrix and x is the matrix representing the source distribution. 45 Chapter 6. Reconstructions 46 This approach produces data without any statistical noise, so random noise vectors may be added to the data to simulate this effect. Attenuation and collimator blurring may be modeled in the system matrix, or the system matrix may be based on the pure Radon transform. Because of practical limitations, scattering effects are generally not modeled in the system matrix. 6.1.2 M O N T E - C A R L O M E T H O D S Monte-Carlo methods use the basic physics of radiation transport to calculate the paths of individual photons through a medium. Random numbers are used to determine the outcome of all probabilistic effects. A large number of photons are tracked in this way, until the effects of statistical noise have been reduced to the desired level. Monte-Carlo simulations require more time to perform than the analytic simulations described above, but they can model the physics of the problem more accurately than other techniques. The Monte-Carlo simulation presented in this was generated using the SIMSET Monte-Carlo program [27]. 6.2 DIGITAL P H A N T O M S 6.2.1 A N A L Y T I C C I R C L E S The phantom that was used in the analytic simulations is shown in Figure 6.1. The light background circle has an intensity of 1/4, the light gray circles have an intensity of 1, the dark gray are of intensity 2, and the black circles have an intensity of 3. The sinograms generated from the analytic projections of this phantom are scaled so as to have an average of 5000 counts per projection. Gaussian distributed noise is added Chapter 6. Reconstructions 47 to the sinogram, with mean 0 and standard deviation corresponding to the square root of the number of counts in each bin. Figure 6.1: Digital Circle Phantom - Light background circle has an intensity of 1/4, light gray circles have an intensity of 1, dark gray are of intensity 2, and black circles have an intensity of 3 6.2.2 M O N T E - C A R L O M C A T P H A N T O M One slice of the activity distribution used for the Monte-Carlo simulations is shown in Figure 6.2. This activity distribution is an approximation of the human heart, generated by the MCAT digital phantom program [28]. The projection data produced by the Monte-Carlo simulation of this phantom (with-out modeling an attenuating medium) has an average of approximately 6100 counts per projection. Because of the nature of the Monte-Carlo simulation each bin of the projec-tion data is corrupted by random Poisson noise with a standard deviation proportional to the square root of the number of counts in the bin, as would be the case with actual radiation. Chapter 6. Reconstructions 48 Figure 6.2: One Slice of the MCAT Phantom's Activity Map 6.3 R E C O N S T R U C T I O N E V A L U A T I O N C R I T E R I A In order to quantitatively compare reconstructions, the two criteria described below were used. However, the final judgment is reached by qualitative comparison. Certain char-acteristics of each reconstruction may be expressed by the numerical criteria described here, but these results must be interpreted by a human. The criteria described below are used to select candidate images among the set of images produced by varying cutoff frequency and regularization parameters. These can-didate images are then compared and discussed. 6.3.1 D I F F E R E N C E F R O M ORIGINAL O B J E C T A N D R E - P R O J E C T I O N M I S F I T TO D A T A The room mean square (RMS) difference between the original object and the recon-structed image is one of the most direct and intuitive ways to evaluate the quality of a reconstruction, but can only be used for reconstructions of simulated data, where the original source distribution is known. The recipe used to calculate the RMS difference between the original object and the Chapter 6. Reconstructions 49 reconstruction is as follows. First, both the reconstructed image and the original source distribution are normalized so that their L 2 norms are one: v = y/VEiVl where the index i runs over all elements in the original object and the reconstruction, x and y. The norm of the difference between u and v is referred to as the RMS difference from the original object. One might expect the optimal reconstruction to be the one where this characteristic value is as low as possible, but this criterium cannot be relied on for experimental studies because the original object is inaccessible in practice. There will be further discussion of this topic in Chapter 7. The second characteristic that is calculated for each reconstruction is the re-projection misfit to the data. The re-projection of the reconstructed image is generated by applying the system matrix to the reconstructed image. The re-projection is then normalized to have the same norm as the sinogram: r = Rx x \\y\\/ \\Rx\\ (6.2.) The re-projection misfit to the data is then: misfit =\\y-r\\/\\y\\. (6.3) This characteristic is more generally applicable because it compares the reconstructed image to the acquired data, which is always available. The expected value of this char-acteristic depends on the noise present in the data, which is known due to the Poisson nature of radiation statistics [23]. What is meant by "expected value" is the misfit that would arise if the reconstructed image were an exact representation of the original object, and the re-projection were Chapter 6. Reconstructions 50 equivalent to noiseless data. In this case, the only factor contributing to the misfit would be the statistical noise present in the acquired data. In the results presented in Chapter 7, the expected value of the misfit is calculated by generating one hundred random noise vectors where each element's standard deviation is the square root of the number of counts in the corresponding sinogram bin. The expected misfit is then calculated: 100 <mts/tt>=^| |ni | | / | |y | | , (6.4) i=i where is the ith noise vector. 6.3.2 NOISE R E S P O N S E The response of each of the reconstruction methods to noise is calculated using the technique presented in [29], where an ensemble of noisy sinograms are reconstructed and the variance and co-variances of each pixel in the reconstructed images can be calculated. Each element of the co-variance matrix C is: [29] C « J = jv3rD/r- /ow-/ i ) (6.5) n=l where / " is the ith pixel value in the nth image from a sample of N noisy images, and n=l Cjj then refers to the co-variance of the ith and jth pixels in the reconstructed image, and Cij would be the variance of the ith pixel. It can be shown that for Gaussian distributed noise the relative error of the variance estimate is y/2/(N — 1) [29]. For this study, the sample size was 200 images, so we would expect a relative error in variance of about 10%. Chapter 6. Reconstructions 51 6.4 RECONSTRUCTION P A R A M E T E R S 6.4.1 S E L E C T I N G C U T O F F F R E Q U E N C Y A N D F I L T E R For all of the reconstruction methods tested here (FRCT, O S - E M and FBP), a Hann filter was used. This selection was made arbitrarily (any filter shape could have been used), but was kept constant for all reconstructions to facilitate comparison. The Hann filter is defined as: in the data should guide the selection of filter shape and cut-off frequency. As can be seen in Figure 2.5, there is a discernible break-point between the constant-amplitude experimental noise and the acquired data. Methods of automatically selecting the cut-off frequency have been investigated (see [30]). In this work, the cut-off frequency is considered to be a free parameter. 6.4.2 R E G U L A R I Z A T I O N P A R A M E T E R S A N D STOPPING C R I T E R I A For the O S - E M reconstructions, the number of iterations is treated as a free parameter. In OS -EM reconstructions, the number of iterations is used (i) as the stopping criteriurri for the iterative algorithm, and (ii) as a regularization technique, mitigating the high-frequency artifacts that appear at higher numbers of iterations. A Hann filter is applied to the reconstructed image after reconstruction, and the filter's cut-off frequency is treated as a second free parameter. For the FRCT reconstructions, the regularization parameter and the minimization algorithm's stopping criteria are separate parameters. The regularization parameter is otherwise (6.7) where / is the input frequency and fc is the cut-off frequency. In practice, knowledge of the imaging system's MTF and the expected nose level Chapter 6. Reconstructions 52 treated as a free parameter, while the stopping criterium is kept constant for a given data-set and reconstruction method. The FR-LS and FR-ML reconstructions of the analytic circle data are considered to have converged when the norm of the objective function's gradient has been reduced by eight orders of magnitude from it's initial value. The FR-LS reconstructions of the Monte-Carlo MCAT data are stopped when the norm of the objective function's gradient has been reduced by six orders of magnitude, and the FR-ML reconstructions of the same data are stopped when the objective functions' gradient has been reduced by seven orders of magnitude. The stopping criteria for the FRCT minimizations were selected by trial and error. The stopping criterium was made more stringent until the gradient could not be reduced any further, and then slightly relaxed. 6.4.3 P A R A M E T E R R A N G E S The analysis described in Section 6.3 was performed for reconstructions of each of the two phantoms with cut-off frequencies, iteration numbers and regularization parameters varying over a range of values, listed in Table 6.1. The values of the RMS difference between reconstruction and original object and the misfit between re-projection and sinogram are plotted versus the reconstruction parame-ters. For the FBP reconstructions this results in two curves, both plotted versus fc, for each of the two phantoms (four curves in all). For of OSEM, FR-LS and FR-ML, where there are two free parameters, surface plots are created when the characteristics discussed in Section 6.3 are plotted versus the reconstruction parameters. In these cases, one curve is extracted from each of the surfaces. For the RMS difference between reconstruction and original object, the number Chapter 6. Reconstructions 53 Analytic Circle Phantom M C A T Phantom FBP fc = 0.1 ... 1 cycles / pixel in steps of 0.1 fc = 0.1 ... 1 cycles / pixel in steps of 0.1 OSEM (four subsets) 1 ... 15 iterations in steps of 1 1 ... 11 iterations in steps of 1 fc = 0.1 ... 1 cycles / pixel in steps of 0.1 fc = 0.1 ... 1 cycles / pixel in steps of 0.1 RC-LS a = lO^.-lO 5 5 on a log scale in steps of IO 0 5 a = 1...107 on a log scale in steps of IO 0 5 fc = 0.2 ... 1.8 cycles / pixel in steps of 0.2 fc = 0.2 ... 1.8 cycles / pixel in steps of 0.2 RC-ML a = 10_1...104 on a log scale in steps of IO 0 5 a — 10~15...105 on a log scale in steps of IO 0 5 fc = 0.2 ... 1.8 cycles / pixel in steps of 0.2 fc — 0.2 ... 1.8 cycles / pixel in steps of 0.2 Table 6:1: Reconstruction Parameter Settings Used in Tests of Reconstruction Methods of iterations (for OSEM) or the regularization parameter (for FRCT) which results in the lowest value at a given cutoff frequency is selected. For the misfit between re-projection and sinogram, the number of iterations or the regularization parameter which results in the value closest to the expected value at a given cutoff frequency is selected. The noise responses of the reconstruction methods are calculated using the sets of parameters which were selected to reconstruct the two phantoms. The norm of the variance image is plotted versus fc. Due to time constraints, the entire co-variance matrix was not calculated - instead, the variance of each pixel is calculated and plotted versus the pixel's mean value for different values of fc. The spatial distribution of the variances is also shown. Finally, the co-variance image of the pixel located at the center of the image is displayed. C H A P T E R 7 R E S U L T S A N D D I S C U S S I O N Tables 7.1 and 7.2 show the figure numbers for each combination of phantom, reconstruc-tion method, and set of results. For the reconstruction methods with two free parameters, the noise response was calculated for the parameter settings that resulted in the least difference from the original object, and for the settings that resulted in the data misfit being as close to the expected value as possible. The surfaces that were generated by plotting the characteristics discussed in Sec-tion 6.3 versus two reconstruction parameters are presented without further discussion in Chapter 9. The features discussed in Sections 7.1 and 7.2 are summarized in Table 7.4 in Sec-tion 7.3. Difference from Original Object Re-projection Misfit from Data Noise Response FBP Circles Fig. 7.1 Fig. 7.1 Fig. 7.15 FBP MCAT Fig. 7.2 Fig. 7.2 Fig. 7.16 Table 7.1: FBP Results 54 Chapter 7. Results and Discussion 55 Least Difference from Original Object Noise Response Expected Re-projection Misfit from Data Noise Response OS-EM Circles Fig. 7.3 Fig. 7.17 Fig. 7.9 Fig. 7.18 FR-LS Circles Fig. 7.4 Fig. 7.19 Fig. 7.10 Fig. 7.20 FR-ML Circles Fig. 7.5 Fig. 7.21 Fig. 7.11 Fig. 7.22 OS-EM MCAT Fig. 7.6 Fig. 7.23 Fig. 7.12 Fig. 7.24 FR-LS MCAT Fig. 7.7 Fig. 7.25 Fig. 7.13 Fig. 7.26 FR-ML MCAT Fig. 7.8 Fig. 7.27 Fig. 7.14 Fig. 7.28 Table 7.2: Two-Parameter Results Chapter 7, Results and Discussion 56 7.1 F I T T O O R I G I N A L O B J E C T A N D F I T T O D A T A Comparing Figures 7.1 and 7.3 makes apparent the streaking artifacts that are typical of FBP reconstructions, especially in the low count areas of the image. In the reconstructions shown in Figure 7.4, it is possible to see similar streaking around the edges of the low intensity background circle. In this case, the streaking is due to the nature of the least-squares fit term, which weights the fit to low-count data bins differently than the maximum likelihood fit term. The streaking is not present in the FR-ML reconstructions shown in Figure 7.5. The streaking in the background of the FR-LS reconstructions is present for both phantoms, with both sets of parameter settings, supporting the hypothesis that it is a characteristic of this reconstruction method. The uniform areas in the reconstructions presented in Figures 7.3 and 7.5 show the difference in noise textures between OS-EM and FRCT reconstructions, especially in the low intensity background circle. The more splotchy low intensity background circle in the FRCT reconstructions may be the main reason that the OS-EM reconstructions achieved a lower difference between the reconstruction and the analytic circle phantom. However, when the MCAT reconstructions of Figures 7.6 and 7.8 are compared it is the OS-EM reconstructions that appear more splotchy. This leads one to conclude that the effect being observed depends on the object being reconstructed. The resolution of the cold spot within the ventricle can be compared in Figures 7.12 and 7.14. It is apparent that OS-EM is not achieving the same degree of contrast that the FR-ML reconstruction is. This observation is supported by a comparison of the small cold spot in Figures 7.3 and 7.5, for example. Over-all these results also show that selecting the reconstruction parameters that correspond to the least difference from the original object result in noisier and higher resolution reconstructions than parameter selections aimed at producing a data misfit as Chapter 7. Results and Discussion 57 close to the expected value as possible. See Figure 7.2, for example, and note that the minimum of the plot on left occurs at a higher fc than the intersection of the two curves in the right hand plot. Comparing the reconstructions of the analytic circle phantom, the author would rank OS-EM first, followed by FR-ML, then FR-LS, and FBP last. Based on the MCAT reconstructions, FR-ML would be ranked first, then OS-EM, followed by FR-LS and FBP. There will be further discussion in Section 7.2. (A) 0.2 cycles / pixel (B) 0.4 cycles / pixel (C) 0.8 cycles / pixel ! |m 0.5 RMS difference from original object vs. f Reprojection misfit to data vs. f t5 0> CO c I 0.4 g 80.35 3= S 0.3 CO 0.25 \ \* \ , x • \ B / 0.2 0.4 , 0.6 . , 0.8 1 in cycles / pixel r in cycles f pixel Figure 7.1: Selected FBP Reconstructions (top), RMS Difference Between FBP Recon-structions and Analytic Circle Phantom vs fc (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and Analytic Circle Sinogram vs fc (bottom, right) Chapter 7. Results and Discussion 58 (A) 0.2 cycles / pixel (B) 0.4 cycles / pixel (C) 0.8 cycles / pixel R M S d i f f e r e n c e f r o m o r i g i n a l o b j e c t v s . f R e p r o j e c t i o n mis f i t t o d a t a v s . f cr- 0.7 0 . 4 0 . 6 0 . 8 f in c y c l e s / p i x e l Figure 7.2: Selected FBP Reconstructions (top), RMS Difference Between FBP Recon-structions and MCAT Digital Phantom vs fc (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and MCAT Monte-Carlo Sinogram vs fc (bottom, right) Chapter 7. Results and Discussion 59 (A) 0.5 cycles / pixel (B) 0.6 cycles / pixel (C) 0.9 cycles / pixel ^ Least R M S Difference From Original Object Corresponding Iterations vs. f. o 0 o £ 0.23 ° 0.225 o © 0.22 o c 2> 0.215 0.4 0.6 0.8 f in Nyquist units (cycles / pixel) 12 CO ^ 8 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) Figure 7.3: Least RMS Difference Between 4 Subset OSEM Reconstructions and Digital Circle Phantom vs fc (left) and Corresponding Iteration No. vs fc (right) - Sample reconstructions use iteration No. corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion GO (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel • Least R M S Difference from original object 0.35 \ Corresponding a vs. f a) o c OJ i _ OJ a co 0.3 rxO.25 0.5 1 1.5 f in Nyquist units (cycles / pixel) s O ) o 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.4: Least RMS Difference Between FR-LS Reconstructions and Digital Phan-tom vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corre-sponding to least difference from original object at selected fc Chapter 7. Results and Discussion 61 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel mL • — n Least R M S Difference from original object 0.351 Corresponding a vs. f o I 0.3 CD io.25 CC 8 ~ ° 1 o 0.5 1 1.5 f in Nyquist units (cycles / pixel) 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.5: Least RMS Difference Between FR-ML Reconstructions and Digital Phan-tom vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corre-sponding to least difference from original object at selected fc Chapter 7. Results and Discussion 62 (A) 0.5 cycles / pixel (B) 0.6 cycles / pixel (C) 0.9 cycles / pixel Least RMS Difference From Original Object Corresponding Iterations vs. f. 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) f in Nyquist units (cycles / pixel) C C Figure 7.6: Least RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs fc (left) and Corresponding Iteration No. vs fc (right) - Sample reconstructions use iteration No. corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 63 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel 0.5 1 1.5 0.5 1 1.5 f, in Nyquist units (cycles / pixel) f in Nyquist units (cycles / pixel) Figure 7.7: Least RMS Difference Between FR-LS Reconstructions and Digital Phan-tom vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corre-sponding to least difference from original object at selected fc Chapter 7. Results and Discussion RA (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Least RMS Difference from original object Corresponding a vs. f. f, in Nyquist units (cycles / pixel) fc in Nyquist units (cycles / pixel) Figure 7.8: Least RMS Difference Between FR-ML Reconstructions and Digital Phan-tom vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corre-sponding to least difference from original object at selected fc Chapter 7. Results and Discussion 65 (A) 0.5 cycles / pixel (B) 0.6 cycles / pixel (C) 0.9 cycles / pixel h ™ -Reprojection misfit to data closest to expected value (0.129) Corresponding Iterations vs. f ca "O a o.i3 I c o '•4—I o CD 0.128 0-0.126 O-CD CC 15 i 1 0 •i—! to CD 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) Figure 7.9: Least RMS Difference Between Re-projection of 4 Subset OSEM Reconstruc-tions and Sinogram Data vs fc (left) and Corresponding Iteration No. vs fc (right) -Sample reconstructions use Iteration No. corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 66 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel 9 tip ^tfftk j^tfhk TO " O Reprojection misfit to data closest to expected value (0.128) £ 0 . w E c 2 0 o cu 15 .14 CD ' CC CO .13 cc 0.5 1 1.5 f in Nyquist units (cycles / pixel) Corresponding a vs. f o cn3 o 0.5 1.5 f in Nyquist units (cycles / pixel) Figure 7.10: Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 67 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Reprojection misfit to data closest to expected value (0.128) Z 0.16 M— in I 0.155 c o "I 0.15 1.0.145 CD CC co 0.14 0.5 1 1.5 f in Nyquist units (cycles / pixel) O 1 0 -1 Corresponding a vs. f 0.5 1.5 f in Nyquist units (cycles / pixel) Figure 7.11: Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 68 (A) 0.5 cycles / pixel (B) 0.6 cycles / pixel (C) 0.9 cycles / pixel 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 fe in Nyquist units (cycles / pixel) f, in Nyquist units (cycles / pixel) Figure 7.12: Least RMS Difference Between Re-projection of 4 Subset OSEM Recon-structions and Sinogram Data vs fc (left) and Corresponding Iteration No. vs fc (right) -Sample reconstructions use Iteration No. corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 69 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel « n o 3=0.12 10.11 o a> 'o Reprojection misfit to data closest to expected value (0.107) Q. CU rr co rr 0.1 0.09 0.5 1 1.5 f in Nyquist units (cycles / pixel) Corresponding a vs. f s 5 dT ° 4 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.13: Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs / c (right) - Sample reconstructions use a corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 70 (A) 0.4 cycles / pixel (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Reprojection misfit to data closest to expected value (0.107) 50.115 | 0.11 c g 1 0.105 Q. <D rr co rr 0.1 0.5 1 1.5 f in Nyquist units (cycles / pixel) 5 4 O ) o ~ 2 1 Corresponding a vs. f 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.14: Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs fc (left) and Corresponding a vs fc (right) - Sample reconstructions use a corresponding to least difference from original object at selected fc Chapter 7. Results and Discussion 71 7.2 RESPONSE T O NOISE The plots of the norm of the pixel variances versus fc for FBP at the top of Figures 7.15 and 7.16 show a continuous relationship between the norm of pixel variances and the cut-off frequency. The same plots for OSEM at the top right of Figures 7.17, 7.18, 7.23 and 7.24 show variations in the norm of the pixel variances that appear to be discontinuous. Figure 7.23 shows more evidenco of the highly unstable nature of OSEM reconstruc-tions. In this case, the parameters were selected to result in the least difference between the reconstructed image and the original object, resulting in higher numbers of iterations than those used in Figure 7.24. Variance images A and B in Figure 7.23 show individual pixels with variances so high that they overshadow the variances of the other pixels. In practice, these artifacts would be avoided by limiting the number of iterations, and would be easy to distinguish from true features if they were present. Similar discontinuities are present to a lesser degree in the "Norm(V) vs. / c " plots of Figures 7.20, 7.21, 7.22, 7.25, 7.26, and 7.27. The degree by which the outlying points are removed from the mean are summarized in Table 7.3. It is appearent from a visual inspection of the OS-EM "Norm(V) vs. f" plots and from Table 7.3 that OS-EM responds in an unpredictable way to the presence of noise in the reconstructed data. The noise response of the other methods is not as varied, making the mean of the norms of the pixel variances at different cut-off frequencies a more meaningful characterization of the reconstruction method's response to noise. In order to better quantify the response of the OSEM algorithm to noise in the acquired data, a larger sample of reconstructions may be required. Another factor that could make the "Norm(V) vs. / c " plots more interpretable would be a finer sampling of the reconstruction parameters. This would enable us to select parameter settings that came closer to achieving the goals of matching the reconstruction Chapter 7. Results and Discussion 72 Mean Max - Mean (Max - Mean) / Mean Fig. 7.15 - FBP, Circles 1.3 x 104 2.8 x 104 2.2 Fig. 7.16 - FBP, MCAT 1.7 x 104 3.8 x 104 2.2 Fig. 7.17 - OSEM, Circles, Least Difference from Original Object 1 x 105 8.4 x 105 8.1 Fig. 7.18 - OSEM, Circles, Expected Misfit to Data 6.1 x 104 5.2 x 105 8.5 Fig. 7.19 - FR-LS, Circles, Least Difference from Original Object 3.4 x 103 1.3 x 103 0.39 Fig. 7.20 - FR-LS, Circles, Expected Misfit to Data 1.9 x 103 6.9 x 102 0.35 Fig. 7.21 - FR-ML, Circles, Least Difference from Original Object 4 x 103 1.8 x 103 0.45 Fig. 7.22 - FR-ML, Circles, Expected Misfit to Data 3.7 x 103 1.1 x 103 0.30 Fig. 7.23 - OSEM, MCAT, Least Difference from Original Object 3.1 x 104 8.2 x 104 2.7 Fig. 7.24 - OSEM, MCAT, Expected Misfit to Data 1.3 x 104 1.1 x 105 8.1 Fig. 7.25 - FR-LS, MCAT, Least Difference from Original Object 2.9 x 103 8.8 x 102 0.30 Fig. 7.26 - FR-LS, MCAT, Expected Misfit to Data 3.3 x 102 7.8 x 102 2.4 Fig. 7.27 - FR-ML, MCAT, Least Difference from Original Object 3.2 x 103 1.9 x 103 0.59 Fig. 7.28 - FR-ML, MCAT, Expected Misfit to Data 1.9 x 102 1.8 x 101 0.45 Table 7.3: Characteristics of "Norm(V) vs. / c " Plots Chapter 7. Results and Discussion 73 to the object or the re-projection to the data, which could result in more consistent pixel variances. Ideally in this type of study, all of the free reconstruction parameters should be allowed to vary over a range of values, but this was impossible due to time constraints. The results presented in Table 7.3 show that the parameters which attempt to generate the expected misfit to the data result in reconstructions which are more stable with respect to noise. This is consistent with the observation that these reconstructions tend to be smoother and have lower resolution, as discussed in Section 7.1. Also, the slightly higher number of counts (and correspondingly lower noise level) in the MCAT data results in slightly lower pixel variances for the iterative reconstruction techniques. The FRCT reconstructions have consistently lower pixel variances than the cor-responding OSEM and FBP reconstructions. The OSEM reconstructions had post-reconstruction filters with lower cut-off frequencies than the cut-off frequencies in the FRCT reconstructions which should give OSEM an advantage in this analysis, but this effect is overwhelmed by the underlying instability of the method. Comparing the shapes of the scatter plots in the second row of Figures 7.17 and 7.19, or any other OSEM and FRCT reconstructions is interesting. In the case of OSEM, there appears to be an linear relationship between the mean pixel value and corresponding variance, while for the FRCT reconstructions the stabilizing effect of the regularizer can be seen in the limited pixel variances. The co-variance images of Figures 7.1 and 7.2 show the relationship between correla-tion lengths in the reconstructed image and cut-off frequency. The co-variance images of Figures 7.26 and 7.28 show larger correlation lengths than any other combination of reconstruction method, phantom and method of parameter Chapter 7. Results and Discussion 74 selection. It has been previously noted that parameters selected to fit the reprojections of the reconstruction to the data result in smoother images which are more stable with respect to noise, and it is reasonable to attribute the longer correlation lengths to the same cause. Taking into account the noise responses of the four reconstruction methods presented here, the author would rank FR-ML first, followed by OSEM; FR-LS and FBP. Chapter 7. Results and Discussion 75 x 10 NormfV) vs. f > § 2 o C 1.5 1 0.5 C / / / / B / t 0.2 0.4. 0.6 I in cycles / pixel 0.8 Variance plot A (0.2 cycles / pixel) 20 40 60 80 Mean Pixel Value Semi-Log NormfV) vs. f 0.2 0.4. 0.6 0.8 f in cycles I f Variance plot B (0.4 cycles / pixel) 40 60 Mean Pixel Value / pixel Variance plot C (0.8 cycles / pixel) 40 60 80 100 Mean Pixel Value 120 140 Variance image A (0.2 cycles / pixel) Variance image B (0.4 cycles / pixel) 50 45 40 35 30 25 20 15 10 Variance image C (0.8 cycles / pixel) Co-Variance image A (0.2 cycles / pixel) Co-Variance image B (0.4 cycles / pixel) 20 40 60 80 100 120 Co-Variance image C (0.8 cycles / pixel) 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 200 150 100 50 |o -50 Figure 7.15: Noise Response of FBP Reconstructions vs. fc for Analytic Circle Data Chapter 7. Results and Discussion 7G Variance image A Variance image B Variance image C (0.2 cycles / pixel) (0.4 cycles / pixel) (0.8 cycles / pixel) 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 Figure 7.16: Noise Response of FBP Reconstructions vs. fc for MCAT Data Chapter 7. Results and Discussion 77 tn 5 0.45 0.4 0.35 0.3 0.25 0.2L Least RMS difference from original object vs f 0.2 0.4 0.6 0.1 f (cycles / pixel) Iteration No. vs f 0.2 0.4 0.6 0.8 1 fc (cycles / pixel) '0, 8 £ 6 I 4 2 0 Norm(V) vs f 0.2 0.4 0.6 0.S f (cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 120 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (0.9 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) : 20 40 60 80 100120 VA- : Variance image C (0.9 cycles / pixel) Co-Variance image C (0.9 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100120 Figure 7.17: Noise Response of 4-Subset OSEM Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 78 0) ffl £ s 0.2 \ 0.18 0.16 0.14 0.12. Reprojection misfit lo data closest to expected value vs f 0.2 0.4 0.6 0.8 I (cycles / pixel) .1 8 CO a> 7 6 Iteration No. vs t V B C \^ i... 0.4 0.6 0.8 1 f (cycles / pixel) Norm(V) vs f 0.2 0.4 0.6 0.8 t (cycles / pixel) Variance plot A (0.4 cycles / pixel) 20 40 60 80 100 120 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (0.9 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Variance image C (0.9 cycles / pixel) 250 200 150 100 50 Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (0.9 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100120 20 40 60 80 100120 Figure 7.18: Noise Response of 4-Subset OSEM Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 79 E g 0) 01 a s § 2 c .£ CO o rr 0.35 0.3 0.25 Least RMS difference from original object vs f c 0.5 1 1.5 (cycles / pixel) l °9 1 0 (a) vs f„ 0.5 1 1.5 f (cycles / pixel) 5000 4000 |-3000 c2000 1000 Norm(V) vs f 0.5 1 1.5 (cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 120 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) 40 30 20 10 0 Variance image B (0.6 cycles / pixel) 80 70 60 50 40 30 20 10 Variance image C (1.6 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) 20 40 60 80 100 120 20 40 60 80 100 120 Co-Variance image C (1.6 cycles / pixel) |40 120 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100120 Figure 7.19: Noise Response of FR-LS Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 80 Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) 40 30 20 10 0 Variance image C (1.6 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) 20 40 60^ 80 100 120 Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 20 40 60 80 100 120 1 0 0 : 20 40 60 80 100120 |10 20 40 5 60 80 0 100 120 10 0 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100120 Figure 7.20: Noise Response of RC-LS Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 81 E 0.35 o "i B 0) 0 c 2? obj 0.3 OJ to 3= c ra 0.25 cn o rr Least RMS difference from original object vs f 0.5 1 1.5 (cycles / pixel) log1Q(a) vs fc 0.5 1 1.5 r (cycles / pixel) 6000 |4000 2000 Norm(V) vs f c c 0.5 1 1.5 (cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 120 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Variance image C (1.6 cycles / pixel) 1 4 0 1 2 0 1 0 0 BO 60 40 20 20 40' 60 80 100 120 Co-Variance image A (0.4 cycles / pixel) 20 40 60 80 100 120 15 10 20 5 40 0 60 -5 80 -10 100 -15 120 1 Co-Variance image B (0.6 cycles / pixel) 20 40 60 80 100 120 mL 20 10 40 5 60 80 0 100 -5 120 Co-Variance image C (1.6 cycles / pixel) 20 40 60 80 100120 60 40 20 0 J-20 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100 120 Figure 7.21: Noise Response of FR-ML Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 82 c o s 3 <D CO •6-° cLS <u — ft Reprojection misfig to data closest to expected value vs f 0.16 0.155 0.15 0.145 0.14 0.5 1 1.5 f (cycles / pixel) 41 log10(<x) vs fc 0.5 1 1.5 (cycles / pixel) 5000 5-4000 £3000 2000 Norm(V) vs f 0.5 1 1.5 f (cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 120 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value 20 40 60 80 100 120 Variance image A (0.4 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100120 1 50 40 30 20 10 0 1 20 1 40 60 o 80 0 100 . -5 120 40 20 0 Variance image B (0.6 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100120 100 80 60 40 20 0 Variance image C (1.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 20 40 60 80 100 120 70 60 50 40 30 20 10 I 40 30 20 10 0 - 1 0 20 40 60 80 100 120 Figure 7.22: Noise Response of RC-ML Reconstructions vs. fc for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 83 e Is si i l II O) 0.82 0.8r 0.78 0.76 0.74 0.72 Least RMS difference from original object vs f 0.4 0.6 0.£ f (cycles / pixel) Iteration No. vs f V C 0.2 0.4 0.6 0.8 1 f. (cycles / pixel) 0.2 0.4 0.6 0.8 f (cycles / pixel) 2.6 t 2 a l i s Variance plot A (0.4 cycles / pixel) 50 100 150 Mean Pixel Value Variance plot B (0.6 cycles / pixel) 450 2500 400 2000 350 300 s • 11500 I250 > 5200 1000 + 150 500 • 100-50, 50 100 150 Mean Pixel Value Variance plot C (0.9 cycles / pixel) 50 1 00 150 Mean Pixel Value Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Variance image C (0.9 cycles / pixel) 350 300 250 200 150 100 50 Co-Variance image A (0.4 cycles / pixel) 20 40 60 80 100 120 V_/\ : 60 30 Co-Variance image B (0.6 cycles / pixel) 20 40 60 80 100 120 60 20 40 40 60 20 80 0 100 -20 120 100 50 0 Co-Variance image C (0.9 cycles / pixel) 20 40 60 80 100 120 100 150 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 Figure 7.23: Noise Response of 4-Subset OSEM Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 84 ssl 0.15 0.14 0.13 0.12 0.11 • 0.1 Reprojection misfit to data closest to expected value vs f ).4 0.6 0.8 1 (cycles / pixel) Iteration No. vs f 0.2 0.4 0.6 0.8 1 f (cycles / pixel) 12 10 £ 8 ¥ 6 o 8 4 2 0 Norm(V) vs f 0.2 0.4 0.6 0.8 f (cycles / pixel) Variance plot A (0.4 cycles / pixel) 50 100 Mean Pixel Value Variance plot B (0.6 cycles / pixel) Variance plot C (0.9 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Variance image C (0.9 cycles / pixel) 100 80 n o 40 20 Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (0.9 cycles / pixel) 20 40 60„ 80 100 120 s r 4 0c 20 40 60 80 100 20 40 60 80 100 120 -2 120 18 I'6 1 20 1 40 I, 60 80 0 100 120 20 40 60 80 100 120 -5 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 Figure 7.24: Noise Response of 4-Subset OSEM Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 85 E o * T3 gf 3= .g ^ .2> W o 5 rr 0.76 0.74 0.72 Least RMS difference from original object vs f C 0.5 1 1.5 (cycles / pixel) 3 4 cn 3 o log10(a) vs fc 0.5 1 1.5 (cycles / pixel) 4000 ~3000 §2000 c 1000 Noim(V) vs f 0.5 1 1.5 (cycles / pixel) Variance plot A (0.4 cycles / pixel) 50 100 150 Mean Pixel Value Variance plot B (0.6 cycles / pixel) 50 100 150 Mean Pixel Value Variance plot C (1.6 cycles / pixel) 50 100 150 Mean Pixel Value Variance image A (0.4 cycles / pixel) so 40 30 20 10 0 Variance image B (0.6 cycles / pixel) Variance image C (1.6 cycles / pixel) co 51) 40 .in 20 10 20 40 60 80 100 120 Co-Variance image A (0.4 cycles / pixel) 20 40 60 80 100 120 20 40 60 80 100 120 30 20 20 40 60 10 80 0 100 -10 120 40 20 n Co-Variance image B (0.6 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100120 140 20 20 40 60 80 0 100 120 40 20 0 Co-Variance image C (1.6 cycles / pixel) 20 40 60 80 100120 40 30 20 10 0 --10 20 40 60 80 100 120 Figure 7.25: Noise Response of FR-LS Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 86 Variance plot A Variance plot B Variance plot C (0.4 cycles / pixel) (0.6 cycles / pixel) (1.6 cycles / pixel) Figure 7.26: Noise Response of RC-LS Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 87 Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 20 40 60. 80 100 120 40 20 0 20 40 60 80 100120 • 60 |40 20 20 40 60 80 100 120 20 40 60 80 100120 J — 20 40 60 80 100120 20 40 60 80 100120 20 40 60 80 100120 20 40 60 80 100120 Figure 7.27: Noise Response of FR-ML Reconstructions vs. fc for MCAT Data, with Parameters Resulting in Least Difference from Original Object Chapter 7. Results and Discussion 88 Reprojection misfit to data closest to expected value vs f o •s a CD CO '5"° Q- S CD „ 0.115 0.11 0.105 0.1 0.5 1 1.5 fr (cycles / pixel) Variance plot A (0.4 cycles / pixel) 50 100 Mean Pixel Value 5 4 3 S 2 1 log10(a)vsfc Norm(V) vs f 0.5 1 1.5 f, (cycles / pixel) Variance plot B (0.6 cycles / pixel) 300 5-200 I100 0 0.5 1 1.5 (cycles / pixel) Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) Variance image B (0.6 cycles / pixel) Variance image C (1.6 cycles / pixel) 20 40 60 80 100 120 Co-Variance image A (0.4 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100 120 Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 20 40 60 80 100120 20 40 60 80 100120 ±Z2 2 1.5 1 0.5 0 -0.5 20 40 60 80 100120 Figure 7.28: Noise Response of RC-ML Reconstructions vs. / c for MCAT Data, with Parameters Resulting in Expected Misfit to Data Chapter 7. Results and Discussion 7.3 SUMMARY 89 Results Features Fig 7.2 Panels A & B and other Figs Selecting reconstruction parameters that result in least difference from the original object results in noisier and higher resolution reconstruction that parameter selection aimed at producing a data misfit as close to the expected value as possible Figs. 7.1 and 7.2 Streaking artifacts typical of FBP reconstructions, especially in low count areas of the image Figs. 7.4 and 7.7 Slight streaking around low intensity edges due to the FR-LS least squares fit term Figs. 7.3 and 7.5 Different noise textures between OS-EM and FRCT reconstructions, especially in the low intensity extended source Figs. 7.3 and 7.5 7.12, 7.14 FR-ML achieves better contrast than OSEM in cold spots Figs. 7.15 and 7.16 Continuous relationship between the norm of pixel variances and FBP cut-off frequency Fig. 7.23 Highly unstable nature of OSEM at high numbers of iterations. Reconstructing more random noise realizations may yield more enterpretable results Figs. 7.17, 7.18, 7.23 and 7.24 Response of OSEM to noise in the acquired data varies eratically with respect to reconstruction parameters, the norm of the pixel variances appears to be discontinuous in "Norm(V) vs. / c " plots. Reconstructing more random noise realizations and using finer resolution in parameter sweeps may yield more enterpretable results Table 7.4: Summary of Results - Part 1 Chapter 7. Results and Discussion 90 Results Features Table 7.3 Parameter settings which attempt to generate expected misfit to data result in reconstructions which are more stable with respect to noise. Table 7.3 FRCT reconstructions have consistently lower pixel variances than corresponding OSEM and FBP reconstructions Figs. 7.17, 7.18, 7.23, and 7.24 Linear relationship between the mean pixel value of OSEM reconstructions and corresponding variance in scatter plots Figs. 7.19-7.22, and 7.25-7.28 Stabilizing effect of Fourier Regularozation can be seen in the shape of the scatter plots Figs. 7.1 and 7.2 Co-variance images show the relationship between correlation lengths in the reconstructed image and FBP cut-off frequency Figs. 7.26 and 7.28 Co-variance images show larger correlation lengths than any other combination of reconstruction method, phantom and method of parameter selection Table 7.5: Summary of Results - Part 2 C H A P T E R 8 C O N C L U S I O N In this work, a new method for SPECT reconstruction has been presented. The method successfully applies theory originally developed in the field of Astrophysics to static SPECT reconstruction. 8.1 G E N E R A L C O M M E N T S The primary advantage of the FRCT reconstruction method is in it's response to noise in the reconstructed data. All reconstruction methods make some form of compromise between image resolution and stability in the reconstructed image. It has been shown that the FRCT reconstruction method achieves a better trade-off between image resolution and stability with respect to noise than the other reconstruction methods studied in this work. Future work might include clinical trial of the FRCT reconstruction method, and further development of the method and the model of the SPECT acquisition system, as described below. 91 Chapter 8. Conclusion 92 8.1.1 F U T U R E W O R K F L U X CONSERVATION T E R M In the field of Baysian image processing, adding a "flux conservation" term to the maxi-mum entropy objective function has been proposed [31]. This method could be applied to FRCT reconstruction, in order to associate physical units such as c™™ta with the reconstructed measurements. T H E S Y S T E M M A T R I X Preliminary results obtained using a system matrix which models fully three-dimensional collimator blurring show that it has a strong stabilizing influence on any iterative recon-struction method. Because the acceptance angle of the SPECT collimator is modeled in this approach, more sinogram elements contribute to each pixel in the reconstructed image, which leads to reconstructions which are less sensitive to statistical variations in the acquired data. Once this model incorporates attenuation effects, it may significantly improve reconstructions of clinical SPECT data. R E G U L A R I Z A T I O N P A R A M E T E R The regularization parameter, «, is a nuissance from the point of view of the person carrying out an FRCT reconstruction. Other fields have employed similar data analysis techniques, and have developed methods for selecting the "natural" or optimal regular-ization parameter for each situation. Examples include Skilling and Gull's work with the Maximum Entropy regularization parameter [32], [33], and the method described in [34]. Ideally, the use of any reconstruction method should be as straight-foward as possible from the user's point of view. One possibility for FRCT would be to allow the user to select a level of stability, and have the algorithm select the highest resolution image that Chapter 8. Conclusion 93 can be reconstructed with the specified stability. D E T E C T O R G E O M E T R I E S A N D IMAGING M E T H O D S FRCT could also be adapted to detector geometries and imaging modalities other than parallel hole collimator SPECT. The application of FRCT to PET would require a weighted least-squares fit term due to the method's Gaussian noise [22] and an appro-priate system matrix. The application of FRCT to electronically collimated "Compton Camera" systems would require an interesting system matrix due to the shape of the back-projection cone. The use of FRCT with list mode data would require a different method of modelling the acquisition system. Such a model might be expressed in the form of a function instead of a transformation matrix. B I B L I O G R A P H Y G. Hounsfield. Computerized transverse axial scanning (tomography): Part I. De-scription of system. Brit. J. Rad., 47:1016-1022, 1973. Gopal B. Saha. Physics and Radiobiology of Nuclear Medicine. Springer-Verlag, New York, 1993. J.A. Sorenson and M.E. Phelps. Physics in Nuclear Medicine. W.B. Sauders Com-pany, 1987. Anna Celler. Private communication. April 2000. J. Radon. Uber die bestmmung von funktionen durch ihre integralwerte langs gewisser mannigfaltigkeiten. Berichte Sachsische Akademie der Wissenschaften. Leipzig, Math. - Phys. KL, 69:262-267, 1917. Shtein. On the application of the radon transform in holographic interferometry. Radiotekh. Elektron., 17:2436-2437, 1972. B. K. Vainshtein and S. S. Orlov. Theory of the recovery of functions from their projections. Kristallografiya, 17:253-257, 1972. C M . Vest. Application of radon transforms to multi-directional interferometry. J. Opt. Soc. Am., 64:468, 1973. A. M. Cormack. Reconstruction of densities from their projections, with applications in radiological physics. Phys. Med. Biol., 18:195-207, 1973. T.F. Budinger, S. E. Derenzo, G.T. Gullberg, W. L. Greenberg, and R. Ft. Huesman. Emission computed assisted tomography with single-photon and positron annihila-tion photon emitters. J. Comput. Assisted Tomog., 1:131-145, 1977. M.E. Phelps. Emission computed tomography. Semin. Nucl. Med., 7:337-365, 1977. G.L. Brownell and S. Cochavi. Transverse section imaging with carbon-11 labeled carbon monoxide. J. Comput. Assisted Tomog., 2:533-538, 1978. T.F. Budinger, G.T. Gullberg, and R.H. Huesman. Emission computed tomography Image Reconstruction From Projections. Springier-Verlag, New York, 1979. S.R. Deans. The Radon transform and some of its applications. John Wiley & Sons, Toronto, 1983. 94 Bibliography 95 [15] A. Lannes, E. Anterrieu, and K. Bouyoucef. Fourier interpolation and reconstruction via shannon-type techniques; part 1: Regularization principle. J. Mod. Opt., 41, 1994. Stephan Blinder. Private communication. January 2000. G. T. Gullberg. The attenuated radon transform: Application to single-photon emis-sion computed tomography in the presence of a variable attenuating medium. Tech-nical report, Lawrence Berkeley Laboratory, Berkeley, CA, 1980. Lee W. Johnson, R. Dean Riess, and Jimmy T. Arnold. Introduction to Linear Algebra. Addison-Wesley, Reading Massachusetts, third edition, 1983. E. Kreyszig. Advanced Engineering Mathematics. John Wiley & Sons, 1993. F. Natterer. The Mathematics of Computerized Tomography. Wiley & Sons, New York, 1986. A.G. Ramm and A.I. Katsevich. The Radon Transform and Local Tomography. CRC Press, Boca Raton, 1996. J.A. Fessler. Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Trans. Med. Imag., 13(2) :290-300, 1994. L.A. Shepp and Y. Vardi. Maximum liklihood reconstruction for emission tomogra-phy. IEEE Trans. Med. Imag., MI-L113-122, 1982. Philip R. Bevington and D. Keith Robinson. Data Reduction and Error Analysis for the Physical Sciences. McGraw Hill, New York, 1992. H. M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag., 13:601-609, 1994. P. Marechal, D. Togane, A. Celler, and J. Borwein. Numerical assessment of the stability of reconstruction processes for computed tomography. IEEE Trans. Nuc. Sci., page 2177, December 1999. Kaplan M.S., Harrison R.L., and Vannoy S.D. Coherent scatter implementation for simset. IEEE Trans. Nucl. Sci., 45(6):3064-3068, 1998. P.H. Pretorius, W. Xia, M.A. King, B.M.W. Tsui, T.S. Pan, and B.J. Villegas. Eval-uation of right and left ventricular volume and ejection fraction using a mathematical cardiac torso phantom for gated pool spect. J. Nuc. Med., 38(10):1528-1534, 1997. Bibliography 96 [29] D.W. Wilson and B.M.W Tsui. Noise properties of filtered-backprojection and ML-EM reconstructed emission tomographic images. IEEE Trans. Nuc. Sci., 40(4):1198-1203, 1993. [30] J.S. Beis, A. Celler, and J.S. Barney. An automated method to determine cutoff frequency based on image power spectrum. IEEE Trans. Nuc. Sci., 42(6):2250-2254, 1995. [31] Nailong Wu. Maximum entropy task for iraf / stsdas. Image Restoration Newsletter, (1):72, 1993. [32] J. Skilling. Classic maximum entropy. Maximum Entropy and Baysian Methods, 45, 1989. [33] S.F. Gull. Developments in maximum entropy data analysis. Maximum Entropy and Baysian Methods, 53, 1989. [34] Simon S. Haykin. Neural networks : a comprehensive foundation. Prentice-Hall, second edition, 1999. C H A P T E R 9 A P P E N D I X A 97 Chapter 9. Appendix A gg High Cutoff Frequency (0.9 cycles / pixel) High Cutoff Frequency (0.9 cycles / pixel) High Iteration No. (8) Low Iteration No. (3) Low Cutoff Frequency (0.5 cycles / pixel) Low Iteration No. (8) Low Cutoff Frequency (0.5 cycles / pixel) High Iteration No. (14) f in cycles / pixel iterations Figure 9.1: Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital Circle Phantom vs fc and Iteration No. - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 99 Low Cutoff Frequency (0.6 cycles / pixel) Low a (10A1) Low Cutoff Frequency (0.6 cycles / pixel) High a (10A4) High Cutoff Frequency (1.2 cycles / pixel) High a (10*5.5) iog 1 0(a) f in cycles / pixel Figure 9.2: Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital Circle Phantom vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 100 Low Cutoff Frequency (0.6 cycles / pixel) Low a (10A-1) Low Cutoff Frequency (0.6 cycles / pixel) High a (10A1.5) 1 IMA \n High Cutoff Frequency (1.6 cycles / pixel) High a (10A3) log10(a) 3 5 4 i s f in cycles / pixel Figure 9.3: Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital Circle Phantom vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 101 Low Cutoff Frequency (0.5 cycles / pixel) Low Cutoff Frequency (0.5 cycles / pixel) High Cutoff Frequency (0 9 cycles / pixel) High Iteration No. (9) Low Iteration No. (5) Low Iteration No. (4) Figure 9.4: Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs fc and Iteration No. - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 102 Low Cutoff Frequency (0.6 cycles / pixel) Low Cutoff Frequency (0.6 cycles / pixel) High Cutoff Frequency (1 2 cycles / Low a (10A1) High a (10*5) High a (10*6.5) Figure 9.5: Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital MCAT Phantom vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 103 Low Cutoff Frequency (0.6 cycles / pixel) Low Cutoff Frequency (0.6 cycles / pixel) High Cutoff Frequency (1 2 cycles / Low a (10A1) High a (10A2.5) High a (10A4) Figure 9.6: Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital MCAT Phantom vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 104 High Cutoff Frequency (0.9 cycles / pixel) High Iteration No. (7) High Cutoff Frequency (0.9 cycles / pixel) Low Iteration No. (3) Low Cutoff Frequency (0.5 cycles / pixel) Low Iteration No. (7) iterations Figure 9.7: Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 105 Low Cutoff Frequency (0.6 cycles / pixel) Low Cutoff Frequency (0.6 cycles / pixel) High Cutoff Frequency (1.6 cycles / Dixel) I rt..F 11AA1 C \ I T : _ I_ / 1 r\ A < r \ .... _ . _ £. • ' Low a ( 0A .5) High a ( 0A4.5) High Cutoff Frequency (1.6 cycles / pixel) Low a (10A2.5) iii n l High a (10A5.5) iog 1 0(a) e f in cycles / pixel Figure 9.8: Surface Plot of RMS Difference Between Reprojection of FR-LS Reconstruc-tions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 106 Figure 9.9: Surface Plot of RMS Difference Between Reprojection of FR-ML Recon-structions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 107 Figure 9.10: Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 108 Figure 9.11: Surface Plot of RMS Difference Between Reprojection of FR-LS Recon-structions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings Chapter 9. Appendix A 109 Figure 9.12: Surface Plot of RMS Difference Between Reprojection of FR-ML Recon-structions and Sinogram Data vs fc and a - Sample reconstructions show effects of faulty parameter settings 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085173/manifest

Comment

Related Items