FOURIER REGULARIZATION APPLIED TO SPECT RECONSTRUCTION By Dylan Dalmar Togane B. Sc. Concordia University, 1997 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E O F M A S T E R OF SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES D E P A R T M E N T O F PHYSICS A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UNIVERSITY O F 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 reconstruction 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 experiments 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 ii Table of Contents iv List of Tables vii List of Figures viii 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 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 1.5.5 Detector Efficiency 1.6 2 Applications : . 11 12 Emission Computed Tomography Image Reconstruction 12 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 3 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 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 4 5 6 Iterative methods 34 3.3.1 34 ML-EM and OS-EM Regularized S P E C T Reconstruction 37 4.1 Classical Regularizers 38 4.2 Fourier Regularization 38 Sensitivity to Noise 41 5.1 FBP 42 5.2 Iterative Methods 43 Reconstructions 45 6.1 Simulation of SPECT data 45 6.1.1 Analytic Simulations 45 6.1.2 Monte-Carlo Methods 6.2 46 Digital Phantoms 46 6.2.1 46 Analytic Circles v J 6.2.2 6.3 6.4 Monte-Carlo MCAT Phantom 47 Reconstruction Evaluation Criteria 48 6.3.1 Difference from Original Object and Re-projection Misfit to Data 48 6.3.2 Noise Response 50 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 8.1 , 91 General Comments 91 8.1.1 92 Future Work Bibliography 94 9 Appendix A 97 vi L I S T 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. / " Plots 72 7.4 Summary of Results - Part 1 89 7.5 Summary of Results - Part 2 90 c vii L I S T OF F I G U R E S 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 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 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 . 17 , . 23 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 Reconstructions and Analytic Circle Phantom vs f (bottom, left) and RMS c Difference Between Re-projection of FBP Reconstructions and Analytic Circle Sinogram vs f (bottom, right) 57 c 7.2 Selected FBP Reconstructions (top), RMS Difference Between FBP Reconstructions and MCAT Digital Phantom vs / c (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and MCAT Monte-Carlo Sinogram vs f (bottom, right) 58 c 7.3 Least RMS Difference Between 4 Subset OSEM Reconstructions and Digital Circle Phantom vs f (left) and Corresponding Iteration No. vs f c c (right) 7.4 59 Least RMS Difference Between FR-LS Reconstructions and Digital Circle Phantom vs f (left) and Corresponding a vs f (right) c 7.5 60 c Least RMS Difference Between FR-ML Reconstructions and Digital Circle Phantom vs f (left) and Corresponding a vs f (right) c 7.6 61 c Least RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs f c (left) and Corresponding Iteration No. vs f c (right) . 7.7 62 Least RMS Difference Between FR-LS Reconstructions and Digital MCAT Phantom vs f (left) and Corresponding a vs f (right) c 7.8 c Least RMS Difference Between FR-ML Reconstructions and Digital MCAT Phantom vs f (left) and Corresponding a vs f (right) c 7.9 63 c 64 Least RMS Difference Between Re-projection of 4 Subset OSEM Reconstructions and Sinogram Data vs f c No. vs f (right) (left) and Corresponding Iteration 65 c ix 7.10 Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs f (left) and Corresponding a vs f c (right) c . . . . 66 7.11 Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs f (left) and Corresponding a vs f c (right) c . . . . 67 7.12 Least RMS Difference Between Re-projection of 4 Subset OSEM Reconstructions and Sinogram Data -vs f c (left) and Corresponding Iteration No. vs f (right) 68 e 7.13 Least RMS Difference Between Re-projection of FR-LS Reconstructions and Sinogram Data vs f (left) and Corresponding a vs f (right) c c . . . . 69 7.14 Least RMS Difference Between Re-projection of FR-ML Reconstructions and Sinogram Data vs f (left) and Corresponding a vs f . . . . 70 7.15 Noise Response of FBP Reconstructions vs. f for Analytic Circle Data . 75 7.16 Noise Response of FBP Reconstructions vs. f for MCAT Data 76 c (right) c c c 7.17 Noise Response of 4-Subset OSEM Reconstructions vs. f c for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object 77 7.18 Noise Response of 4-Subset OSEM Reconstructions vs. f c for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data . . . 78 7.19 Noise Response of FR-LS Reconstructions vs. f for Analytic Circle Data, c with Parameters Resulting in Least Difference from Original Object . . . 79 7.20 Noise Response of RC-LS Reconstructions vs. f for Analytic Circle Data, c with Parameters Resulting in Expected Misfit to Data 80 7.21 Noise Response of FR-ML Reconstructions vs. f for Analytic Circle Data, c with Parameters Resulting in Least Difference from Original Object . . . x 81 7.22 Noise Response of RC-ML Reconstructions vs. f for Analytic Circle Data, c with Parameters Resulting in Expected Misfit to Data 82 7.23 Noise Response of 4-Subset OSEM Reconstructions vs. f for MCAT Data, c with Parameters Resulting in Least Difference from Original Object . . . 83 7.24 Noise Response of 4-Subset OSEM Reconstructions vs. f for MCAT Data, c with Parameters Resulting in Expected Misfit to Data 84 7.25 Noise Response of FR-LS Reconstructions vs. f for MCAT Data, with c Parameters Resulting in Least Difference from Original Object 85 7.26 Noise Response of RC-LS Reconstructions vs. f for MCAT Data, with c Parameters Resulting in Expected Misfit to Data 86 7.27 Noise Response of FR-ML Reconstructions vs. f for MCAT Data, with c Parameters Resulting in Least Difference from Original Object 87 7.28 Noise Response of RC-ML Reconstructions vs. f for MCAT Data, with c Parameters Resulting in Expected Misfit to Data 9.1 Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital Circle Phantom vs f and Iteration No c 9.2 99 c Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital Circle Phantom vs f and a 100 c 9.4 Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs f and Iteration No c 9.5 98 Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital Circle Phantom vs f and a 9.3 88 101 Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital MCAT Phantom vs f and a 102 c xi 9.6 Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital MCAT Phantom vs f and a 103 c 9.7 Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs f and a c 9.8 Surface Plot of RMS Difference Between Reprojection of FR-LS Reconstructions and Sinogram Data vs f and a 105 c 9.9 104 Surface Plot of RMS Difference Between Reprojection of FR-ML Reconstructions and Sinogram Data vs f and a 106 c 9.10 Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs f and a c 107 9.11 Surface Plot of RMS Difference Between Reprojection of FR-LS Reconstructions and Sinogram Data vs f and a c . . : 108 9.12 Surface Plot of RMS Difference Between Reprojection of FR-ML Reconstructions and Sinogram Data vs f and a c xii 109 CHAPTER 1 INTRODUCTION 1.1 A I M OF THIS W O R K Single Photon Emission Computed Tomography (SPECT) is a diagnostic medical imaging 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 F R E C T and other SPECT reconstruction methods on SPECT data, allowing comparison. 1.2 STRUCTURE OF T H E 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 regularization is intended to complement the existing techniques. Chapter 2 presents the Fourier Slice Theorem and its relevance to E C T 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, collimator blurring correction, and scatter correction are outlined with varying degrees of detail. A variety of formulations of the image reconstruction problem are presented in Chapter 3. The principles behind some existing reconstruction methods (FBP, ML-EM, OSEM) are described. The general formulation of image reconstruction as an inverse problem is also presented here. In Chapter 4, Fourier regularization for E C T reconstruction is presented, as well as three classical regularization terms. The Least Squares and Maximum Likelihood fit functions, 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^are also discussed. Additional results are presented without discussion in Chapter 9. 1.3 N U C L E A R MEDICINE IMAGING Diagnostic medical imaging began in 1895 when Wilhelm Roentgen produced X-Ray images 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), single 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 201 T1 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 D E C A Y AND INTERACTIONS OF PHOTONS WITH MATTER Radioactive decay is a process in which an unstable [parent] nucleus transforms into a more stable [daughter] nucleus by emitting particles and/or photons 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 summarized 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 radioactive 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 P E T 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 A N N I H I L A T I O N 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 = I e-* "M* ', 1 Q (1.2) where I is the initial beam intensity, and the integration is carried out over the path 0 that the radiation travels. [3] 1.5 T H E A N G E R 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-POSITION SIGNAL CATHODE RAY T U B E Y-POSITION SIGNAL GATING CIRCUITRY PULSE-HEIGHT ANALYZER Z-PULSE POSITION LOGIC CIRCUITS Q P M - T U B E ARRAY NoI(TJ0CRYSTAL COLLIMATOR PATIENT 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 C R Y S T A L 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 P M T s , E L E C T R O N I C S 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 f'^'. [4] c 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 Organ Brain and Spinal Cord 13 Radio-pharmaceutical Tc-Sodium Pertecnetate ( 99m 99m T c . D T p 99m Tc04 ) A I-Isopropyl-p-iodoamphetamine (IMP) 123 Tc-Hexamethylpropylene Amine Oxime (HMPAO) Tc-Ethyl Cyteinate Dimer (ECD) F-Flouordeoxyglucose (FDG) 99m 99m ls 15 0 - H 0 , Xe In-DTPA 13& 2 ni Thyroid I-Sodium Iodide (and other iodine isotopes) Tc-Sodium Pertechnetate Tc-Macroaggregated Albumin ( Tc-MAA) X e and Xe gas Tc-Aerosol 131 99m Lung 99m Diagnosis blood flow, tumor, brain death cerebral perfusion and metabolism cerebral perfusion (stroke, infarction, tumor) cerebral perfusion glucose metabolism oxygen consumption blood flow leakage of cerebro-spinal fluid (CSF), hydrocephalus thyroid function and structure morphology perfusion defects 99m i33 12Y 99m Liver 1-Rose Bengal 131 T c - I D A derivatives (-HIDA, -PIPIDA, -BIDA, -DISIDA) Tc-Sulfur Colloid and Tc-Albumin Colloid Ga-gallium Tc-Sulfur Colloid and Tc-Albumin Colloid Tc-labeled Red Blood Cells Se-Selenomethionine I- or I-Orthoiodohippurate (Hippuran) Tc-Mercaptoacetylglycylglycine (MAG3) Tc-DTPA Tc-Gluceptate Tc-DMSA 99m 99m 99m e7 Spleen 99m 99m 99m Pancreas Kidney 75 123 99m 99m 99m airway obstructions abnormal distribution in bronchial diseases liver disease, biliary obstruction gallbladder and biliary tree imaging tumor detection abscesses and tumors morphology morphology digestive enzymes 131 filtration rate and effective renal plasma flow (EPRF) 99m Table 1.1: Nuclear Medicine Radio-pharmaceuticals and Applications - part 1 [2] Chapter 1. Introduction Organ (cont'd) Skeleton 14 Radio-pharmaceutical (cont'd) Tc-Phosphate and Phosphonate complexes ( Tc-MDP and Tc-HDP) T1-Thallous Chloride Tc-Sestamibi Tc-Teboroxime Rb-Rubidium Chloride N-Ammonia F-Fluorodeoxyglucose (FDG) Tc-Pyrophosphate l n - and. I-labeled Antimyosin Antibody Tc-RBC 99m 99m Heart 99m Diagnosis (cont'd) bone blood flow, and bone formation rate 201 yym perfusion studies yym S2 13 18 99m U 1 PET perfusion studies PET metabolic studies myocardial infarct 131 yym gated equilibrium cardiac blood pool Table 1.2: Nuclear Medicine Radio-pharmaceuticals and Applications - part 2 [2] Name Nuclear Notation 2Ay A !Ll ->• Vr X v* X A v A /3~ Emission and z + l Ay (/3~,7) Emission A Y ?1 A v z Z-l X £ i^Y* 4 _ Y Positron (/?+) A and (/3 ,7) decay + A Z Z Electron Capture (EC) and (EC,7) decay ix iX Decay by a emission Y X 4 i_ v Description The e~ is called a beta particle (/?"). The daughter nucleus is in an excited state and then emits a photon. The daughter nucleus is sometimes in an excited state [(P ,j) decay]. + An orbital e~~ is captured by the nucleus, transforming a proton into a neutron. E x 4 i^Y* 4 i_ Y E x The daughter nucleus is sometimes in an excited state [(EC/y) decay]. The a-particle consists of two neutrons and two protons (a ^He nucleus). Table 1.3: Types of Radioactive Decay [3] CHAPTER 2 T H E FOURIER SLICE T H E O R E M AND COMPUTED TOMOGRAPHY 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- < ' >dr, i27r / k r (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 i27r<k ' dk. r) (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 -fieldswhere it has been applied include Astronomy, 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 9 - and passing at a distance p from the origin (L in Figure 2.1). s 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 , the / [(the Radon 2 transform of /)] is determined by the line integrals of / . And if / is defined on R then / is determined by the surface integrals of / over two-dimensional 3 planes. [14] Chapter 2. The Fourier Slice Theorem and Computed Tomography 2.2 19 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 BackProjection 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 m )e- ^dp, i2 P / (2.4) •inf where k = k cos(9) and k — k sin(9). x p y p 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,k ) is known on a star-shaped domain p 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-logi plot of f(k ) when the Radon 0 p 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-logi plot of 0 f{k ) when no noise is present, and the Fourier transformed random noise vector. p 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 uniformly 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 magnitude 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 thefidelitywith 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 where I max and I min = TJ—^7~T' ( 5 ) 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) = C (k)/C (k), out in (2.6) where the test pattern shown in Figure 2.6 is used to determine the M T F 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 TRANSFORM AND T H E SYSTEM MATRIX 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 0 50 100 150 200 25 250 nz = 4904 Figure 2.8: Sparsity Structure of the Radon Matrix 2.6 COLLIMATOR 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 x +z 2 PSF(x, z) = exp 2 2(v + * Ay))\ 2 2 ( 2 ' 7 ) where cr is the intrinsic response of the scintillation crystal, PMTs and electronics, and t <T is the depth dependent collimator response. [16] c detector source Figure 2.9: Collimator Blurring Geometry 2.7 ATTENUATION 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, ,e) = exp [ - J ^ ^ M 0 * ( p - < * , O ) # ] P r U(e,p) = [nMe, ) = fA(tP,o)f{£)6(jp-(o,z))dt. P where 0 L is a unit vector perpendicular to 0, pointing along L in Figure 2.10, and £' • 0 > £ • 6 denotes the segment of L in between the point £ and the physical camera 1 L 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. CHAPTER 3 SPECT 3.1 FILTERED RECONSTRUCTION BACK-PROJECTION METHODS (FBP) 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 without 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 fC> (B) (A) 0.02 0.015 0.01 0.005 20 40 60 80 100 120 20 40 60 80 100 120 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 F B P 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 It can be shown ([19]) = ju(9,p)v(d,p-q)dq. (3.4) that a convolution in the spatial domain corresponds to a multiplication in the Fourier domain: inf U(k)V(k)e * dk, i2 / •inf {Kr) (3.5) Chapter 3. SPECT Reconstruction Methods where $ _ 1 30 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 transforms of the convolution kernel and acquired data with respect to the radial variable p, and <3>~ represents the inverse Fourier transform with respect to the radial variable p. 1 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 experimental 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: •(K*g)*f = n*(g*Hf) (3.7) where * is a normal 2-dimensional convolution (as opposed to *, the convolution with respect to p in the Radon domain). 3.2 R E C O N S T R U C T I O N AS A N 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. x = argmin{o(x) = e(i?x,y) x G R ,x G n (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 W which may contain positivity constraints), and e is a 1 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 AND 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 L norm of the difference between Rx and y: 2 e (Rx, ) = || Jfcc - y ls y || = ]T([ifcc]i - yi) , 2 (3.9) 2 i also known as the squared Euclidean norm or the least squares function. The gradient (needed by quasi-Newtonian minimization algorithms) of ei with respect s to x is: V e {Rx, y) = 2R* \\ [Rx] - y ||. x ls (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-AO'E-^y-ifc) (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 likelihoodfitterm 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 LIKELIHOOD 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: y = [Rx\i = V \ Ri,j j x t 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]: lnP(y,p) = Y^iVilnPi] - ^ P i ~ J^Ms/iO- (3.13) 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 ; : e i(Rx,y) m = -lnP(y,Rx) - J^fn^!) = ]P[[-Rx]; - yiln[Rx\i . (3.14) Taking the gradient of this function with respect to x (for use with quasi-Newtonian algorithms) gives us [23]: (3.15) 3.3 I T E R A T I V E 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 AND 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 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 dominates 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 O S - E M 8 subsets, 1 iteration O S - E M 8 subsets, 4 iterations 36 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. CHAPTER 4 REGULARIZED SPECT RECONSTRUCTION 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 = argmin |o(x) = ^e(Rx,y) + ^ap(x) x € W >, 1 (4.1) 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 38 REGULARIZERS 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 p (x) bs = ^Xjlnxj, p (x) kl = y^Xjln^-, j=i °i j=i Pii(x) = In the definition of p^i, the ||x|| , 2 XQJ'S • x p (x) q = (4.2) (x,Qx). are the components of some reference model of the object. In the definition of p , Q is a symmetric positive semi-definite matrix, which q 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 minimization algorithms) are [26]: V p (x) x b5 = lnx, V p i(x) x fc = lnx-lnxo, (4.3) V /3tt(x) x = 2x, V p,(x) x = 2Qx, with the In function understood to be component-wise (lnx = [Inxi, lnx , lnx 2 4.2 lnx j]). 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 sampled at all. The Fourier domain is partitioned (see Figure 4.1) in order to ensure that Chapter 4. Regularized SPECT Reconstruction •r^— -iq 1 i fi 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 optimization problem: / ss argmin j o(f) = -e{Tlf, s * /) + 11| (1 - S)F \ l f = f(v)eV,TER s }. (4.5) In a discrete setting, this becomes x = argmin |o(x) = ^e(Rx, C y) + | | | B x | | s s 2 X € K , x e 7> j , n (4.6) where C is the matrix representation of the s * • operation and B is the matrix reps s 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. l|2 II The Fourier regularizer p/ = B x || is a special case of the quadratic regularizer p s q with Q = B*B . The gradient of pj is therefore: S V x P / = 2B*B x. s (4.7) CHAPTER 5 SENSITIVITY TO NOISE In this chapter the effect, on the reconstructed image, of experimental noise in the acquired data will be analyzed for a selection of reconstruction methods. For some methods the sensitivity to noise can be studied through mathematical analysis, 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) = 0). x Using sensitivity analysis, a quantitative measure of a reconstruction method's sensitivity 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 response to noise (such as the covariance matrix) is desired. All methods can have their response to noise studied empirically by repeatedly reconstructing 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 C calculated in Chapters 6 and 7 is related to the x sensitivity matricies described below by: C = SCyS* , (5-1) X where C is the sinogram co-variance matrix. [26] y 5.1 FBP The discrete FBP reconstruction can be described as x = R*&- g* y (5.2) 1 p where R*, ^ p , g, and <fr are matrix representations of the elements making up the 1 p 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 = S 6y, fbp with S = iii**- ^. 1 fbp (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 5.2 43 ITERATIVE METHODS 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 ) = ^-V e(R5t, y) + | V p(x) = 0. x x (5.5) x Replacing y with y + Sy, the same condition must hold with x replaced by the perturbed solution x + Sx: Y V e(i?[x + Sx],y + Sy) + | V p(x + 8x) = 0. x (5.6) x Keeping only the first order terms in the development of Equation (5.6) yields (using Equation (5.5)) [R*H R + aH ] Sx = -R*H Sy, £ e (5.7) £ i in which we have denoted H and H the Hessian (matrix of second partial derivatives) £ e of £ at (Rx,y) and the Hessian of g at x, respectively. We therefore obtain Sx = SSy with S = -[R+HcR + aHg]' R*H . 1 £ (5.8) The sensitivity matrix S depends on the Hessians of e and p, listed in the table below: Chapter 5. Sensitivity to Noise 2i?*|i?x-y|| V e (Rx,y) x 44 ls V ;f /(Z?x,y) = £^[1-(*/[«*],)] V p (x) = lnx V pfci(x) = ln(x/x ) VxPti(x) = x m x bs x 0 ' 2x V Pg(x) 2Qx V P (X) 2B*B x x X S e, H m = 2R*R = Vi/x] H Pbs = 1/x H p k l = x /x H P U 0 2 2Q 8 Table 5.1: Gradients and Hessians of e and p = 2/3;^ CHAPTER 6 RECONSTRUCTIONS 6.1 SIMULATION O F 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 MONTE-CARLO METHODS 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 CIRCLES 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 MONTE-CARLO MCAT PHANTOM 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 (without 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 projection 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 RECONSTRUCTION EVALUATION CRITERIA In order to quantitatively compare reconstructions, the two criteria described below were used. However, the final judgment is reached by qualitative comparison. Certain characteristics 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 candidate images are then compared and discussed. 6.3.1 DIFFERENCE FROM ORIGINAL O B J E C T AND RE-PROJECTION MISFIT TO DATA The room mean square (RMS) difference between the original object and the reconstructed 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 norms are one: 2 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 characteristic 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||, i=i where 6.3.2 (6.4) is the i h noise vector. t NOISE RESPONSE 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 i th pixel value in the n th image from a sample of N noisy images, and n=l Cjj then refers to the co-variance of the i th and Cij would be the variance of the i th and j th pixels in the reconstructed image, 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 6.4 51 RECONSTRUCTION PARAMETERS 6.4.1 SELECTING CUTOFF FREQUENCY AND FILTER 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: (6.7) otherwise where / is the input frequency and f is the cut-off frequency. c In practice, knowledge of the imaging system's MTF and the expected nose level 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 cutoff frequency have been investigated (see [30]). In this work, the cut-off frequency is considered to be a free parameter. 6.4.2 REGULARIZATION PARAMETERS AND STOPPING CRITERIA For the O S - E M reconstructions, the number of iterations is treated as a free parameter. In O S - E M reconstructions, the number of iterations is used (i) as the stopping criteriurri for the iterative algorithm, and (ii) as a regularization technique, mitigating the highfrequency artifacts that appear at higher numbers of iterations. A Hannfilteris 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 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 PARAMETER RANGES 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 parameters. For the FBP reconstructions this results in two curves, both plotted versus f , for c 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 FBP OSEM (four subsets) 53 Analytic Circle Phantom f = 0.1 ... 1 cycles / pixel in steps of 0.1 1 ... 15 iterations in steps of 1 f = 0.1 ... 1 cycles / pixel in steps of 0.1 a = lO^.-lO on a log scale in steps of IO f = 0.2 ... 1.8 cycles / pixel in steps of 0.2 a = 10 ...10 on a log scale in steps of IO f = 0.2 ... 1.8 cycles / pixel in steps of 0.2 c c RC-LS 55 05 c RC-ML _1 4 05 c M C A T Phantom f = 0.1 ... 1 cycles / pixel in steps of 0.1 1 ... 11 iterations in steps of 1 f = 0.1 ... 1 cycles / pixel in steps of 0.1 a = 1...10 on a log scale in steps of IO f = 0.2 ... 1.8 cycles / pixel in steps of 0.2 a — 10~ ...10 on a log scale in steps of IO f — 0.2 ... 1.8 cycles / pixel in steps of 0.2 c c 7 05 c 15 5 05 c 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. variance image is plotted versus f . c The norm of the 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 f . c 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. CHAPTER 7 RESULTS A N D DISCUSSION Tables 7.1 and 7.2 show thefigurenumbers for each combination of phantom, reconstruction 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 Section 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 Section 7.3. FBP FBP Circles MCAT Difference from Original Object Fig. 7.1 Fig. 7.2 Re-projection Misfit from Data Fig. 7.1 Fig. 7.2 Table 7.1: FBP Results 54 Noise Response Fig. 7.15 Fig. 7.16 Chapter 7. Results and Discussion Least Difference from Original Object OS-EM FR-LS FR-ML OS-EM FR-LS FR-ML Circles Circles Circles MCAT MCAT MCAT Fig. Fig. Fig. Fig. Fig. Fig. 7.3 7.4 7.5 7.6 7.7 7.8 55 Noise Response Fig. Fig. Fig. Fig. Fig. Fig. 7.17 7.19 7.21 7.23 7.25 7.27 Expected Re-projection Misfit from Data Fig. Fig. Fig. Fig. Fig. Fig. Table 7.2: Two-Parameter Results 7.9 7.10 7.11 7.12 7.13 7.14 Noise Response Fig. Fig. Fig. Fig. Fig. Fig. 7.18 7.20 7.22 7.24 7.26 7.28 Chapter 7, Results and Discussion 7.1 FIT T O ORIGINAL OBJECT 56 A N D FIT T O DATA 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 FRLS 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 f than the intersection of the two curves c 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 F B P 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 !| t5 Reprojection misfit to data vs. f \ 0> CO c I m RMS difference from original object vs. f 0.5 (C) 0.8 cycles / pixel 0.4 \* g x 80.35 \ 3= S 0.3 CO 0.25 • , \ 0.2 B / 0.4 , 0.6 . , 0.8 1 in cycles / pixel r in cycles f pixel Figure 7.1: Selected F B P Reconstructions (top), RMS Difference Between FBP Reconstructions and Analytic Circle Phantom vs f (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and Analytic Circle Sinogram vs f (bottom, right) c c Chapter 7. Results and Discussion (A) 0.2 cycles / pixel R M S difference 58 (B) 0.4 cycles / pixel from original object vs. f (C) 0.8 cycles / pixel R e p r o j e c t i o n misfit to d a t a v s . f cr- 0.7 0.4 0.6 f in c y c l e s / p i x e l 0.8 Figure 7.2: Selected FBP Reconstructions (top), RMS Difference Between FBP Reconstructions and MCAT Digital Phantom vs f (bottom, left) and RMS Difference Between Re-projection of FBP Reconstructions and MCAT Monte-Carlo Sinogram vs f (bottom, right) c c Chapter 7. Results and Discussion (A) 0.5 cycles / pixel 59 (B) 0.6 cycles / pixel ^ Least R M S Difference From Original Object o (C) 0.9 cycles / pixel Corresponding Iterations vs. f. 0 o £ 0.23 12 CO ° 0.225 o ^ © 0.22 o 8 c 2> 0.215 0.4 0.6 0.8 f in Nyquist units (cycles / pixel) 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 f (left) and Corresponding Iteration No. vs f (right) - Sample reconstructions use iteration No. corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel GO (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel • Corresponding a vs. f Least R M S Difference from original object 0.35 \ a) o c OJ i_ OJ 0.3 s a co O) o rxO.25 0.5 1 1.5 f in Nyquist units (cycles / pixel) 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.4: Least RMS Difference Between FR-LS Reconstructions and Digital Phantom vs f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 61 (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel mL n • — Corresponding a vs. f Least R M S Difference from original object 0.351 o I 8 0.3 ~° CD o io.25 1 CC 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 Phantom vs f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.5 cycles / pixel 62 (B) 0.6 cycles / pixel (C) 0.9 cycles / pixel Least RMS Difference From Original Object 0.2 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) C Corresponding Iterations vs. f. C 0.2 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) Figure 7.6: Least RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs f (left) and Corresponding Iteration No. vs f (right) - Sample reconstructions use iteration No. corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 63 (B) 0.6 cycles / pixel 0.5 1 1.5 f, in Nyquist units (cycles / pixel) (C) 1.6 cycles / pixel 0.5 1 1.5 f in Nyquist units (cycles / pixel) Figure 7.7: Least RMS Difference Between FR-LS Reconstructions and Digital Phantom vs f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c RA Chapter 7. Results and Discussion (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) f in Nyquist units (cycles / pixel) c Figure 7.8: Least RMS Difference Between FR-ML Reconstructions and Digital Phantom vs f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.5 cycles / pixel 65 (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 15 ca "O a I c o o.i3 i 0.128 1 0 •i—! to '•4—I CD o CD 0-0.126 OCD CC 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 Reconstructions and Sinogram Data vs f (left) and Corresponding Iteration No. vs f (right) Sample reconstructions use Iteration No. corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 66 (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel tip ^ t f f t k ^jtfhk 9 TO Reprojection misfit to data closest to expected value (0.128) Corresponding a vs. f "O £ 0 . 15 w E c 2 0.14 o cn3 o cu o C D .13 ' CC CO cc 0.5 1 1.5 f in Nyquist units (cycles / pixel) 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 f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 67 (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Reprojection misfit to data closest to expected value (0.128) Z M— Corresponding a vs. f 0.16 in I 0.155 c 1 o "I O 0.15 1.0.145 0 -1 CD CC 0.5 1 1.5 co 0.14 f in Nyquist units (cycles / pixel) 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 f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.5 cycles / pixel (B) 0.6 cycles / pixel 0.2 0.4 0.6 0.8 1 f in Nyquist units (cycles / pixel) e 68 (C) 0.9 cycles / pixel 0.2 0.4 0.6 0.8 1 f, in Nyquist units (cycles / pixel) Figure 7.12: Least RMS Difference Between Re-projection of 4 Subset OSEM Reconstructions and Sinogram Data vs f (left) and Corresponding Iteration No. vs f (right) Sample reconstructions use Iteration No. corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 69 (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Reprojection misfit to data closest to expected value (0.107) Corresponding a vs. f n« o 3=0.12 s 10.11 dT °4 o a> 'o Q. CU rr 5 0.1 co rr 0.09 0.5 1 1.5 f in Nyquist units (cycles / pixel) 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 f (left) and Corresponding a vs / (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion (A) 0.4 cycles / pixel 70 (B) 0.6 cycles / pixel (C) 1.6 cycles / pixel Reprojection misfit to data closest to expected value (0.107) Corresponding a vs. f 50.115 | 5 4 0.11 c g 1 0.105 O) o Q. <D rr co rr 0.1 ~2 0.5 1 1.5 f in Nyquist units (cycles / pixel) 1 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 f (left) and Corresponding a vs f (right) - Sample reconstructions use a corresponding to least difference from original object at selected f c c c Chapter 7. Results and Discussion 7.2 71 R E S P O N S E T O NOISE The plots of the norm of the pixel variances versus f for FBP at the top of Figures 7.15 c and 7.16 show a continuous relationship between the norm of pixel variances and the cutoff 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 reconstructions. 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. / " plots of c 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. / " plots more interpretable would c 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 Fig. 7.15 - FBP, Circles Fig. 7.16 - FBP, MCAT Fig. 7.17 - OSEM, Circles, Least Difference from Original Object Fig. 7.18 - OSEM, Circles, Expected Misfit to Data Fig. 7.19 - FR-LS, Circles, Least Difference from Original Object Fig. 7.20 - FR-LS, Circles, Expected Misfit to Data Fig. 7.21 - FR-ML, Circles, Least Difference from Original Object Fig. 7.22 - FR-ML, Circles, Expected Misfit to Data Fig. 7.23 - OSEM, MCAT, Least Difference from Original Object Fig. 7.24 - OSEM, MCAT, Expected Misfit to Data Fig. 7.25 - FR-LS, MCAT, Least Difference from Original Object Fig. 7.26 - FR-LS, MCAT, Expected Misfit to Data Fig. 7.27 - FR-ML, MCAT, Least Difference from Original Object Fig. 7.28 - FR-ML, MCAT, Expected Misfit to Data Max - Mean (Max - Mean) / Mean 1.3 x 10 1.7 x 10 2.8 x 10 3.8 x 10 2.2 2.2 1 x 10 8.4 x 10 8.1 6.1 x 10 5.2 x 10 8.5 3.4 x 10 1.3 x 10 0.39 1.9 x 10 6.9 x 10 0.35 4 x 10 1.8 x 10 0.45 3.7 x 10 1.1 x 10 0.30 3.1 x 10 8.2 x 10 2.7 1.3 x 10 1.1 x 10 8.1 2.9 x 10 8.8 x 10 0.30 3.3 x 10 7.8 x 10 2.4 3.2 x 10 1.9 x 10 3 0.59 1.9 x 10 1.8 x 10 1 0.45 4 4 5 4 4 5 4 3 3 3 5 3 2 3 3 4 4 3 2 3 2 3 4 5 2 2 Table 7.3: Characteristics of "Norm(V) vs. / " Plots c 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 corresponding 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 correlation 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 NormfV) vs. f x 10 Semi-Log NormfV) vs. f C / > § 2 o / C 1.5 B 1 0.5 t 0.2 / / / 0.4. 0.6 I in cycles / pixel 0.8 0.2 Variance plot A (0.2 cycles / pixel) 20 0.4. 0.6 f in cycles /Ipixel f Variance plot B (0.4 cycles / pixel) 40 60 80 Mean Pixel Value Variance plot C (0.8 cycles / pixel) 40 60 Mean Pixel Value Variance image A (0.2 cycles / pixel) 0.8 40 Variance image B (0.4 cycles / pixel) 60 80 100 Mean Pixel Value 120 140 Variance image C (0.8 cycles / pixel) 50 45 40 35 30 25 20 15 10 Co-Variance image A (0.2 cycles / pixel) Co-Variance image B (0.4 cycles / pixel) Co-Variance image C (0.8 cycles / pixel) 200 150 100 50 |o -50 20 40 60 80 100 120 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.15: Noise Response of FBP Reconstructions vs. f for Analytic Circle Data c Chapter 7. Results and Discussion Variance image A (0.2 cycles / pixel) 20 40 60 80 100 120 7G Variance image B (0.4 cycles / pixel) Variance image C (0.8 cycles / pixel) 20 40 60 80 100 120 20 40 60 80 100 120 Figure 7.16: Noise Response of FBP Reconstructions vs. f for MCAT Data c Chapter 7. Results and Discussion 77 Least RMS difference from original object vs f Iteration No. vs f Norm(V) vs f 0.45 '0, 0.4 8 0.35 £ I 0.3 tn 5 0.25 6 4 2 0.2L 0 0.1 0.2 0.4 0.6 f (cycles / pixel) 0.2 20 40 60 80 100 Mean Pixel Value 1 0.2 0.4 0.6 0.S f (cycles / pixel) c Variance plot B (0.6 cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 0.4 0.6 0.8 f (cycles / pixel) 120 50 100 Mean Pixel Value Variance plot C (0.9 cycles / pixel) 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) 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 VA- : 20 40 60 80 100 120 20 40 60 80 100120 Figure 7.17: Noise Response of 4-Subset OSEM Reconstructions vs. f for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion Reprojection misfit lo data closest to expected value vs f 78 Iteration No. vs t Norm(V) vs f 0.2 \ 0) ffl £s 0.18 .1 8 V CO 0.16 a> 7 0.14 0.12. 6 0.2 0.4 0.6 0.8 I (cycles / pixel) 40 60 80 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) C i... 0.4 0.6 0.8 f (cycles / pixel) Variance plot A (0.4 cycles / pixel) 20 B \^ Variance plot B (0.6 cycles / pixel) 120 1 0.2 Variance plot C (0.9 cycles / pixel) 50 100 Mean Pixel Value Variance image B (0.6 cycles / pixel) 0.4 0.6 0.8 t (cycles / pixel) 50 100 Mean Pixel Value 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. Circle Data, with Parameters Resulting in Expected Misfit to Data f for Analytic c Chapter 7. Results and Discussion Least RMS difference from original object vs f E g 79 l ° 9 ( a ) vs f„ Norm(V) vs f 10 0.35 5000 4000 0) 01 as §2 c .£ CO o 0.3 |-3000 c2000 c 0.25 rr 0.5 1000 1 1.5 (cycles / pixel) 0.5 1 1.5 f (cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 Mean Pixel Value 0.5 Variance plot B (0.6 cycles / pixel) 120 Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image B (0.6 cycles / pixel) Variance image A (0.4 cycles / pixel) 40 1 1.5 (cycles / pixel) Variance image C (1.6 cycles / pixel) 80 70 30 60 20 40 10 20 50 30 10 0 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 40 60 120 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100120 Figure 7.19: Noise Response of FR-LS Reconstructions vs. f for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion Variance image A (0.4 cycles / pixel) 80 Variance image B (0.6 cycles / pixel) Variance image C (1.6 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 40 30 20 10 0 Co-Variance image A (0.4 cycles / pixel) 20 |10 40 20 40 60^ 5 60 0 100 80 80 100 120 120 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100 120 10 0 10 0 : 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100120 Figure 7.20: Noise Response of RC-LS Reconstructions vs. f for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data c Chapter 7. Results and Discussion 81 Least RMS difference from original object vs f E o log (a) vs f 1Q Norm(V) vs f c c 0.35 c 6000 B obj "i 0) 0 c 2? to 3= c 0.3 |4000 OJ ra cn o 0.25 2000 rr 0.5 1 1.5 (cycles / pixel) 0.5 r 20 0.5 Variance plot B (0.6 cycles / pixel) Variance plot A (0.4 cycles / pixel) 0 1 1.5 (cycles / pixel) 40 60 80 100 120 Mean Pixel Value Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value 50 100 Mean Pixel Value Variance image B (0.6 cycles / pixel) Variance image A (0.4 cycles / pixel) 1 1.5 (cycles / pixel) Variance image C (1.6 cycles / pixel) 140 120 100 BO 60 40 20 Co-Variance image B (0.6 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) 15 20 10 40' 5 60 0 80 -5 100 -10 120 20 40 60 80 100 120 20 40 60 80 100 120 -15 60 20 20 10 40 60 60 5 80 20 80 0 100 120 1 40 40 0 100 -5 20 40 60 80 100 120 120 20 40 60 80 100120 J-20 mL 20 40 60 80 100120 20 40 60 80 100 120 Figure 7.21: Noise Response of FR-ML Reconstructions vs. f for Analytic Circle Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion 82 Reprojection misfig to data closest to expected value vs f c o s3 <D CO •6-° cLS <u — ft log (<x) vs f 10 Norm(V) vs f c 5000 0.16 0.155 5-4000 41 0.15 £3000 0.145 0.14 2000 0.5 1 1.5 f (cycles / pixel) 0.5 Variance plot A (0.4 cycles / pixel) 0 20 40 60 80 100 Mean Pixel Value 1 1.5 (cycles / pixel) 0.5 1 1.5 f (cycles / pixel) Variance plot B (0.6 cycles / pixel) 120 Variance plot C (1.6 cycles / pixel) 50 100 Mean Pixel Value Variance image A (0.4 cycles / pixel) 50 100 Mean Pixel Value Variance image B (0.6 cycles / pixel) 1 50 Variance image C (1.6 cycles / pixel) 70 100 60 40 80 30 60 20 40 50 40 30 20 10 20 0 0 10 Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) 1 1 20 40 60 o 80 100 120 20 40 60 80 100120 I 40 20 30 40 20 60 10 80 0 100 .-5 120 40 20 0 20 40 60 80 100120 Co-Variance image C (1.6 cycles / pixel) 0 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100120 20 40 60 80 100 120 Figure 7.22: Noise Response of RC-ML Reconstructions vs. f for Analytic Circle Data, with Parameters Resulting in Expected Misfit to Data c -10 Chapter 7. Results and Discussion 83 Least RMS difference from original object vs f Iteration No. vs f 0.82 Is e si il II O) 0.8r 0.78 V 0.76 0.74 C 0.72 0.4 0.6 0.£ f (cycles / pixel) 0.2 Variance plot A (0.4 cycles / pixel) 0.4 0.6 0.8 f. (cycles / pixel) 1 0.2 Variance plot B (0.6 cycles / pixel) 0.4 0.6 0.8 f (cycles / pixel) Variance plot C (0.9 cycles / pixel) 450 2500 400 350 2000 2.6 • s ta 2 11500 lis > 300 I 5200 1000 + 250 150 500 • 10050, 50 100 Mean Pixel Value 150 50 Variance image A (0.4 cycles / pixel) 100 150 Mean Pixel Value 50 Variance image B (0.6 cycles / pixel) 1 00 150 Mean Pixel Value Variance image C (0.9 cycles / pixel) 350 300 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) 60 40 100 20 40 150 60 20 0 20 40 60 80 100 120 V_/\ : 20 40 60 80 100 120 60 30 20 40 60 80 100 120 20 40 60 80 100 120 -20 80 100 120 100 50 0 20 40 60 80 100 120 20 40 60 80 100 120 Figure 7.23: Noise Response of 4-Subset OSEM Reconstructions vs. f for MCAT Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion 84 Reprojection misfit to data closest to expected value vs f Iteration No. vs f Norm(V) vs f 0.15 12 0.14 10 £ 0.13 ssl 8 ¥ 6 o 0.12 8 0.11 • 4 2 0.1 0 ).4 0.6 0.8 (cycles / pixel) 1 0.2 Variance plot A (0.4 cycles / pixel) 0.4 0.6 0.8 f (cycles / pixel) 1 0.2 Variance plot B (0.6 cycles / pixel) 50 100 Mean Pixel Value 50 Variance image A (0.4 cycles / pixel) 0.4 0.6 0.8 f (cycles / pixel) Variance plot C (0.9 cycles / pixel) 100 Mean Pixel Value 50 Variance image B (0.6 cycles / pixel) 100 Mean Pixel Value Variance image C (0.9 cycles / pixel) 100 80 no 40 20 Co-Variance image A (0.4 cycles / pixel) Co-Variance image B (0.6 cycles / pixel) I' 1 Co-Variance image C (0.9 cycles / pixel) 6 20 20 40 40 60„ 60 80 80 100 100 120 -2 20 40 60 80 100 120 s r 4 0c 20 40 60 80 100 120 1 I, 0 120 18 20 40 60 80 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 Figure 7.24: Noise Response of 4-Subset OSEM Reconstructions vs. f for MCAT Data, with Parameters Resulting in Expected Misfit to Data c Chapter 7. Results and Discussion 85 Least RMS difference from original object vs f log (a) vs f 10 Noim(V) vs f c E 4000 o * T3 gf 0.76 3= . g ^ .2> W o 5 0.74 rr 3 cn o ~3000 4 3 §2000 c C 0.72 0.5 1000 1 1.5 (cycles / pixel) 0.5 100 Mean Pixel Value 0.5 Variance plot B (0.6 cycles / pixel) Variance plot A (0.4 cycles / pixel) 50 1 1.5 (cycles / pixel) 150 50 Variance plot C (1.6 cycles / pixel) 100 150 Mean Pixel Value 50 Variance image B (0.6 cycles / pixel) Variance image A (0.4 cycles / pixel) 1 1.5 (cycles / pixel) 100 150 Mean Pixel Value Variance image C (1.6 cycles / pixel) co so 51) 40 40 30 .in 20 20 10 10 0 Co-Variance image B (0.6 cycles / pixel) Co-Variance image A (0.4 cycles / pixel) 20 30 40 20 60 10 80 100 0 120 -10 Co-Variance image C (1.6 cycles / pixel) 40 20 20 140 40 60 20 80 100 0 120 20 40 60 80 100 120 20 80 10 0 120 20 40 60 80 100120 n 60 100 40 20 0 40 20 20 40 60 80 100 120 30 40 20 40 60 80 100120 20 40 60 80 100120 --10 20 40 60 80 100 120 Figure 7.25: Noise Response of FR-LS Reconstructions vs. f for MCAT Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion Variance plot A (0.4 cycles / pixel) 86 Variance plot B (0.6 cycles / pixel) Variance plot C (1.6 cycles / pixel) Figure 7.26: Noise Response of RC-LS Reconstructions vs. f for MCAT Data, with Parameters Resulting in Expected Misfit to Data c Chapter 7. Results and Discussion Co-Variance image A (0.4 cycles / pixel) 87 Co-Variance image B (0.6 cycles / pixel) Co-Variance image C (1.6 cycles / pixel) • 60 20 20 |40 40 40 60. 80 100 80 100 120 40 20 0 60 20 120 J— 20 40 60 80 100120 20 40 60 80 100120 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. f for MCAT Data, with Parameters Resulting in Least Difference from Original Object c Chapter 7. Results and Discussion 88 Reprojection misfit to data closest to expected value vs f o •s a CD 10 0.115 5 0.11 4 CO '5"° QCD log (a)vsf S„ 300 5-200 3 0.105 S 0.1 0.5 1 1.5 f (cycles / pixel) r Variance plot A (0.4 cycles / pixel) 50 100 Mean Pixel Value Norm(V) vs f c I100 2 1 0 1 0.5 1.5 f, (cycles / pixel) 0.5 Variance plot B (0.6 cycles / pixel) 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) 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 2 40 1.5 60 1 80 0.5 100 0 120 -0.5 20 40 60 80 100120 20 40 60 80 100 120 20 40 60 80 100120 ±Z2 20 40 60 80 100120 20 40 60 80 100120 Figure 7.28: Noise Response of RC-ML Reconstructions vs. / for MCAT Data, with Parameters Resulting in Expected Misfit to Data c Chapter 7. Results and Discussion 7.3 89 SUMMARY Results Fig 7.2 Panels A & B and other Figs Figs. 7.1 and 7.2 Figs. 7.4 and 7.7 Figs. 7.3 and 7.5 Figs. 7.3 and 7.5 7.12, 7.14 Figs. 7.15 and 7.16 Fig. 7.23 Figs. 7.17, 7.18, 7.23 and 7.24 Features 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 Streaking artifacts typical of FBP reconstructions, especially in low count areas of the image Slight streaking around low intensity edges due to the FR-LS least squares fit term Different noise textures between OS-EM and FRCT reconstructions, especially in the low intensity extended source FR-ML achieves better contrast than OSEM in cold spots Continuous relationship between the norm of pixel variances and FBP cut-off frequency Highly unstable nature of OSEM at high numbers of iterations. Reconstructing more random noise realizations may yield more enterpretable results 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. / " plots. Reconstructing more random noise realizations and using finer resolution in parameter sweeps may yield more enterpretable results c Table 7.4: Summary of Results - Part 1 Chapter 7. Results and Discussion Results Table 7.3 Table 7.3 Figs. 7.17, 7.18, 7.23, and 7.24 Figs. 7.19-7.22, and 7.25-7.28 Figs. 7.1 and 7.2 Figs. 7.26 and 7.28 90 Features Parameter settings which attempt to generate expected misfit to data result in reconstructions which are more stable with respect to noise. FRCT reconstructions have consistently lower pixel variances than corresponding OSEM and FBP reconstructions Linear relationship between the mean pixel value of OSEM reconstructions and corresponding variance in scatter plots Stabilizing effect of Fourier Regularozation can be seen in the shape of the scatter plots Co-variance images show the relationship between correlation lengths in the reconstructed image and FBP cut-off frequency 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 CHAPTER 8 CONCLUSION 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 GENERAL COMMENTS 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 8.1.1 FUTURE 92 WORK F L U X CONSERVATION TERM In the field of Baysian image processing, adding a "flux conservation" term to the maximum entropy objective function has been proposed [31]. This method could be applied to FRCT reconstruction, in order to associate physical units such as ™™t with the c a reconstructed measurements. T H E SYSTEM MATRIX 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 reconstruction 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. REGULARIZATION 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 regularization 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 GEOMETRIES AND 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 appropriate 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. BIBLIOGRAPHY G. Hounsfield. Computerized transverse axial scanning (tomography): Part I. Description 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 Company, 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 annihilation 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 emission computed tomography in the presence of a variable attenuating medium. Technical 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. CRC Press, Boca Raton, 1996. The Radon Transform and Local Tomography. 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 tomography. 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. Evaluation 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 MLEM reconstructed emission tomographic images. IEEE Trans. Nuc. Sci., 40(4):11981203, 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. CHAPTER 9 APPENDIX A 97 Chapter 9. Appendix A High Cutoff Frequency (0.9 cycles / pixel) High Iteration No. (8) gg High Cutoff Frequency (0.9 cycles / pixel) 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 f and Iteration No. - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A Low Cutoff Frequency (0.6 cycles / pixel) Low a (10 1) A 99 Low Cutoff Frequency (0.6 cycles / pixel) High a (10 4) A iog (a) 10 High Cutoff Frequency (1.2 cycles / pixel) High a (10*5.5) f in cycles / pixel Figure 9.2: Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital Circle Phantom vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A Low Cutoff Frequency (0.6 cycles / pixel) Low a (10 -1) A 1 IMA\n 100 Low Cutoff Frequency (0.6 cycles / pixel) High a (10 1.5) A log (a) 10 35 4 is High Cutoff Frequency (1.6 cycles / pixel) High a (10 3) A f in cycles / pixel Figure 9.3: Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital Circle Phantom vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A Low Cutoff Frequency (0.5 cycles / pixel) High Iteration No. (9) 101 Low Cutoff Frequency (0.5 cycles / pixel) Low Iteration No. (5) High Cutoff Frequency (0 9 cycles / pixel) Low Iteration No. (4) Figure 9.4: Surface Plot of RMS Difference Between 4 Subset OSEM Reconstructions and Digital MCAT Phantom vs f and Iteration No. - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A Low Cutoff Frequency (0.6 cycles / pixel) Low a (10 1) A 102 Low Cutoff Frequency (0.6 cycles / pixel) High a (10*5) High Cutoff Frequency (1 2 cycles / High a (10*6.5) Figure 9.5: Surface Plot of RMS Difference Between FR-LS Reconstructions and Digital MCAT Phantom vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A Low Cutoff Frequency (0.6 cycles / pixel) Low a (10 1) A 103 Low Cutoff Frequency (0.6 cycles / pixel) High a (10 2.5) A High Cutoff Frequency (1 2 cycles / High a (10 4) A Figure 9.6: Surface Plot of RMS Difference Between FR-ML Reconstructions and Digital MCAT Phantom vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A High Cutoff Frequency (0.9 cycles / pixel) High Iteration No. (7) 104 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 f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A 105 Low Cutoff Frequency (0.6 cycles / pixel) Low (10 I rt..Fa 11AA11.5) C \ Low Cutoff Frequency (0.6 cycles / pixel) I T : _I_ a / 1 r\ 4.5) r\ High (10 A A A < High Cutoff Frequency (1.6 cycles / Dixel) .... _ 5.5) . _ £. • ' High a (10 A High Cutoff Frequency (1.6 cycles / pixel) Low a (10 2.5) A iii n l iog (a) 10 e f in cycles / pixel Figure 9.8: Surface Plot of RMS Difference Between Reprojection of FR-LS Reconstructions and Sinogram Data vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A 106 Figure 9.9: Surface Plot of RMS Difference Between Reprojection of FR-ML Reconstructions and Sinogram Data vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A 107 Figure 9.10: Surface Plot of RMS Difference Between Reprojection of 4 Subset OSEM Reconstructions and Sinogram Data vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A 108 Figure 9.11: Surface Plot of RMS Difference Between Reprojection of FR-LS Reconstructions and Sinogram Data vs f and a - Sample reconstructions show effects of faulty parameter settings c Chapter 9. Appendix A 109 Figure 9.12: Surface Plot of RMS Difference Between Reprojection of FR-ML Reconstructions and Sinogram Data vs f and a - Sample reconstructions show effects of faulty parameter settings c
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Fourier regularization applied to spect reconstruction
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Fourier regularization applied to spect reconstruction Togane, Dylan Dalmar 2000
pdf
Page Metadata
Item Metadata
Title | Fourier regularization applied to spect reconstruction |
Creator |
Togane, Dylan Dalmar |
Date Issued | 2000 |
Description | Fourier Regularized Emission Computed Tomography (FRECT), a reconstruction method for Emission Computed Tomography (ECT), is presented and compared with other reconstruction 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 experiments 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. |
Extent | 21698290 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085173 |
URI | http://hdl.handle.net/2429/10986 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
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
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085173/manifest