UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Scatter correction in Positron Emission Tomography Sirota, Dimitri 1995

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

Item Metadata


831-ubc_1995-0644.pdf [ 7.55MB ]
JSON: 831-1.0085133.json
JSON-LD: 831-1.0085133-ld.json
RDF/XML (Pretty): 831-1.0085133-rdf.xml
RDF/JSON: 831-1.0085133-rdf.json
Turtle: 831-1.0085133-turtle.txt
N-Triples: 831-1.0085133-rdf-ntriples.txt
Original Record: 831-1.0085133-source.json
Full Text

Full Text

S C A T T E R C O R E C T I O N I N P O S I T R O N E M I S S I O N T O M O G R A P H Y by D I M I T R I S I R O T A B.Sc. McGil l University 1992 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M . S C . i n THE FACULTY OF GRADUATE STUDIES in T H E D E P A R T M E N T O F P H Y S I C S We accept this thesis as conforming to the required standard T H E UNIVERSITY OF^RITISH C O L U M B I A October 1995 ©Dimitr i Sirota 1995 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 The University of British Columbia Vancouver, Canada Date ficl /3 hS DE-6 (2/88) Abstract Abstract: Positron Emission Tomography (PET) aims to identify metabolic functioning through imaging the in situ physiological progression of a metabolic analogue tagged with a positron emitting radio-isotope. Following ejection, the positron promptly annihilates with a surrounding electron giving rise to two oppositely directed y-rays of equal energy. By coincidentally detecting the collinear y-rays, PET is able to characterize the underlying metabolite distribution and concentration. Physical scattering of one of the two coincident y-rays inside the medium supporting the metabolite will result in a mispositioning of the decayed radio-tracer. The set of all scattered events, following acquisition, will contribute background to the set of unscattered or primary events. Camera sensitivity to this activity concentration is improved by acquiring events over all possible coincident detector pairs. This kind of volumetric or 3D PET imaging also increases the likelihood that the camera will accept scatter events. The present research investigates a spatially variant convolution type method to extricate the scatter component from the total measured events in 3D imaging. In the scope of the. present work a spatially variant kernel will be derived which upon convolution with measured 3D PET image data will return a matching scatter component. Based on this, kernel a computer algorithm is developed, implemented and performance tested on known activity distribution sources. It is found that the present algorithm performs better or of equal merit in comparison to,, pre-existing scatter correction algorithms. Moreover, the position variable algorithm is optimized to perform within a clinical setting. ii Table of Contents T a b l e o f C o n t e n t s : A B S T R A C T T A B L E O F C O N T E N T S L I S T O F T A B L E S L I S T O F F I G U R E S 1, INTRODUCTION _ l 2. P O S I T R O N E M I S S I O N T O M O G R A P H I C I M A G I N G 5 2.1. P H Y S I C S 5 2.1.1. TRACER KINETICS 5 2.1.2. EVENT DETECTION • . 6 2.1.3. VOLUMETRIC IMAGING 7 2.1.4. E C A T CAMERA CONFIGURATION 8 detector 8 2D Vs 3D Plane Definitions 9 2.i:4.3. 2D vs. 3D Sensitivity Differences 10 Signal-to-Noise Performance Gain in 3D Over 2D .12 2.2. R E C O N S T R U C T I O N 13 2.2.1. EVENT ACQUISITION AND REGISTRATION • . 1 3 2.2.2. DATA ORDERING - T H E SINOGRAM 14 2.2.3. Sinogram Inversion - Image Reconstruction. * 16 3. S C A T T E R ' •• ; - ' ' 18 3.1. A Q C U I S I T I O N M I S R E G I S T R A T I O N 18 3.2. P H Y S I C A L O R I G I N O F S C A T T E R 19 iii Table of Contents 3.3. MANIFESTATIONS OF SCATTER BASED ERROR 23 3.3.1. RANDOM EVENTS Y 24 . 3.3.2. ATTENUATION ./• 25 3.3.3. MISREGISTRATION !- 27 3.4. SCATTER IMPACT 29 4. S C A T T E R C O R R E C T I O N • • 30; 4.1. SCATTER SUPPRESSION 30 4.2. SCATTER SUBTRACTION 31 4.2.1. SEGMENTED ENERGY WINDOW METHODS 31 4.2.2. INTEGRAL TRANSFORMATION 34 deconvolution 36 Iterative Convolution Subtraction ' 37 convolution subtraction " ' 3 7 5. M E T H O D F O R D E T E R M I N I N G T H E K E R N E L H 40 5.1. THEORY 40 5.2. DERIVING RELATIONSHIP BETWEEN S A N D H. 41 5.3. IMPLEMENTATION 41 5.3.1. MATERIALS . 42 5.3.2. SAMPLING PROCEDURE , ;. - ; 44. 5.4. PLANAGRAM CONSTRUCTION 45 6. I M P L E M E N T A T I O N 47 6.1. SELECTING SM O D E L " . . 47 6.1.1. POSITION VARIABLE SCATTER RESPONSE IN T H E PROJECTION . • '. 47 Motion Along L O R - • %. 48 Displacement Perpendicular to the L O R (Transverse Direction) 52. Axial Variation Perpendicular To L O R 54 6.2. DETERMINING THE M O D E L SCATTER FUNCTION 56 6:2.1. CONSTRAINTS ON CANDIDATE FUNCTION 56 - . 6.2.2. CANDIDATES • , . . 57 6.3. ALGORITHM 59 6.3.1. OPTIMIZATION CRITERION . . 5 9 6.3:2. SELECTION OF MASKS : , 62 7. E X E C U T I O N 67 7.1. SELECTING THE M O D E L 67 7.2. CALCULATING H 72 7.2.1. A, RESULTS Y - - 73 7.2.2. A 2 RESULTS 75 . 7.2.3. A 3 RESULTS \ 76 Table of Contents 7.2.4. A , R E S U L T S 7.2.5. A 5 R E S U L T S 7.3. G E N E R A T I N G T H E C O M P L E T E K E R N E L 7.3.1. F O R M I N G T H E P A R A M E T E R DISTRIBUTIONS 7.3.2. A , 7.3.3. A 2 7.3.4. A 3 -7.3.5. A 4 7.3.6. A 5 7.4. I N T E R P R E T A T I O N 8. T E S T O F C O R R E C T I O N 8.1. A L G O R I T H M 8.1.1. N O R M A L I Z A T I O N 8.2. A L G O R I T H M D E M O N S T R A T I O N 8.3. C O M P A R I S O N O F T H E C O R R E C T I O N E F F I C A C Y 9. C O N C L U S I O N V List of Tables . L i s t of T a b l e s : Table 2-1: Properties of PET Radionucleides Table 6-2 Selection of Inner Mask Table 7-1 yf Comparison For Fit on Central Source Table 7-2 yf Comparison For Fit on Offset Source List of Figures L i s t of F i g u r e s : Figure 2-1: 2D Projection of 3D Source ' - - ' • ' .... • ' • ' 7 Figure.2-2 2D ECAT Plane Definitions - ' ^  - - / • - 10; Figure 2-3 Axial Sensitivity: 2D Vs 3D Imaging •' - • ^_ 11 Figure 2-4: Demonstration of a Tomographic Projection and Sinogram Construction 15. Figure 3-1 Misregistration Due to Scatter. • . - 18 Figure 3-2 Relative Contribution of Photoelectric Affect, Compton Scatter, and Pair. Production as a Function of • • Energy1 • • •''. - ' ' . . 20 Figure 3-3 Energy Distribution as a Function of Scattering Angle for 511 KeV Photons • • .. 21 • Figure 3-4 Mashed Projection of a Centralized Point Source (259,760 counts) 28 Figure 5-1 Source, Holder and Phantom . ' . 44 Figure 5-2 Sinogram Rebinning Into Planagrams ' 45 Figure 6-1' Chord Specification In Object Space ; • • 48 Figure 6-2 Simulation44 • • - • • • ; 49 Figure 6-3 Transverse Comparison ' . ' 50 Figure 6-4 Imaging Configuration ' " • ' ', : • 51 Figure 6-5 Configuration for Transverse Comparison • 52 Figure 6-6 Transverse Cross-Section , - , . . 5 3 Figure 6-7 Configuration for Axial Comparison ' ._ ' 55 Figure 6-8 Axial Cross-Section •- 55 Figure 6-9 Magnified View of Scatter From Centralized Point Source ' 63 vi i List of Figures Figure 6-10 cross-section of Radially Offset Source, With and Without End Mask • ___66 Figure 7-1 Scatter Planagram Distributions From Central Source . • •' 68 • - Figure 7-2 Transverse Cross-Section "•' • ' ". ' • _^ ___ . 69 Figure 7-3 Transverse Cross-Section ;. 71 Figure 7-4 Scatter Planagram Distributions From Offset Source ': . - —72 Figure.7 5 a, Distribution . - ' • - • " : : __]_—_74, Figure 7-6 Comparison of Measured Counts as a Function of Radial and Axial Projection Position_ '.'75 •Figure 7-7 a. Distribution _J - __ : — : 76• Figure 7-8 a3 Distribution 1 - - ." • - . '_ , '_ : : 77 : .Figure 7-9a, Distribution • • ; • . ; . . —; .78 Figure 7-10 a5 Distribution " ._ : ' 79 Figure 7-11 Look Up Table Implementation Of Assigning h -»: ; . 80 , Figure 7-12 Planagram Distribution For at ' '• 1 •. '. 83 Figure 7-13 Planagram Distribution For a2 - ••' . • - • '• • 84 Figure 7-14 Planagram Distribution For a, ' • • • : . • -85 Figure 7-15Planagram Distribution For'a4 '" • ' • ' ' •"' ' • • ' 86 Figure 7-16 Planagram Distribution for as ! . - ' • ' _87 Figure 7-17 Scatter Fraction Distribution on the Planagram Domain ".• • • ' ' ' 8 8 Figure 8-1 Convolution Algorithm •' ' - • - • 91 • Figure 8-2 SF Measurements From The Offset Source ' . • • •'- • ' 94 Figure 8-3 SF Sampling Procedure • ':' • ' • /• 96 Figure 8-4 Test Phantoms " ' • • • . , '. • • 102 Figure 8-5 Planagram View of Scanned Central! lot Spot Object : '" •. ' 103 • Figure 8-6 Variation in yf and Ax' with Successive Iterations •. : '" 104 Figure 8-7 Evolution in Scatter and Primary Count Planagram Distributions '' 105 Figure"8-8 Transverse View Through Planagram: Original, Corrected, and Scatter 106 ' Figure 8-9 Axial View Through Planagram: Original, Corrected. Scatter •_ 107 • • ' ' '• •. ;.viii' \ List of Figures Figure 8-10 ROI Positioning Inside the Test Phantoms ; : _109 Figure 8-11 Activity Contrast in the Hot Spot Source : : 110 Figure 8-12 Hot Spot to Background Activity Contrast in Hot and Cold Spot Source 111 Figure 8-13 Cold Spot to Background Activity Contrast in Hot and Cold Spot Source • 112 ix Acknowledgements A c k n o w l e d g e m e n t s : Fde like to thank the people of PET: Vesna, thanks for all the help and most of all for hanging on (almost to the end of the millennium no less) and not losing faith in me; Tom and Dick thanks for not thinking me too much a flake (I r . hope); To the little people: Scott, Ken, Theresa and the remaining coterie - I owe you (nothing too expensive mind you). Thanks to everybody. later! dimitri Introduction 1. INTRODUCTION Modern medical imaging has greatly advanced the practice of medicine. Current imaging techniques permit the clinical evaluation of biological structure and function without the need for surgically invasive procedures. Most importantly the ability to render form and action has proven panacea to the greater understanding of human pathology. Equipped with this tool some ischemia, disease and cancer may be.quickly and cost effectively characterized, possibly delaying morbidity and extending life. The efficacy towards both clinical diagnosis and scientific investigation visited by the new visualization schemes is validated by the increasing role assumed by the imaging laboratory inside the hospital setting. Each imaging modality aims to reconstruct an in situ diagram of anatomy or physiological functioning based on either the transmission or emission of radiation (possibly non-ionizing such as acoustical) through the scanning subject. A smaller class of these modalities, notably CT, SPECT and PET reconstruct their internal diagrams based on a series of detected radiation samples or measurements taken at discrete points surrounding the subject. This class of detection and reconstruction comes under the general rubric of tomographic projection imaging. Both PET1 (Positron Emission' Tomography) and SPECT2 (Single Photon Emission Computerized Tomography) measure the tomographic projections from randomly emitted photons (gamma - "y" radiation), emanating either directly from the decay of a radio-tagged pharmaceutical (SPECT) administered, to the host or by way of an ancillary e +-e ': arinihilation process (PET). Since the biodistribution of the 1 Introduction ; . radiopharmaceutical is known to follow either a specific metabolic path or concentrate in an organ with some characteristic metabolic rate the temporal collection of emitted radiation facilitates the investigator to glean insight into physiologic functioning within a precise pathway or organ. In the case of PET, the collected tomographic data can upon reconstruction, be transposed into quantitatively meaningful physiologic parameters. In exact terms, PET measures in vivo metabolic activity based on the detection of coincident 511 KeV fs emitted from the annihilation of a positron (fj +) i.e. a positively charged electron ejected from a positron emitter labeled source internally administered to the subject (through inhalation or injection). Reconstructing the source activity distribution that gave rise to the measured coincident events therefore permits the medical investigator to assess physiological function non-invasively, in-situ. For example, lsO-labeled water may be used to investigate regional brain function (activation studies) based on the differential blood (oxygenation) requirements of separated brain regions. However, PET unlike the related modality of SPECT or other functionally capable imaging modalities, will determine quantitatively meaningful physiological parameters such as blood flow in mL.min^.lOOg1 of tissue or biochemical receptor density and affinity of binding KDof a particular neuroreceptor.3 This ability for quantitation, due to the ability to correct for medium attenuation, permits the derivation of exact biochemical function models when PET images are used in conjunction with blood chemistry and biomathematical paradigms: 82Rb labeled RbCl used in myocardial perfusion imaging indicates heart muscle injury; 2 Introduction the metabolic uptake of glucose in the. heart and brain using the 1 8 F labeled analogue Fluorodeoxyglucose (FDG) will indicate organ metabolism. At the shared U B C / T P v I U M F PET facility, research is focused exclusively on the brain. The local effort is devoted to the investigation of neurodegenerative diseases with particular emphasis on Parkinson's disease. This disease, characterized by a deficit of the neurotransmitter dopamine, is investigated through the targeted application of radioligands such as 1 8F-dopa and nC-Raclopride that concentrate in the striatum, a crucial part of the brain involved in the motor control system. The radioligands address the presynaptic vesicle dopamine storage capabilities and the postsynaptic D 2 dopamine receptor concentration. Characterization of the disease requires the detection of differential radioligand uptake in the striatal substructure, caudate and putamen. This level of differentiation demands high resolution and high sensitivity imaging capabilities. The highest sensitivity is achieved by volumetric imaging capabilities (see section 2). High sensitivity is desired also in scanning situations where the radioactivity dose administered to patients is of special concern (infant scanning; repeat studies on the same subject) or when reduced scanning times are crucial to the operation of the facility, as for example in a purely clinical setting. A description of volume imaging is conveyed in the proceeding section. However, worth noting here is that attended with the perfunctory execution of a scan is the losseous (in relation to the collected signal i.e. coincident counts) effect of scatter, whereby otherwise detectable 7 events are scattered from their emitted trajectory by the biological material surrounding the radiation source (herein the "source"). Introduction Invariably. the scatter degrades the statistical quality of the collected events. The subsequent image derived from the data fails to be quantitative. For this reason removing the scatter is a key component of the quantitation protocol. The scatter correction method expounded within this dissertation is an evolutionary construction on earlier work in the field. The deviation from previous scholarship lies in the attempt at accounting for the spatial variability of the scatter in the two dimensional image projections generated in volumetric (3D) scans. A correction method is developed and implimented whereby a unique scatter response function is determined for each point in the 2D tomographic projection. Upon convolution with the measured data (or estimate thereof) this "position variable" kernel returns the full scatter distribution across the projection. Following a review of PET (chapter 2) and the general form of the scatter problem in it (chapter 3) , the type of scatter correction ultimately adopted in the course pf this research is revealed in the context of correction methodologies in general (chapter 4).. . Chapter 5 develops the theory behind deriving the projection scatter profiles from the position variable kernel. The actual derivation of the kernel is left to chapter 6 and 7. Chapter 8 examines the practical performance capability of the position variable scatter correction algorithm in comparison to alternate scatter correction methods. It is found that the position variable correction performs equal to better than the contrasted methods. 4 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING 2. Positron Emission Tomographic Imaging 2.1. P H Y S I C S 2.1.1. T R A C E R K I N E T I C S PET is based on the visualization of a radiopharmaceutical that is inside the individual being imaged. Below is a list of the common PET radionucleides:2 Table 2-1: Properties of PET Radionucleides "C 1 3N 1 5o ' 1 8 F. 6 8Ga 82Rb HalMife 20.4 min 10 min 2 min 110 min 68 min 75 sec E m V M e V ) .96 1.20 1.74 .633 1.899 3.4 A l l the above elements possess half-lives of under 2 hours thus minimizing ionizing radiation exposure to the patient during circulation (following injection of the radiopharmaceutical the usual bodily clearance,facilities metabolize and excrete the radiolabeled agents).4 Upon introduction into the patient, the radio-compound follows the metabolic pathway of their naturally occurring analogue.4 During the course of its progression along the metabolic pathway radio-label decays by the emission of a positron. The energetic positron rapidly slows by coulomb interaction with surrounding electrons to thermal energies in situ. The distance the positron strays from the generating decay event is proportional to its ejected energy. Therein the positron will annihilate with a 5 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING neighboring electron giving rise to two oppositely directed 511 KeV Ys. Two physical factors limit the resolvability of a PET camera. The path length of the positron between its origin and annihilation places an emitter dependent limit on the spatial resolution. Secondly, the co-linearity of the emitted Ys is not exact. In the laboratory reference frame the e+ - e will conserve the diminished although non-zero momentum from the ejected (3+. In the lab frame the opposingly directed Ys will thus deviate from complete 180° separation by p . 2.1 0 s ' * m . c e —e which is typically on the order of .25°. 2 The combined affects confer a lower resolvability limit for high resolution PET scanners of 2mm-3mm.5 2.1.2. E V E N T D E T E C T I O N The PET camera is comprised of y-detectors arranged in either multiple circular.or hexagonal rings positioned contiguous to each other around a central aperture. Typically the detection material is comprised of scintillating crystal with high Z to favour the photoelectric effect cross section, such as bismuth germanate (BGO), thallium doped sodium iodide Nal(TI), barium fluoride (BaF2) or cesium fluoride (CsF). Of these BGO is most prevalent because of the material's relatively high intrinsic efficiency defined as the fraction of impinging radiation giving rise to a signal pulse, and relatively high energy resolution (21% at 511 KeV). Also, BGO unlike the alkali halide salts like Nal is not hygroscopic - as this will tend to distort the crystal color and so attenuate the emitted optical light energy.2 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING The selection of crystal material and topology determines the camera efficiency and resolution. Small crystal size raises the intrinsic camera resolution at the possible expense of decreased sensitivity (more empty space between crystal elements). The crystals are optically coupled to photomultiplier tubes (PMTs)2 which transform the optical photo-pulse, generated inside the crystal by an incident y, into a voltage pulse proportional in size to the deposited energy of the yi 2.1.3. V O L U M E T R I C I M A G I N G Given a tracer distribution f(x,y,z) in the aperture, a volumetric imaging protocol will generate a set of 2D shadow projections p(xr,y^0,^) (see Figure 2-1) of the activity by acquiring coincident y-events across all opposing detector pairs. Figure 2-1: 2D Projection of 3D Source 6 Object: f ^ V ^ A x Projection: pOc^^e) S = .X. 9 cose sine The collected event or counts at a given 6,cp represents the integral of activity through the source at that orientation. If events, through physical interdiction by y-opaque 7 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING ! barriers called septa or acquisition logic, are selectively restricted to a small sub-set of all the possible planes intersecting the object, the corresponding imaging program is said to be slice-oriented 2D. The current art in PET however makes use of full volume 3D imaging in order to optimize the detection of the available source activity signal.1 The improved sensitivity afforded by 3D over 2D imaging acts to augment the signal-to-noise ratio of the reconstructed source image.7 Particularly, under scanning protocols where source activity levels become a medical concern, the utility of 3D over 2D imaging becomes most evident (see. section j The proceeding research is conducted for 3D imaging. A review of the equipment used in the course of the current research is conducted below. 2.1.4. E C A T CAMERA CONFIGURATION DETECTOR \ ! The current work is performed on the Siemens/CTI ECAT 953B neurocamera. The camera is comprised of a series of y-sensitive BGO detectors arranged into two adjacent rings, 76.5cm in diameter each, divided into 12 buckets of 4 detector crystal blocks 5.0cmx5.4cmx30.0cm(deep) in dimension. ;In order to improve spatial resolution, each block is saw cut (with .6mm crystal gaps of various depths) i into an 8x8 array of individual crystals elements - 5.66mm wide and 6.15mm long. Each crystal block module is optically coupled to four Hammamatsu photomultiplier tubes. The many-to-one ratio of crystal elements to PMTs requires that correct event assignment (to the scintillating crystal element) be based on. the 8 ! Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING differential light output from the crystal at the input of the PMTs (as light leakage between neighboring crystal elements is unavoidable)8 An Anger Logic9 light positioning scheme is invoked on the ECAT in order to select the correct scintillating crystal segment.. 2 D Vs 3 D P L A N E D E F I N I T I O N S I .' • ' i The ECAT is capable of acquiring events in either slice-oriented 2D or volumetric 3D configurations. The 10.6 cm axial field of view (FOV) consists of the 16 crystal rings (8 crystal rings/block x lblock/rihg x 2rings) generating 16 transaxial direct planes (formed from opposing detector pairs). In 2D mode the ECAT supports 15 interdicting 9 cm long annuli of Tungsten septa that extend into the camera FOV between the 16. crystal rings. The y-opaque septa restrict axial activity sampling to rigidly defined transverse planes. By limiting the angular acceptance angle, the number of scattered events (see section 3 for a discussion of y ray interaction with matter) is kept to a minimum at the expense of overall camera sensitivity. However, rather than limit acquisition to the strictly defined direct planes, improved axial sampling and statistics is achieved by combining adjacent planes to form 31 effective axial'planes consisting of alternating direct and cross-planes (see Figure 2-2) axially separated by .3375cm. The y rays detected in a strictly defined plane together with those detected in the crystal rings seperated by. ring difference 2 and axially centered on the direct plane will define new direct acquisition planes. The coincident events where y rays are detected in crystal rings; seperated by ring 9 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING '• . . ! ' . " difference 1 and 3 will form cross acquisition planes that are.axially centered between two direct planes. The 31, camera defined, 2D direct and cross-plane definitions are illustrated below: v I, Figure 2-2 2 D E C A T Plane Definitions ! jiiriiiiiiiji] i i i i iK i i i i i ; < * X K X iiijiiiiiiix ^ i i i i i i i i l l i i l i i l i l i l i K i i i i i \ \ \ • 1 \ : . [ i . .i i / . / / "t I. • 1 7 I 1 . j \ I \.\ 1 f t J J < i. s / // i t 1 / 1 f 1 I 1 ' ...Cross Plane Definition i\ hoi Plane 28 i \ f • \ V Direct Plane Definition For Plane 13 j / / ; i j 1 \ 1 \ t 1 f • I - \ \ \ f / / t t • i • > \ 4 \ \ - v 1 • i • 1 \ \ \ ! * > * 31 — - Z o The corresponding camera resolution in 2D has been measured to be in the neighborhood of 6mm in the radial direction and 5mm in the axial direction (see reference 1 0 for exact figures) although the value varies with source position radial offset (the camera response to a point source becomes increasingly more elliptical and wider as the source is moved radially outward). 1 1 I • • * . . . . • ,. For 3D acquisition on the E C A T , the 2D plane defining septa are retracted from the F O V , permitting acquisition across all opposing axial ring combinations up to a ring difference of 15. The resultant 162 = 256 . unique plane combinations are summarized by their polar (0) inclination relative to the transverse plane 0=0: 2 D vs. 3 D SENSITIVITY D I F F E R E N C E S :"' . .. • . *• . '• • .  •" • . • 1' ' '. - • . . . . . ' ; Retracting the septa from the camera aperture increases the Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING • ' ! ' • \ acceptance angle of the detector crystals. The net affect is an increase in the detection efficiency. In the particular case of the E C A T , septa retraction results in an efficiency gain from ,51% to 3.2%.12 A simple rationalization is demonstrated in the axial sensitivity comparison for a hypothetical three.ring camera (here each of the three planes is separated from the other by an interdicting black line -septa) with retractable septa : . Figure 2-3 Axial Sensitivity: 2D Vs 3D Imaging Sensitivity representative of the s. \ • / j / xxxxxx: I X X X X X X X X X X X X X X X X X ^ K K X X X X X X X X X X : II II X •/• \ ; • / • v . - / \ B X K X X X X X X X X X X X X X X X X X X : ' • r V / \ t\ A / I " . • . \ / \ \ X X X X K X X X X X X X X J C X X X K X X X X M 3D No Septa i 2D |With:Septa The improvement in the sensitivity defined as the number of counts detected in a unit time for each unit activity in the source, stems from gain in the geometrical efficiency for the camera achieved by retracting the septa. With the septa retracted, the (camera geometry dependent) out of plane acceptance angle on the E C A T 'rises from .4° to 7.5°. The absolute sensitivity when measured for a radially centered ; 8 F line consequently grows from 5.1xl03counts-sec'^MBq1 to 3.2xl04 counts-sec"1 •MBq 1. 1 2! \ Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING 2.l.4.4.SiGNAL-To-NoiSE PERFORMANCE G A I N IN 3 D O V E R 2 D The increased collection of true coincident events is mitigated by the simultaneous increase in the collection of random (defined in section 3.3.1) and scatter events (no septa to prohibit oblique cross-plane scattering). Only primary events contribute meaningful information regarding the spatial distribution and intensity of the activity source. Volumetric imaging therefore represents a trade-off between the signal enhancing benefit of the increased acceptance of primary counts contrasted with the signal degrading impact of scatter. A compact gage of the impact of 3D acquisition on image quality is garnered through the use of hoise-equivalent-counts (NEC) which incorporate the noise like impact of scatter and randoms.13 The image quality of a particular scan is given by its signal-to-noise ratio.13 This can be shown to be equivalent to the NEC parameter which incorporates noise affects introduced into the image due to the subtraction of random events and scatter count rates from the full measured count rate.13 Simulation studies on the ECAT reveal a peak NEC performance at .5 uCi/ml for the 3D imaging of a 68Ge filled source. The matching slice-oriented NEC distribution for the identical source and phantom will not peak until 2.5 uCi/ml activity.14 However, the exact performance gain is wholly contingent upon the particular size of object being imaged, activity distribution and count rate. At high count rates in large objects there is. no significant performance gain between the two imaging modes of 2D and 3D.3 '* !' - . • 12 I Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING At relatively low radioactivity concentrations (activity <.5 jiCi/ml) the increase in the primary events outweighs the competing desultory affect of scatter and random events. The net increase in sensitivity has demonstrated an improvement for data quality in scans involving 18F-labeled tracers (FDG and Flourodopa) as well as n C labeled compounds (Deprenyl, Flumazenil, Diprenorphine).7 The corresponding kinetic modeling of the active neurotransmitters and receptors involved in the uptake of these radioligands is thus also improved. 2.2. R E C O N S T R U C T I O N 2.2.1. E V E N T A C Q U I S I T I O N A N D R E G I S T R A T I O N Opposing crystal elements are connected through a coincident circuit, y rays simultaneously coregistered by opposing detector elements, within a preset coincidence time, window (-12 ns), constitute a single count: reflecting one decay event occurring along the ray connecting the detector pair. The ray thus delineated by the opposing crystals designates a line of response or LOR. Multiple counts acquired over a predetermined imaging time signifies the shadow of activity (activity ratexscan time) • • • ' • i through the scanned subject along the trajectory spanned by the LOR. The sum of counts detected in one detector channel (a single combination of two detector elements) is therefore proportional to the integrated tracer activity through the source along that LOR path. . 13 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING-2.2.2. D A T A O R D E R I N G - T H E S I N O G R A M The recorded sum of events registered along a given LOR represents the integral of activity through the source being scanned at the trajectory of the LOR. Thus each detector channel (opposing detector pairs) defines a sample of activity proportional to the integral of radio-tracer concentration through the source at the orientation of the LOR intersecting the channel. The collection of counts along a given LOR orientation thus defines the heretofore expressed angular projection: where the angular view of the projection is denoted by the angular coordinates 0,(p of Figure 2-1. The e ' - e + annihilations occur with a Poisson distributed probability, inside the source of distribution (f(x,y,z) I x,y,z-x). Therefore, event detection at the scintillating crystal occurs with a Poisson distribution bearing an intensity P at the channel corresponding to the projection (6,(1)) at position {(xr,y)=V\ in the 2D projection. 2.2 P( l) a o = J W(x I l) e,, * PSF(x 11)M * f (x)dx + R(l> The PSF is the Point Spread Function of the tomograph i.e. it's! response to a point source. The functional W(x 1I) correspond to the survival probability of an annihilation at x inside the source being detected at I . The accidental coincidences or random events (events where 2 y-rays are detected simultaneously but originate in separate annihilations) R contribute additively over the image space X to the projection space L . The Poisson-distributed incident counts are recorded at the output of the PMTs by the real-time coincidence/sorter electronics determined with a coincidence window ' The coordinate transformation is stated i n Figure 2-1. 14 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING technique15 which simultaneously assigns true coincident and delayed random coincidences to their appropriate detector channel (see section 3.3.1). The detected set of events P(xf/yl)e^ may then be ordered into a sequence of parallel ray-sums (same Q,<\>), generating shadow projections of the activity distribution through the object. This is most easily demonstrated for the example of 2D imaging where only the single azimuthal parameter <j) is of relevance. Moreover since the projection is now single dimensional, the projection coordinate system of xr,yr reduces to x-r. The transformation between the object space coordinate sytem x,y and that of the projection at § is also given in Figure 2-4: Figure 2-4: Demonstration of a Tomographic Projection and Sinogram Construction Projection; p(xr,<M Object: ffr.y) tf S i n o g r a m : s(Xr4) -**— single projection sine wqyelraced out by a points (x^ y^ ) The.fx^ histogram of the acquired data is referred to as the sinogram and is the basic format in which data are stored by the front end acquisition system. Notice that the 15 Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING . ' ., \ detected set of events P(xr,§) were assumed to give a faithful projection of the source unlike the case of equation 2.2. For 2D imaging the complete sinogram set will consist of 31 odered sinograms planes, each referencing a different acquisition plane.11 In 3D imaging, the added degree of freedom inherent in the polar (oblique) angle expands the 2D sinogram set to 256 sinogram planes in total. ' 2.2.3. S i n o g r a m Invers ion - Image R e c o n s t r u c t i o n The mathematical theory of image reconstruction from line integrals was first expounded in- the seminal 1917 paper of Radon.16 The most common reconstruction employs a filtered-backprojection (FBP) of the tomographic projections: re y ' • • 2.3 f(x,y,z) = j"d4) Jdecos9-*F2-1^>(ylii,vL)e,+ *G(v li,v l y) e i t] • 0 - V ; where P(vxr, v^ e*. is the two dimensional Fourier transform of the 2D projection and the arc length -\|/<9<\)/ corresponds to the polar extent swept by the tomographic camera. The apodizirig window function G controls the amplification of the noise' power spectrum at higher spatial frequencies - effectively restoring the power spectrum components of the image, degraded by the non-constant modulation transfer function of the camera.17 ; Although a large number of alternate but equally valid reconstruction algorithms are available18 the major advantage of FBP is its ease of implementation and speed of execution. The FBP has been extended to a practical 3D image reconstruction by the implementation of a 2D filter developed by Colsher19, and the two step Chapter 2: POSITRON EMISSION TOMOGRAPHIC IMAGING i reconstruction procedure developed by Kinahan and Rogers.20 This 3D reconstruction method is currently being used in the U B C / T R I U M F PET facility. 17 Chapter 3: Scatter 3. S C A T T E R 3.1. AQCUIS IT ION MISREGISTRATION The emitted y rays will not always transport through, the object material unobstructed. Interaction with outer valence electrons in the biological tissue surrounding the distributed activity source may deflect the photon from its incipient path while simultaneously degrading its energy. At y. energies of 511 KeV, the nearly exclusive mechanism for interaction in situ is Compton Scatter (see section 3.2). The systematic redirection of photons gives rise to an' error in the spatial assignment of events in the tomographic projections. The subsequently recorded events along the newly delineated LOR no longer intersects the progenitor tracer position; the nascent event is misregistered.. Figure 3-1 Misregistration Due to Scatter The diagram in Figure 3-1 reveals that if one of two events originating in the ' ' ' • • • ' • • ' . i • • . • annihilation, scatters into an alternate detector channel, the subsequent LOR 18 Chapter 3: Scatter desginated by / references a wholly different angular view from the incipient LOR designated by i. The event assigned the projection defined by the orientation of LORf is therefore specious. When the affect is distributed over the entire scattering medium, the resultant projections are found to be corrupted with a spatially non-uniform distribution of scatter events. The scattered events effectively contribute a superfluous (they do not retain any spatial information regarding the original source activity) background to the simultaneously recorded primary or unscattered events. For ID projections (2D imaging), the scatter events contribute roughly (depending on camera design and object composition) 10%-15% of the total number of recorded events. The additional. degrees of freedom in acquisition permitted by volumetric 3D imaging combined with the absence of septa however augments the scatter portion to a number in excess of 35%. The magnitude of the error in the consequent reconstruction will summarily increase. Thus the quantitative validity of the reconstructed source activity (its resemblance to the originating activity) is correspondingly lessened. 3.2. P H Y S I C A L O R I G I N O F S C A T T E R For human tissue whose density is comprable to that of water non Compton scatter modes, namely photoelectric absorption bear a negligible probability or cross-section for 511 KeV y-rays (Figure 3-2). 19 Chapter 3: Scatter Figure 3-2 Relative Contribution of Photoelectric Affect, Compton Scatter, and Pair Production as a Function of Energy 2 h-z UJ o LU U_ 111 o 0 n < ^ Z LU 1 DC < LU z . 4..IT Water 0.001 0-01 h i&jpton b a t t e r 0.1 1.0 10 PHOTON ENERGY (MeV) Clearly Compton scatter is the. principle mode of y-interaction. Compton scatter is a fully quantum mechanical behavior describing the energy degrading deflection of an energetic photon off a free electron occurring when E Y ~ m ec2. The probability of scattering a 511 KeV y into an angle (6) about its original trajectory due to its interaction with the, biological materials' atomic electrons is closely approximated* by. the Kleih-Nishina differential cross-section (for unpolarized photons):21 , ' : • 3.1 do(0) T e 2, E f da 2 n\.c E; E f E ; The cross section is i n strictest, terms accurate for free electron syterns only ' 20 Chapter 3: Scatter Commensurate with the angular deflection is an energy degradation for the y given U 21 3.2 Ly(e) = ((hv) 2^ ™sr.) 1+COS0 ,l+-^(l-sin0) The probability of a photon scattering into an angle (0) decreases with increased (0) at all y energies. ' Figure 3-3 Energy Distribution as a Function of Scattering Angle for 511 K e V Photons c 500.0 o o CL 400. Q TJ _ 300.0 it: c LU • I ! r -1 _ -:A'ng]e of Photon S c a t t e r i n g Compton Sca t te r ing . -o f 5 i i keV Photons Evidently the larger the camera aperture, the greater the scatter component in the projections 21. Chapter 3: Scatter The likelihood of a photon with an energy E traversing a dense medium, such as water without scatter is fully characterized by that, material's coefficient of linear attenuation [i(x,y,z;E) for the material distribution f(x,y,z). The probability of a photon originating at a. point rt inside the material and surviving to a point r2, without Compton collision is described by the probability amplitude: 3.3 P(rpr2;r^) = e x p j ^ where the coeffient of linear attenuation \i is implicitly allowed to vary across f(x,y,z). The result holds true when scatter is the sole mechanism available for removing photon flux along the original photon trajectory (cross-section for photo-electric is negligible for 511 KeV y's in water ; see Figure 3-2) i.e. the scatter component at r must be equal to the negative instantaneous change of photon flux at F: 3.4 ^ P ^ - W O Given an original activity distribution A(r) within the material distribution f(r), the number of events recorded at the detector element at r / (assuming perfect detector efficiency) over a scan time T will be dependent upon the attenuating material distribution and path length between activity source and detector: N(r') = Jdt'J df* K J *g(F-r)*A(r)*P(F'„r ;E) T rov 4 * \ * . The aperture function g(r'-r) delineates the detector geometry accepting events; the larger the aperture the greater the number of recordable events. The events (r'-r) 3.6 T FOV ^"r 22 \dt f df*-^r4T*A(r')*g(r'-r)*(l-P(r',r;Ey)) ' ' Chapter 3: Scatter do not reach the detector surface at / . The portion of emitted y's, represented above, may undergo scatter away from the originally prescribed trajectory at either a single or possibly multiple site locations. As previously demonstrated, the correspondingly delineated count is either fully lost (scattered outside detector aperture) or if detected, misplaced into a spurious projection ray-sum. The assumption that the measured projection is equal to the shadow of activity therefore falls short. 3.3. M A N I F E S T A T I O N S OF S C A T T E R B A S E D E R R O R The scattering of the photon will result in several, possibly coincident, aliasing affects. If the scattered photon is registered with its non-scattered partner, the . assignment of the consequent event to a specious projection bin will effect an error of misregistration. Alternatively, a scattered photon may contribute a random count to the projection ray-sums; if the scattered photon is detected in coincidence with an unrelated, possibly also scattered photon, the subsequent event bears no correlation to the underlying activity source distribution. Furthermore, independent of whether the scattered event is in fact detected and so misregistered, a primary event is always lost from a projection ray-sum. This implies a depletion or reduction in count number from the collected projection ray-sums. Thus the imaged activity source will appear attenuated in strength. Failing a post-collection normalization of the ray-sums for object attenuation, a quantitatively accurate reconstruction of the original activity proves impossible. 23 Chapter 3: Scatter While all these heretofore listed artifacts originate in the scattering process, they impact differently on the projections and so are usually handled in separate ways. Although the issue of scatter represents the focus of the present work, the attended issues of attenuation and random events will be discussed briefly due to their obvious connection to scatter. 3.3.1. R A N D O M E V E N T S Two photons originating in separate decay events will contribute a random event R to the projection ray-sum if detected in coincidence. The rate of random coincidence is conditional on the joint rate of photon detection within each of the ray-sums' detector pair defined bins (indexed by i,j) and the coincident time window v. 3.7 ' R,, =2l;XNJ The relative number of random events may be succinctly expressed as a fraction of total recorded events. For the ECAT these numbers may become very significant: 10.6% random fraction per |i.Ci/cc of activity in 2D or 76% / ^iCi/cc in 3D within the FOV.22 On the ECAT the random events are processed simultaneous to the coincident events by the real-time sorter but through a 128ns delayed coincident circuit (implying that the coincidentally recorded events originated in different decays). The offending random contribution is subsequently excised by subtracting the measured random events from the measured coincident data, sinogram bin by bin. The procedure of subtracting a statistically noisy random events measurement from a noisy emission scan will result in the summing of the total image noise. The affect of additional noise 24 Chapter 3: Scatter due to random event subtraction has not yet been explored in 3D PET(although it has in 2D). 2 3 Iterative reconstruction methods avoid the problem of noise propagation in the random removal program by the inclusion of random events as an argument in the generating recursion expression. 3.3.2. A T T E N U A T I O N The deflection of a photon outside of the detecting solid angle will deplete or attenuate the measured count-sum in the channel originally intersecting the emitted fs. The loss of events to attenuation reveals itself as an effective error of omission. The omission is not uniform across the projections, but is instead dependent (for medium based scattering) on the depth of scattering material separating the detector element and emitting source; for PET where Ys are measured in coincidence, the length of the whole ray intersecting the object determines the ultimate attenuation in a particular detector channel. As observed in equation 3.6, the attenuation increases with increased path length traversed by the emitted photon pair. Thus, two identically strong activity sources at disparate depths inside the scattering material will project identical activity distributions onto the same channel if they fall on the same axis. Still in PET, the attenuation of the projected activity will depend upon the particular density distribution and geometry of the objects being scanned. Hence a single activity source may be viewed differently by different detector channels separated by some angular arc through the tomograph. Several methods exist for correcting for the lost counts in a projection. They all rely on some estimated or measured attenuation map over all the coincident rays . . . 25 . Chapter 3: Scatter intersecting the FOV. Each sampling bin in the projection corresponds to a L O R at the inclination of the projection view displaced from the radial center by |r| = Jxf+rf . Each bin can thus be independently corrected by reversing the equation for attenuation along each LOR. 3.8A/(F)= \dt \ dJ'*-^^*g(7-7)*A(T)*P(T\T\E)*zxr>U T FOV 47t |F-r | [i The attenuation correction factors (i can be estimated in several ways. One method is to calculate them directly by assuming the attenuation factors to be those of water. The problem of attenuation correction therefore reduces to estimating the path length of a given L O R through the subject. This method, referred to as calculated attenuation (intended to both avoid exposing the patient to any additional ionizing radiation and introduce further noise associated with the added scan), makes use of thresholding algorithms applied to the imaged object in order to determine its contoured boundary from the emission scan and thus enabling the derivation of the sundry L O R path lengths. The principal problem with the method lies in its sweeping assumptions regarding the attenuation factors. While brain density is roughly uniform (with an attenuation density approximately equal to that of water), it .becomes discontinuous at the brain/cranial interface. Moreover, the inherent assumption of uniform density is completely invalid in torso imaging due to both bone deposits and air deposits (lungs). A common technique is to perform an ancillary attenuation transmission scan before the emission scan using any number of positron emitting sources (three rotating opposing 6 8Ge rod sources forming an equilateral triangle on the E C A T ) outside the 26 Chapter 3: Scatter . -scanned object, near the detectors. Traditionally the implementation first involves recording a transmission scan using the external source with no object in the FOV. This is termed a blank scan. Then the subject is placed in the F O V and scanned with the same external source. The attenuation correction factors are then calculated from the comparison of the two scans normalized for duration. Before the projections are corrected for attenuation, a scatter correction will be necessary though. Since the L O R chord to which a scattered photon pair is assigned is spurious, the offending projection count should not be scaled. 3.3.3. MISREGISTRATION . For a typical brain tomograph 60% of the scattered events are scattered less than 30° . 2 4 However, even small angle scattering will result in noticeable errors in the assignment of LORs. The gross affect of scatter based misregistration is to introduce a low spatial frequency background noise into the projections. In an idealizing simplification, the misregistered events S, here on referred to as scatter, will be regarded as a spatially varying noise added onto the spatially varying primary (i.e. unscattered) true count distribution P: 3-9 ' M( l x ) l y )=P(l x , l y )+S(L,L) where M is the measured projection generated by primary events P and scatter events . 5. In order to demonstrate the salient properties of scatter, consider a 6 8 Ge point source, secured centrally within a 10 cm radius cylinder filled with water (see section 5.3.1 for details). The imaged source is* centrally positioned inside the F O V and will 27 Chapter 3: Scatter therefore project to the same position in each angular projection. Subsequently, the individual angular projections may be summed (in this case over all 192 angular views) or mashed in order to statistically augment the projection measurement. The mashed projection, formatted into a 2D planagram (see section 5.4) is displayed below. For clarity, the vertical axis is logarithmic. Figure 3-4 Mashed Projection of a Centralized Point Source ( 259,760 counts) The central peak of the mashed projection, corresponds to the projected point source position. The spread in the peak is attributed to the PSF of the imaging camera. It reflects the blurring effects introduced by the finite detector size; positron range following positron decay; non-parallel flight of the photon pair; and variable photon penetration into off-center crystal elements.25 It is noteworthy that it differs slightly . 2 8 ' Chapter 3: Scatter between the two projection axis referred to above as r, z instead of yr (since yr is identically the axial z axis in a planagram). Outside the projection peak the measured events are exclusively attributable to scatter. Underneath the projection peak, the events are a combination of scatter and primary events. Beyond the point source, no further activity sources exist. The counts in the projection lying outside the projection peak are misregistered events that were deflected from different projections - their number decreasing with increased distance away from the projection peak. Since the fan of coincident crystals extends to cover the whole F O V in PET the scatter is not restricted to the projected region containing scattering material, but instead extends over the whole F O V smoothly - truncated in the above diagram to 40cm in diameter radially by 10.6cm axially. Several factors will affect the distribution and size of the scatter. The magnitude of the scatter is proportional to the initial activity present inside the object. The density distribution of the material object and in particular its relation to the source activity position will determine the distribution of the scatter. 3.4. S C A T T E R I M P A C T The scatter fraction (SF), projectionbins _A • : • 3.io S F = . , . .•: projectionbins - . I K •''} 29 Chapter 3: Scatter an index conveying scatter contribution to the measured projection, increases from 16% to 41% for a 68Ge line source centered inside a 20cm diameter right-angle cylinder, when imaged in 2D vs. 3D mode respectively, on the 953B. 4 . S C A T T E R C O R R E C T I O N Correction may take either the more approximative form of compensation or else full removal. This requires that the scatter events be either rejected during collection using physical interdiction or energy discrimination or otherwise estimated and then subtracted off. Volumetric imaging precludes the use interdicting inter-slice septa barriers. Energy discrimination has limited use due to the limited energy resolution of the detectors. The following chapter will review the common post-collection (i.e. software) scatter correction methods.. Particular emphasis is given to the methods broadly categorized as integral transformation (convolution subtraction, deconvolution) although energy based methods will also be discussed for the sake of completeness. 4 . 1 . S C A T T E R S U P P R E S S I O N The presence of scatter in the projection may be accounted for by adopting an effective attenuation coefficient along each LOR, that incorporates both the presence of attenuation and scatter. Effective attenuation values that include scatter have been variably measured between .085 and .088 cm'1 in contradistinction to the usual attenuation factor of .097 cm'1 for attenuation in water.26 The effective attenuation ' 30 . Chapter 3: Scatter . value will, due to the inclusion of camera dependent scatter performance, vary according to camera geometry. The method however has never been widely adopted because of its failure to remove the scatter events from the measured projection set. This results in an excessively smoothed image (due to superposition of smooth scatter background) and. so decreased contrast. The decreased contrast hinders quantitative analysis of the reconstructed activity distribution.26 • 4 . 2 . S C A T T E R S U B T R A C T I O N Correction for scatter subsumes the additive noise model for the projection ray-sums. A scatter distribution in the projection is estimated from the observed data. The determined scatter distribution is then subtracted from the total measured projection ray-sums, bin by bin. Several methods of estimating the scatter have been proposed. In broad classification they come under the aegis of either energy window sampling, or integral transformation. 4 . 2 . 1 . S E G M E N T E D E N E R G Y W I N D O W M E T H O D S All scattered events are degraded in energy (see Equation 3.2). Optimally, it would be ideal to excise the scatter from the acquired events by energy discrimination in the detector. However, crystalline detectors in general, and BGO based detectors27 such as the ECAT in particular have poor energy resolution (about 25%)/ This fact limits the use of energy discrimination as an exclusive means of scatter The exact energy resoluation varies for each B G O crystal element inside a detector block (see [27]). • 31 Chapter 3: Scatter suppression. The lower level energy discriminator setting as a result becomes a compromise between rejecting low energy scatter events without losing poorly resolved primary events. On the ECAT 953B, this concern has led to the use of 250 KeV and 380 KeV for lower settings in 2D and 3D imaging respectively. Alternate use of energy information for scatter correction has also been investigated. The schemes generally fall under the rubric of energy windowing methods. An illustration of the method is effected by the example of the Dual Energy Window acquisition.28 • Segmenting the energy spectrum into non-overlapping sections or windows permits the simultaneous acquisition over independent portions of the spectrum. This method is appealing since it takes into account the scatter originating from radioactivity outside the FOV. It does however require acquisition of larger data sets and still suffers from geometrical approximations.29 Events acquired in each of the lower and upper windows Mumr(lower) can be expressed as a sum of scatter S and primary P events: _ - • ' ^^upper ^upper ^upper M = S, + P lower '-'lower 1 lower The two equations can be solved simultaneously for the primaries in the upper window:28 4.2 P =\ \*M - J — • • —= \*M, \ " /^upper ' / "upper J V /^upper. / 'upper J The various upper-lower count ratios are derived from a line source measurement in a water-filled phantom and in air. However, the method does not take 32 Chapter 3: Scatter into account possible spatial variations of the ratios measured with the phantoms. This method has been found to perform fairly well for objects of homogenous density, but poorly otherwise.30 Another proposal has been to place two energy windows of differing widths below the photopeak.31 Comparing the ratio of collected events in the two windows first during a uniform source calibration scan and then during a patient scan will return a modifying parameter associated with the degree of deviation from source uniformity in the. patient. This may then be used to compensate the primary estimate under the photopeak for source inhomogeneity in the patient study. A.newly expounded method32 eliminates the previous problem of dependence on phantom measurement without relying on multiple lower energy windows. Instead, the method involves estimating the trues (as opposed to the scatter) underneath the photopeak by using acquisition information from the upper spectrum of received counts. The events lying at the outer right hand extreme of the observed spectrum correspond to primary events alone. Their location to the right of the 511 KeV photopeak mean, stems from poor energy resolution inside the. detector. Establishing two energy windows, one.inside the other, the primary events in the window covering the photopeak from Ej to E 2 can be found based on a predetermined ratio function between the primary events in this window and those in a narrow upper window between E u p (peak for 511 KeV y< E u p < E2) and E 2. The major drawback of this method is the fact that it estimates the primary events from a smaller and thus statistically noisier sample. It is also dependent on the stability of the energy 33 Chapter 3: Scatter discrimination threshold. Energy windowing methods have an advantage over other scatter correction schemes. One is that they can be used to reject scatter from outside the FOV. In fact, they can be used to eliminate all energy degraded events regardless whether they originate in the object, detector or elsewhere. While this issue is relatively relevant in neuro-imaging3 3, it does become a significant concern in whole body PET where activity is distributed through out the circulation.34 Conversely the methods demand larger data sets and possibly, new energy gating electronics capable of rhultispectral acquisition. Since the primary events are obtained by a subtraction of measured events from total events, their statistical noise content is compounded. Filtering techniques are being used to decrease the noise content in the unscattered events. 4 . 2 . 2 . I N T E G R A L T R A N S F O R M A T I O N The last most widely adopted method of scatter compensation renders a linear systems perspective to the interaction of the emitting source with the surrounding dense material. In the usual formulaic statement, the scatter constituency in the projection S(lx,lJ is expressed as:: 4.3 S(lx,ly)= jP(lx',ly')®h(lx',ly';lx,ly) projection The kernel h, is a scatter distribution function giving the projected scatter amplitude at the position lxly for a point source at the position lx',ly'. Since the scatter response is presumed to act linearly with activity/the following must hold for two independent 34 Chapter 3: Scatter activity sources indexed by 1 and 2: 4.4 s(ix,iy) = slaxjy) + s2(ixjy) = •' J{/?(/;,v) + / , 2(/;^ ,)}®w x',/;;/ x,g projection Underlying the linear model of scatter response is the assumption that the scatter is separable from the primary true data. i.e. the projection distribution is the additive sum of a true (non-scattered) count distribution and a scatter count distribution. It is also assumed that changes of the source position along the projection direction do not affect S. Hence the measured' counts Md^lJ in the projection may be expressed in terms of S(lx,lJ and P(lx,lJ • 4-5 . , M(l x , l y )=S(l x , l y ) + P(lx,ly)-This model is analogous to the usually adopted model of signal noise as additive, wherein the role here of signal is the primary count distribution P while the scatter evokes a semblance to noise. In this sense, the scatter while not statistical in origin admits a similar treatment to that,of unwanted noise. O n the condition that a distribution for the scatter is derivable, then the true count distribution may be extricated. The complete solution will be 4.6 P( l x , l y )=M(l x , l y ) -S(L, l y ;P) The most glaring failure of the heretofore expressed method (Eq. 4-3) is the endemic inability to isolate the true count density. The solution for the true count density will involve an integral equation. The tractability in recovering P from arbitrary M is not immediately clear. Several simplifying mechanisms exist to generate an approximate solution. They may be broadly classified as either that of convolution subtraction35, or deconvolutidn3 6 , 3 7 Both methods share a number of common features, most notable being the 35 Chapter 3: Scatter requirement of estimating an analytical scatter distribution in the projections. From herein the methods diverge. DEGONVOLUTION The method of deconvolution begins with the normative relation, already elaborated, of scatter being the convolution of the primary count distribution with that of a scatter. kernel h. In its simplest incarnation the kernel is heuristically selected to be a position independent mono-exponential38 or similar function bearing the property of even symmetry (raised cosine39; gaussian). The choice of function is based on its observed consistency with the available data. In 2 dimensions (2D projections) the kernel may be on first approximation prescribed as a radially symmetric, position invariant function on the projection xr,yr coordinate system.36 However, the scatter profile will rarely be either radially symmetric or even. It is in fact skewed with the degree of eccentricity dependent upon the location of the source inside the projection itself - and so position dependent. kernel form across the whole projection space. Shift invariance in the kernel facilitates the recovery of the true count distribution through Fourier Inversion Methods: r For reasons of simplicity the kernel has been presumed to be spatially invariant i.e. 4.8 36 Chapter 3: Scatter where F corresponds to the 2D Fourier transform41. C O N V O L U T I O N SUBTRACTION Alternatively the primary activity distribution may be calculated through an iterative convolution method where by the initial estimate of P as is required by the method of deconvolution, is approximated on a first pass as M , which may then be used to recursively generate improved estimates of P based on the intermediate estimation of the scatter. The iterative procedure is concisely expressed as follows: 4.9 M - S ( P n i ) ; = P n ; P ° '= M The recursion may be repeated until a stopping criterion is met (such as no meaningful statistical change between consecutive estimates of P. CONVOLUTION SUBTRACTION There exists an alternate possibility for the recovery of P. This is based on the ansatz that the scatter component has a valid and demonstrable formulation as a linear functional of the full, measured projection distribution. 4.10 S'(l x,l y)= j M G x',I y ' ) ® h , ( l x ' , l y , ; l x f L ) projection No prior estimate of the primary count density is necessary. in implementing the correction. Avoiding the iterative correction eases computation (fewer calculations) and removes the necessity of establishing a stopping criterion for terminating the 37 Chapter 3: Scatter calculation. The method espoused by Eq. 4-8 is termed convolution subtraction. The method of convolution subtraction was first proposed in 1982 by M . Bergstrom42 in order to avert the inverse problem encountered in the method of direct deconvolution. Originally it was applied to the correction of scatter events in slice-oriented PET. The scheme postulates the existence of a kernel h' whose action on the measured projection data M, returns a scatter distribution S' indistinguishable from the actual scatter profile S. The primary thesis of this method is that there exists a linear, correlation between the measured count distribution and the component scatter enmeshed within it. The previously proposed scatter kernel h (see is no longer suitable in the current paradigm, as its action was restricted to the primary count density. The implication of the hypothesis is that the modified kernel K (bearing the same functional form as h) must satisfy the constraint equation: 4.11 P ® h = M ® h ' Although the method of deconvolution is not founded on any known physical principals, it proves valid when applied over the full projection data set.43 The method of discovery for h' follows from the preliminary identification of h. The subsequent form of K is assumed to share the same functional form as does h but owning new parametric values. The new parameter values are determined through a second minimization of the convolution product M®h' to the observed scatter distribution. It turns out this procedure depends on the first discovery of h (see section 7). The accuracy of both correction schemes is limited by the preliminary accurate determination of h - this measurement becomes the linchpin of integral transformation 38 Chapter 3: Scatter correction paradigms. A robust kernel h must account for the variable source distribution in the projection as well as the commensurate eccentricity in the distribution - which may vary differently in the two coordinate-ordinate directions inside the projection. Hence the statement of the functional kernel model must bear sufficient parametric freedom to accurately mirror the scatter data in the measured projections. Satisfying this in conjunction with the further assumption on the form and extent of the primary distribution (having previously modeled the scatter) admits the possibility of correctly extricating the unknown position variable kernel h . The strength of the convolution subtraction correction will be limited by the accuracy with which the kernel h is known. While it is in fact observed that sensitivity in the functional form of the kernel will not play a decisive role in the performance of the correction,43 the magnitude of the function will strongly influence end performance. In brain imaging where the object is relatively uniform in composition it has been suggested that the scatter kernel will not vary significantly across the projection space.36 The proceeding chapter will examine methods used in solving h with the end goal of implementing an iterative convolution subtraction . 39 Chapter 5: Method for Determining The Kernel h 5. M E T H O D F O R D E T E R M I N I N G T H E K E R N E L H 5.1 . T H E O R Y The fundamental ansatz of convolution subtraction (see maintains that the scatter distribution S in the projection coordinate system, spanned by xr,yr herein referred to as lx,l is expressible as: 5.1 • s ( u y ) = p ( W < 8 > n ( u * ) However, following the collection of coincident counts, only the total measured count distribution is immediately accessible; the contributing primary P and scatter S subcomponents are inextricable from the total - at least in the usual case of a distributed source. Conversely, if a localized point activity source with total activity A is used such that it's projection is equal to: 5 2 Mil l ) - \ A l*=l*'ly=ly x , y [0 elsewhere the problem of separating the scatter from the total count distribution simplifies. The counts outside l'x,l' must.originate in scattering processes, since they do not arise from the direct projection of the activity source P. Based on the possibility that S is separable from P in the projection, the kernel response of the scattering object to the initial known source activity becomes quantifiable. The validity of this claim is revealed below. 40 Chapter 5: Method for Determining The Kernel h 5.2. D E R I V I N G R E L A T I O N S H I P B E T W E E N S A N D H . Assuming an idealized point source inside a scattering medium and an idealized detection system (thus neglecting blurring introduced by finite resolution), the matching scatter response^  effected by the kernel h, will be given by: 5.3 S(lxJyX,l.=A.h(Ufj;j;) It is therefore possible to estimate the functional form of the scatter kernel from the scatter distribution in the projection bins not containing the source projection. Due to the finite resolution of the tomograph the source profile will have a width different from zero associated with it. The functional form of the scatter kernel can therefore be evaluated from projection regions far enough from the source profile. Subsequently by integrating both sides of equation 5.3 and rearanging the integrated kernel is seen to be equal to the scatter fraction: ' - • * 5.4 jh(lx,ly;lx',ly')=."**- A , • ( projection The final result above establishes that the kernel h is derivable from the scatter S, according to a simple paradigm. The experimental execution of this paradigm is examined in the following section. 5.3. I M P L E M E N T A T I O N A complete characterization of the scatter response kernel h over the entire projection, reduces according to the previous prescription, to the preliminary characterization of the scatter effected by. a point source at disparate positions within the projection. The characterization of the scatter for a projected source position If,I' is 41 Chapter 5: Method for Determining The Kernel h accomplished by fitting a parametrized model function ^modeil; • d >l >fli'ao >• • • ' a )onto the observed scatter data. The parameters a1-an introduce sufficient degrees of freedom into Smodel so that it may accurately estimate S. The number of parameters n and hence the ultimate profile of Snadel is determined by examining the salient features of the scatter distribution (see Section 6.1). Upon the selection of the analytical form for the scatter distribution fits are performed of Snmde, onto the scatter derived from each projected point source position lx',lz'. The recorded scatter samples fix the parameter settings a1-an of Smodel for each point //,//. The correspondingly assembled set of parameters now indexed by lx',lz' thus fully describes the form of the scatter as the point source position is varied across the projection space. By virtue of the connection established in equation 5.3 between the kernel h and the scatter, the kernel is similarly everywhere discerned (up to a scaling factor equal to the integrated projection count density). 5 . 3 . 1 . M A T E R I A L S The point source idealization, heretofore used to develop the program for determining h, was realized in the course of the experiment through the use of a 68Ge point source encased inside a cylindrical steel sheath 3mm long. Since the focus at the UBC/Triumf PET facility is brain imaging exclusively, the object used to represent the scattering material was water enclosed within a 20cm diameter Plexiglas cylindrical phantom, which is the standard phantom used in this type of study. 42 Chapter 5: Method for Determining The Kernel .h The use of a 20cm diameter water filled phantom is not universally accepted for modeling the head. 4 4 The typical male or female human head is elliptical in shape averaging 19.6cmxl5.5cm and 18.5cmxl4.5 respectively.45 However, the scatter distribution has been shown not to vary significantly as a function of object. size.39 Furthermore, the cylindrical phantom offers symmetry benefits. The object will bear reflective symmetry in the axial and transverse plane about the central axial and radial axis. Ultimately this will be useful in estimating the scatter kernel parameters fl3-fl„. over the whole 2D projection (the parameter model must conserve the system symmetry; see section Z.3.1). The use of water as the attenuating medium is accepted based onthe fact that non-osseous tissue has a density very similar to water. A supporting holder for the point source was constructed from a hollow plastic brace sandwiched between two Styrofoam cushions - which were subsequently wedged at the base and ceiling face of, a 20cm diameter- 26 cm long polyethylene cylinder phantom. The selection of low density plastic and Styrofoam as construction material was demanded by the requirement that the support not attenuate the emitted Y's.; Furthermore, the plastic brace permitted easy displacement of the point source along its major axis. The entire apparatus is displayed below: 43 Chapter 5: Method for Determining The Kernel h Figure 5-1 Source, Holder and Phantom • Phantom 5.3.2. S A M P L I N G P R O C E D U R E The s amp l ing prescr ip t ion adopted here in w i l l be to take three axia l ly d i sp laced measurements of an 8cm r ad ia l ly offset 6 8 G e poin t source. This w i l l generate 3*96 angular projections each p ro f i l i ng an alternate source pos i t ion ins ide the projection. Each subsequent 2D projection is identif iable b y the projected source pos i t ion / / , / / w i t h i n it. The m o d e l scatter funct ion , n o w indexed b y lx',lz', is fit to each projection. This p r o g r a m w i l l generate a characteristic set of m o d e l parameters indexed to source pos i t ion . The f inal result is a look-up table of scatter responses (fully descr ibed b y their salient parameter settings), where each element i n the table corresponds to the response * The 192 angular views are summed by 2 44 Chapter 5: Method for Determining The Kernel h of the object to a source at that projection position. Each time a fit is performed for Sn at the projected source position / / , / / , a single set of fit parameters at-as is generated. 5.4. P L A N A G R A M C O N S T R U C T I O N The scatter correction is not performed directly on the sinogram of acquired events. Prior to analysis, the raw sinogram data are reformatted into a set of true 2D shadow projections termed planagrams according to the method expounded by Bailey et al.35 The resulting structure is delineated by a metric consisting of the radial lx coordinate-ordinates 128 bins long and the axial ly coordinate-ordinate 31 bins in extent. The methodology of the rebinning program is summarized in Figure 5-2: Figure 5-2 Sinogram Rebinning Into Planagrams 256 , 1 8 0 - - 1 0 - 0 • 1 0 . 3 0 1 3 < . ' : • 2 0 r ( cm) 2 0 - 5 . 3 z ( cm) 5 . 3 - 2 0 r ( cm) 2 0 S e l e c t 6 - 0 D i r e c t P l a n e s O n l y 45 Chapter 5: Method for Determining The Kernel h The planagrams comprised of planes characterized by a ring difference greater than 0 are not assembled. In this implimentation, the scatter distribution doesn't vary significantly across the sweep of a small polar angle.36 For the E C A T , with an axial width of 10.6 cm and aperture diameter of 76 cm, the maximal polar range is 7 .9° . The measured activity distribution due to path length differences through the object across different polar views is negligible. Before planagram formation from sinogram data, the sinograms are scaled by a relative efficiency factor to account for detection efficiency variation in the detector structure (see section The characterization of the scatter and commensurably the kernel h is thus performed on the planagram projection space. From herein the primary count distribution is easily calculable (see Eq. 4.6) over all the planagrams. The result is then reconstituted into the original sinogram form for reconstruction processing. 46 Chapter 6: Implementation 6. I M P L E M E N T A T I O N 6.1. S E L E C T I N G S M 0 D E L 6.1.1. P O S I T I O N V A R I A B L E S C A T T E R R E S P O N S E I N T H E P R O J E C T I O N The scatter response h must be fully characterizable in the projection for the correction scheme to work. Differences in scatter response in the projection have their origin in the variable distance traversed by the scattering events when the source is displaced inside the scattering-attenuating medium. The types of variation are separable into three classes. Given an angular 2D projection of an activity point source inside a scattering medium-a scatter profile change can be effected by displacing the point source along the L O R (chord AB) in the Figure 6-1, perpendicular to it in the transverse direction (chord DC), or perpendicular to it.in the axial direction (chord FG). 47 Chapter 6: Implementation Figure 6-1 Chord Specification In Object Space The impact of each variation to the ultimate scatter profile is detailed below. M O T I O N A L O N G L O R When the point source is displaced along the LOR, the scatter profile (ID transverse slice parallel to LOR through 2D projection) may on first approximation be unchanged due to the coincident method of detection in PET. However, at source positions near the edge of the phantom, phantom geometry may be expected to alter the scatter slightly in comparison to a centrally positioned source. In 2D slice-oriented PET this is indeed observed.42 This result is confirmed for 3D PET by 48 Chapter 6: Implementation - . the analytic simulation work of Barney for a radially displaced point source inside an 8.5cm radius water cylinder.46 Below is a cross-section . through the 2D scatter projection generated from consecutive scan simulations at Ocm, 3.5cm and 7cm radial offsets. It is evident that the projection becomes more.peaked as the source position inside the.phantom nears the edge of the phantom. Figure 6-2 Simulation46 18 16 0} 14 § 12 10 3 4 1 c ' ' ' ~r—'—r—~~1—— 1 L CO, 0 ,0) — - 0 * 5 , 0 , 0 ) — , r.. ; 17 ,0,0) — • / \ * -1 1 1— _ J 1 -20 -IS -10 - 5 i o 15 ao fx' *W Both the total scatter as well as the scatter fraction increase a small amount. The relative increase in the scatter between a source at the phantom center and at 7 cm is 6 0/ ; - • 1°'- ' A small variation is confirmed empirically in the plariagram domain - of primary interest to this thesis. The comparison was performed between first a plana gram. (259,760 counts) image of an axially and radially centered point-source, in the manner of section 5.3.1, and second--a planagrahv(24,229 counts) projection of an axially centered and radially offset source (selected from the sampled planagram set of section 5.3.2). Since the centered point source projected to nearly the same projection position in all views, the counting statistics were improved by adding together Chapter 6: Implementation (mashing) all angular planagram projections. The net result was displayed in Figure 3-4. For the sake of direct comparison, the contrasted planagrams are normalized to their respective total number of counts. Along the transverse coordinate of Figure 6-3 (displayed on a vertical log plot) the profile of the radially centered and offset source nearly coincide as is expected (the lx ("r") offset between the two images is due to the lack of a planagram sample corresponding to the projected source position of the central source): Figure 6-3 Transverse Comparison 4.0x10" c o o E 3.0x10 o cn c u CL -4 o -4-' o 2.0x10. "a CD _N "5 E ° 1.0x10 - 4 c O o I ; MI ii ; ; i .I / \ : . ;' W '/'! ! "ii !! • i ^ i • • • !! .j f : : : !,y; i : l ! I: •l, i , ! ' i /V . A-ii:/: v. - J — J — — _ ± — i : i—i L _ -20 -10 10 • ,• o r (cm) r(\eft)= 0.46,8750cm: r(YiahO = 0.781250cm: z= 0.00000cm 20 The dashed and dotted lines in Figure 6-4 correspond to the following imaging configurations (as viewed transaxially): 50 Chapter 6: Implementation Figure 6-4 Imaging Configuration 'dotted' "dashed' r (cm) r -20 20 r (cm) -20 20 B (.781cm,0.000crri,7.34cm) ly • Ix /' • .(.469cm,0.00cm,0.00cm) V (lx,ly,n) However, like, the simulation results of Barney, there is revealed a difference in magnitude when the scatter fractions are compared. In order to assay the scatter, a model function of .the form expounded in section 6.2.2 (n=2) was fit to the planagram data according to the method of 6.3.1. This revealed a scatter fraction of .412 and .480 for the radially offset and non-offset planagram projections respectively - a difference of 7% consistent with the simulated data. Were the projected scatter responses to differ greatly when a source is displaced along its L O R in the object space the proposed integral convolution scatter correction becomes impossible to implement. No method exists for uniquely identifying a single projection to a definite progenitor source position. Herein the kernel h will be assumed to be equal along all source positions along the same LOR. 51 Chapter 6: Implementation D I S P L A C E M E N T P E R P E N D I C U L A R T O T H E L O R ( T R A N S V E R S E D I R E C T I O N ) Alternatively the point source may be displaced perpendicular to the L O R in the transverse direction. This configuration is summarized in Figure 6-5: Figure 6-5 Configuration for Transverse Comparison "dotted' "dashed' r (cm) r (cm) -20 20 -20 20 c JD (7.34cm,0.00cm,0.00cm) Ix (.469cm,0.00cm,0.00cm) N s (lx,ly,n) The corresponding transverse cross-section is revealed inFigure 6-6 (on a log scale): 52 Chapter 6: Implementation Figure 6-6 Transverse Cross-Section* 4.0x10 ;-4 c O o E 3 .0x10" o c D -2 2 .0x10 4 O -E c 1 . 0 X 1 0 " 4 -1—I—I—I—1—I—I—I—I- I . i ; — i — r — i — i — ^ 1—i—p— i—i—i—i—j-: I 1 - i — i — i — r -: l,. H i I II I i V l' i 1/ V / " 1 ' _L J — i — i i i i _ -20. •-10 0 10 r (cm) r(lefV>= -7.34375cm: rCriahO= 0.468750cm: z= 0.00000cm The planagram from which the offset cross-section is derived contained 147,179 counts. It is therefore plain to see that the peak position corresponding to the projected activity source will shift in its position along the projection. More significantly the scatter distribution gains an asymmetry since the scattering material on either side of the source, as viewed from the projection, will no longer be equal. In the simulation studies of Barney, displacing the source in the transverse plane perpendicular to the LOR, resulted in the skewing of the scatter profile. As the point source is displaced towards 2C "r(left)" refers, to the projected lx' of the dotted source projection and "r(right)" to the dashed source projection position. 53 Chapter 6: Implementation the edge of the phantom the scatter distribution evolves a compression on the side nearest the edge (see Figure 6-6). A candidate function used to model an arbitrary scatter distribution must therefore contain sufficient parametric freedom to reproduce the observed asymmetry. In practice, the freedom is enshrined by the inclusion of position sensitive (in the . projection) titrable parameters within the model scatter function (see 6.2.2). Unlike the problem of scatter variation for positions along the same LOR, the situation entertained here is completely describable in the projection space - no additional assumptions are needed. Accordingly, the scatter function may be indexed by source position., AXIAL VARIATION PERPENDICULAR To L O R The same compression is evidenced in the scatter tails when the source approaches the transverse edge of the tomograph. The scatter variation due to source shifts in the planagram axial / direction will be less pronounced than in the transverse direction but identical in profile. Hence given the following configuration: 54 Chapter 6. Implementation . . Figure 6-7 Configuration for Axial Comparison* "dashed" "dotted" r (cm) r (cm) -20 20 (8.469cm,3.71cm,8.88cm) Ix (.469cm,8.88cm,8.88cm) (lK,ly,h) The axial cross-section through the resultant planagram is: Figure 6-8 Axial Cross-Section 0.1 000 F g"-0.01 00 CX" o 0.0010 0.0001 I \ /I /1 . /1 . i / I i J _ / V / A ' -LA I I . S i / _i__vL -4 -2 0 • 2 4 ' r (cm) r= 0.781250cm: z(lef.t~)= 0.00000cm: z(YiahO= 3.71 250c'm • Mot ion along n corresponds to motion along the L O R . . 55"' Chapter 6: Implementation When the point source is shifted axially, the peak position is similarly shifted in the projection. Furthermore, increased scatter count and eccentricity are both observed, as in the parallel affected transverse case. The importance of the finding is that the axial scatter may be treated with the same, model scatter function as was the transverse scatter distribution in the planagram. The proceeding chapter will examine the model scatter function selection process. 6 . 2 . D E T E R M I N I N G T H E M O D E L S C A T T E R F U N C T I O N 6 . 2 . 1 . C O N S T R A I N T S O N C A N D I D A T E F U N C T I O N The criteria made of the scatter function is that it accurately account for the most obvious features of the measured scatter distribution. The candidate function must bear enough freedom in the form of tunable parameters to render the most salient features of the scatter distribution i.e. the amplitude, eccentricity and peak position. As already discussed, the scatter will exhibit elliptical symmetry in the 2D projections. Moreover along each direction there must be an accounting of the deviation of the scatter tail due to the off-center displacement of the source. At minimum, a five parameter function (one for amplitude, and two each for eccentricity and peak position along each planagram direction) will be necessary if separability in the lx and / projection coordinates is desired. ^ 56 Chapter 6: Implementation • 6.2.2. C A N D I D A T E S Previous efforts at deriving a convolution scatter kernel have emphasized exponential models either spherically36 or elliptically symmetric.38 Other researchers have maintained gaussian47 or cosine39 models. Fundamentally the choice of model function is determined heuristically. There does not yet exist any analytical derivation that will give insight into the form of the scatter function. The model herein proposed supports the coordinate separability while providing the minimum five parametric degrees of freedom. But rather than commit to a single analytic form i.e. gaussian or exponential it tests six analytic possibilities. 6.1 SmM(K,ly) = ^*{zM-a2*(lx ~ /,)") + a, Hlx-lx')Y {exp(-a4*(/z - / / )") + a5 *(Z, - lz' n =1.0;1.2;1.4;1.6;1,8;2.0 The model parameters a3 and as generate an asymmetric linear shift in the scatter tail slopes centered on the peak position in each of the lx and / y coordinates respectively. The multiplier terms dx-lx') and dy-l') respectfully are chosen to insure a smooth variation in the lx, ly tandem with opposite directionality in the affected shift of the scatter as one envisions it from the left to the right of the peak. The shift expression effectively pulls up one side of the scatter while simultaneously pushing down the scatter on the opposite side of the peak. This is consistent with the scatter compression evidenced in chapter 5 for point sources shifted away from the coordinate center. The parameters a2 and a4 are inversely related to the width of the scatter distribution in the lx ,ly direction respectively. The smaller the respective parameters 57 . Chapter 6: Implementation are, the faster the decay of the scatter distribution towards the edges of the planagram. If the scatter is concentrated in the vicinity immediately surrounding the peak, a rapid decay in the scatter profile will result in a larger parameter value. If the distribution conversely is broader, the scatter fall-off determines a commensurably smaller relative parameter value. The parameter a, simply designates the amplitude of the scatter in absolute terms relative to the peak magnitude (in counts). Prior to implementing the kernel within a convolution subtraction the parameter a1 will be replaced by an overall normalization factor related to the effective scatter fraction along each LOR. A discussion of this issue is postponed until section 8.1.1. In order to assess the best model power n, two comparison tests aire performed. A centralized point source, both radially and axially, will generate a symmetric scatter distribution about the projection peak. This permits the adding of the planagrams, thus augmenting the number of counts and so the statistical accuracy in the single summed planagram. The symmetric data will reveal the preferred scatter model in the sense that one model will more closely resemble the scatter profile than do the others. When the source is displaced off-center in the projection, the scatter distribution will become skewed. Whichever of the six possible functions best accommodates the position modified scatter, in the sense of optimizing the fitting routine, represents a candidate scatter model. The function that scores best under the combined regimes will be deemed the optimal function. 58 Chapter 6: Implementation 6 . 3 . A L G O R I T H M 6 . 3 . 1 . O P T I M I Z A T I O N C R I T E R I O N Each candidate function is fitted to the planagram scatter distribution using a maximum-likelihood based paradigm. 4 1 For ease of implementation an unconstrained weighted least squared merit function is used. The resulting chi-square statistic measures the cumulative error in modeling the observed scatter SlxJ, with the analytic model Smode,.41 , _ 2 'V V 1 ( \ . l i ~ ^ m o d e ; ( ^ ' ^ ' a i > a 2 » a 3 ' a 4 ' a 5 ) ^ 6.2 X =LX — However since only the region outside the projected peak contains scatter exclusively, the calculation of %2 must exclude the region coinciding with the projected point source as not to bias (with primary counts) the minimization. Furthermore, the region at the extreme ends of the scatter distribution will prove deleterious to the fitting procedure due to the statistical noise (see section 6.3.2). Both problems are remedied by using a. mask over select regions in the planagram. This has the affect of reducing the number of fitted points in exchange for improved algorithm performance. The selection df the mask is reviewed in the proceeding section (see section 6.3.2) The weights cte(y correspond to the stochastic counting error for each.'S,^ , , .Tn the least-squares formulation of maximum-likelihood these are implicitly gaussian.41 The %2 statistic of equation 6-2 therefore estimates the normalized (to unit variance) deviation of the model from the modeled data. The deviant or measured error oy, present in the planagram is not however in reality normally distributed. The positron 59 Chapter 6: Implementation ' • * decay and subsequent e + - e' annihilation occur randomly with a Poisson probability distribution. The experimental uncertainty in the planagram counts is thus Poisson distributed and equal to: Under large counting regimes the substitution of a Poisson deviate for a normally distributed deviate is validated based on the Central Limit Theorem which establishes the convergence of a discrete Poisson distributed statistic to a continuous Gaussian statistic in the limit of large counts. However the convergence is not uniform. Particularly in low count regimes (per bin) as is the case for the planagram scatter tails, the implicitly gaussian %2 minimization will underestimate the weight assigned to less likely events - those distributed near the tail of the scatter distribution. Nevertheless, the maximum likelihood method was adopted based on the availability of a routine optimized for dealing with non-linear model functions of the kind encountered here. Secondly, the minimization returns an easily contrastable relative measure of the goodness of fit i.e. the %2 . The value of the statistic will therefore be to gage the relative merit of each of the candidate functions at fitting the observed data rather than an absolute measure of a particular candidates equality to the scatter - the trial functions are based on heuristic observation. The selection of a specific minimization algorithm was strongly influenced by the non-linearity of all the candidate functions. The chosen. Levenberg-Marquardt method4 1 is optimized to converge to an absolute minimum in the neighborhood of local minima. For non-linear model functions, slight displacements along a parameteric trajectory may significantly displace the %2 since it is itself a non-linear 60 ' Chapter 6: Implementation functional of the parameters. The Levenberg-Marquadt Method avoids this difficulty by alternately switching between two different minimization regimes depending on the proximity to the minimum. Far from the minimum the algorithm assumes a gradient based steepest descent method. As the absolute minimum is approached, the gradient descent is replaced by a more precise inverse Hessian (approximate the minimum as a quadratic surface). However since the minimization is formulated on a Gaussian model of the error distribution, planagram bins containing zero counts pose particular problems. The observed count sum in a given bin is Poisson distributed and hence has a non-zero probability of registering zero. The gaussian model premised on a continuous count distribution does not assign any probability to a no count. Since the only estimate of the count deviation in a bin is given by 6.3 the lx,l contribution to the %2 is identically singular i.e. %] l ~ . Subsequently all planagram bins with zero counts are omitted.* The algorithm is implemented via a C coded program which accepts sinogram image files (truncated to the FOV). Following the normalization of the data to the crystal ring efficiencies, the program extracts the direct sinogram planes and reformats these extant planes into the heretofore detailed planagram projections. After suitably masking of the peak and neighboring region as well as the bins containing zero counts the merit function of Eq. 6-2 is enumerated. The first calculation of %2 required the ^Alternately , replacing the 0 count wi th a small positive number (less than one) wou ld similarly introduce a poorly specified contribution (the deviation of the count w i l l exceed the count).. 61 • Chapter 6: Implementation initialization of the parameter settings. Speedy convergence was discovered with, the following selection: 6.4 a, = .2*peak a2,a4=(.01)n a3,a5 = 0. Following the premier calculation of %2 the gradient of the merit function determines a measured adjustment to the model parameter settings. If the subsequent %2 enumeration shows an increase, the adjustments are reversed (adjust parameters in the opposing direction). The calculation of "i l s then repeated until it's change is no longer statistically significant i.e. it varies by less than a single degree of freedom. At this point a final Hessian type minimization of the %2 functional and last parameter determination is performed. 6 . 3 . 2 . S E L E C T I O N O F M A S K S The choice of masking affects the quality of the fit. Underestimating the extent of the peak region biases the fit to the primary counts. Including the outlying regions of the planagram poses the difficulty prejudicing the analytic fit by noise. The number of events registered for each planagram measurements range in the total number of counts recorded. For the current experiment this figure lies between 25,000 and 40,000 (with a mash of 3) in the scans performed here. Of the recorded counts, the bulk fall in a neighborhood surrounding the peak. The trailing scatter in the measured distribution tails registers fewer counts and therefore succumbs to large noise affects. The prominence of this affect is clear in the image texture generated by the typical scan. The jagged pattern is evident even in the 62 Chapter 6: Implementation ; • following blow up of the mashed central point source, where the counting statistics are relatively good as can be seen in Figure 6-9: Figure 6-9 Magnified View of Scatter From Centralized Point Source 100.-The selection of a mask was achieved in two ways. First ah optimal size for the central mask was determined based on a comparative %2 valuation with various windowing realizations. Secondly the qualitative success or failure of a fit was. tested with and without an.end mask. In formulaic terms, given a projected point source at lx',l' On a 128 by 31 planagram projection, the inner mask will extend for a region of 2lx X2ly i.e.: 6.5 j i _ i bound < Y < Y ' _|_jbound 5 X 1 . • ' X 1 X X . j [ " ' i , i i bound < ; i < [ ] • _ i _]b° u n d-y y ~ " A y ~ y h .63 Chapter 6: Implementation • Where•lbm,nd corresponds to the outer bin (as measured from the peak bin) delimiting the peak region. The outer mask is not well determined in the sense of the inner mask. It's exact size does not have an impact on the performance of the fit. A convenient limit was found to be ± 45 bins transaxially and ± llbins axially. -For the selection of the central mask (no outer mask) a gaussian fit was performed on the planagram of figure Figure 6-9. The delimiting bounds lxb°"" and l*""* were each allowed to vary independently from between 3 and 12, and 2 and 8 bins respectively. The resultant %2 (normalized to the number of degrees of freedom; the number of data points minus the number of parameter constraints) is tabled below: . Table 6-2 Selection of Inner Mask Ix 3 6 12 2 - 1.272 (1836) -* - 1.477 (1792) -5 3.130 (1834) 1.139 (1770) 1.485 6 - 1.058 (1748) -8 - . 1.038 (1704) -where the resident degrees of freedom (in parenthesis) in the planagram are given by the number of points used in the fit minus the number of parameters in the model. Note that the normalized %2 is not expected to approach one as this would only be expected in the case of a normally distributed random variable. According to the model above the optimal fit was achieved with larger bounds (8) along the axial coordinate. However, an axial bound greater than five while larger than the PSF in that diction was felt to simultaneously exclude too many points 64 Chapter 6: Implementation directly contiguous to the projection peak. This is important since the points near the peak do not accommodate any of the proposed analytic models. Specifically they are flatter than either an exponential or gaussian functional model. This means, that excluding these points would tend to exaggerate the success of the fit at the expense of overestimating the scatter directly underneath the peak. For this reason the central mask was selected (as a compromise) to be a square patch of 6lx5ly. To determine the necessity of an outer patch, the same guassian fit was applied to the previously elaborated central source projection. The resultant fits were however inconclusive in either %2 (1.139 (1770) with and 1.129 (3847) without) or appearance. This suggests two things: first that the presence of a mask is not critical for a centrally located source and secondly that the outlying points do not inveigh important information towards the fit (low count points imply large variances). This, is consistent with the earlier assumption that the exact bound for the outer bound is non-critical. When the same comparison is made for a radially offset source (in the projection), the above inconclusive result no longer holds. As an examination of the following diagram reveals, the fit is wildly altered by extricating the outlying points from the planagram: 65 Chapter 6: Implementation Figure 6-10 Cross-section of Radially Offset Source, With (dotted) and Without (dashed) End Mask where the source position in the planagram is given by the r and z coordinates underneath the plot. Although the %2 changes from only 1.06 to 1.17 between the "good" and "bad" fit (reflecting the minimal merit of %2 as an absolute indicator of the goodness of fit) the qualitative difference is overwhelming. The scatter fit exhibited with masking more clearly resembles the intuited form of a scatter profile as a distribution peaking underneath the projected activity source peak. Therefore from herein the mask is adopted in all further fittings. 66 Chapter 7: Execution 7. E X E C U T I O N 7.1. S E L E C T I N G T H E M O D E L The six trial functions expounded in equation 6.1 were contrasted on a central and radially offset source planagram projection both previously elaborated in section 6.1.1. The efficacy of any single model is explored on two counts: first in its ability to account for the salient features of the scatter distribution as the source is displaced in the projection, this is reflected in a minimal %2. Secondly the function representing the scatter should vary smoothly (i.e. continuous first derivative) across the planagram; consistent with the measurement process. The six scatter planagram distributions generated by fitting the trial functions to the data from the scan of a central source below (normalized to the total number of counts in the planagram) is seen in Figure 7-1. 67 Chapter 7: Execution Figure 7-1 Scatter Planagram Distributions From Central Source Clearly the magnitude of the scatter (all the graphs are scaled to the peak height of n=l) diminishes with increasing n as one proceeds from the upper left n=l (exponential) distribution to the lower right n=6 (gaussian). This effect is concurrent with a decrease in the width of the distribution for the lower n functions as reflected in the increased kurtosity (cuspiness). The sharply peaked non-analytic scatter distributions for lower n's run contrary to the expected smoothness in a measured scatter distribution. Moreover, the lower n functions fall prey to overestimating the scatter underneath the projection peak. This is based on the previous insight of section 6.3.2 wherein the increase in the scatter as the projection peak is approached was noted to taper off. This affect is made more obvious by examining the transverse cross-section through the projection planagram with the six fits superimposed for clarity: 68 Chapter 7: Execution Figure 7-2 Transverse Cross-Section - 2 0 - 1 0 0 . 1 0 . 20 r (cm) r= . 0.468750cm. z= 0.00000cm The former insight is confirmed by a comparison of the %2 valuation for all six Table?-]/ %2 Comparison For Fit on Central Source n-2 n=3 rt=5 n-6 (1.0) (1.2) (1.6) (1.8) (2,0) 1.256 1.184 1.144 1.127 1.127 1.139 models.* Note the value of the matching exponent is listed in parenthesis underneath. * A full summary of the each fit and derived parameter settings is given in appendix 1 + In order to gleen significant information from the %2 value, a precision of 3 decimal places was kept. Since the degrees of freedom is on the order of IO3 and a statistically significant x 2 difference is consititued by a variation of 1, a precision to 1/10 3 was maintained. 69 Chapter 7: Execution According to this last comparison,, the optimal scatter model lies in the neighborhood of n-A:n-5. When the same %2 consideration is applied to the offset source the result is: Table 7-2 Y 2 Comparison For Fit on Offset Source • • i i i if iViii i i l l MiM l i i l A l i l l r .4 . i i i i i i i i i i i i i i i i i i i i i i i i i i i l l l i > iu.rtii i 1111111 iin.i . i . . M.II IIIIIIIIIIIM.IIIIUIIIIIIH.I,,,, m iv=2 n-3 n=4 m**5 n=6 <1,0) (1,2) <1,4) (2.0) • ti 1.130 1.138 1.145 1.153 1.16.0 1.167 In this latter case the lower n functions were more successful at fitting the scatter tails. How accurate an indication of goodness of fit the %2 statistic is, is not clear since the data points adjacent to the peak do not contribute to it. Regardless, the maximum difference in the %2 value between functions is minor (2% for the offset source and 10% for the central source). The reason stems from their roughly equal success at fitting the central portion of the scatter tail (excluding the beginning portion at the peak and end portion at the scatter base). This is evident from a comparison of a cross-section through the offset source planagram: * See appendix 2 for detailed information. 70 Chapter 7: Execution Figure 7-3 Transverse Cross-Section 50 - i — i — I — r - 2 0 ~I 1 1 1 I 1 - ~ i — i — i — r ~ . / V W \ . A , - 1 0 10 . • 0 r (cm) -7 .34375cm. z= 0.00000cm where the model functions range from n=l (exponential) to n=6 (gaussian) as one 20 r= proceeds down through the six fits. . The matching planagram projections are displayed in Figure 7-4: 71 Chapter7:-Execution . Figure 7-4 Scatter Planagram Distributions From Offset Source Although clear differences emerge, in the qualitative feel of the disparate scatter distributions, the ^determinant does not absolutely disqualify one function from another. The preference for a smooth scatter distribution suggest the use of a model with higher n. According to this paradigm the gaussian with n=2 will herein be adopted as the model function. 7.2. C A L C U L A T I N G H Having selected the scatter model the calculation of the scatter parameters on the planagram domain can be derived permitting the implementation of a position sensitive convolution scatter correction using either the kernels horh'. The kernel h is derived from the planagram projections generated at each angular view. By fitting the 2D scatter distribution, each planagram effectively delivers a set of fitted measurements for the scatter model parameters, indexed to the projected source Chapter 7: Execution position in that planagram. The set of these measurements over all available angular views represent a sampling of the parameter values on the whole planagram. By further hypothesizing a characteristic parameter distribution for each parameter, the samples may be used as input to an interpolative fit. The corresponding parameter distributions (which are now explicitly position variant over the planagram domain) form the bases for the position dependent convolution correction. The proposed strategy is alluded to in Figure 7-11. Following the application of the gaussian fit to each planagram projection, the derived parameter settings are indexed to the radial i.e. transverse lx and axial ly planagram coordinates. Each axially separated source position: z=0.00cm, z=1.69cm, z=3.71cm, corresponds to a matching Zy position in the planagram. In the diagrams that follow, each scan defines a series of parameter samples at several values of lx. Each sample originates in a distinct planagrams. Again for the benefit of improved statistics six planagrams were added together resulting in thirty-two distinct planagrams for each scan. The result of the parameter sampling is summarized below: 7.2.1. A , R E S U L T S The parameter at determines the scale of the scatter distribution. Each of the three plots in Figure 7-5 below corresponds to a different axial source position (relative to the axial center z=0.00cm) within the object space. Each of the a t samples is taken at sequential angular views (within each view the source projects to an alternate lx, denoted "r" below, point inside its planagram). 73 Chapter 7: Execution Figure 7-5 a1 Distribution a Variation i ; with Transverse Source Position in the Planagram projected source position Since the cylindrical phantom supports a reflective axial symmetry about the central axial plane the three distributions correspond to an effective sampling of five axial source positions (translate the outer two sampled distributions through the central z=0.00cm direct plane). The relatively higher magnitude of the parameter a1 at larger transverse (r) offsets and at larger axial offsets (z) stems from the underlying magnitude of the measured counts at these positions as is witnessed in Figure 7-6 below: 74 Chapter 7: Execution Figure 7-6 Comparison of Measured Counts as a Function of Radial and Axial Projection Position 6000 4000 2000 0 6000 4000 2000 0 6000 4000 2000 0 0 r (cm) 0 ' (cm) • (cm) The root of this difference lies in differential amount of attenuation encountered by the emitted radiation when the source positioned at the phantom edge projects to a central projection bin and alternately to an edge projection bin. 7.2.2. A 2 R E S U L T S The coefficient a2 is proportional to the inverse of the width of the scatter distribution in the transverse lx (r) direction. The corresponding determination of its value as a function of r and z is given in Figure 7-7: 75 Chapter 7: Execution Figure 7-7 a 2 Distribution a Variation 2 with Transverse Source Position in the Planagram •0.003 . . . . . . . . . . . . . . . . . . . . . . r . . As expected the value is everywhere negative (recall that in the formulation of the scatter model the sign is incorporated into the parameter). The second notable feature of the disparate axial distributions is that they seemingly coincide. This means that the parameter is almost exclusively determined by the radial source position in the planagram. The Value is seemingly indifferent to any axial variation of the projected source position. 7.2.3. A 3 R E S U L T S The parameter a3, one recalls was originally invoked in order to compensate the model scatter for the asymmetry observed in the scatter data when the projected source position moves away from center. Since the parameter's multiplier is the relative distance from the projection peak, one expects the parameter to change 76 Chapter 7: Execution ' ' sign as it passes through the radial center of the planagram. This is indeed observed in Figure 7-8: Figure 7-8 a3 Distribution a Variation 0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 with Transverse Source Position in the Pianagram z=0.00cm lm ~' z=1.69cm A - . . . z=3.71cm >* • 1 ' • km* • • / - ' X '. 1 1 1 i i i i i i —i i i t 1 1—1 1 1 L _ r i i i r t- , (cm) projeoled source \uYntt<iti The reflective symmetry is a consequence of the inherent tomographic symmetry - a reflection about r=0 is equal to viewing the source from an angle displaced by 180°. Again, the axial variation in the parameter measurements seems minimal. 7.2.4. A4 R E S U L T S The radial variation of the parameter a4 (which determines the axial width of the scatter distribution on each planagram) is displayed in Figure 7-9: 77 Chapter 7: Execution Figure 7-9 a4 Distribution a Variation • . • 4 with Transverse Source Position in the Pianagi 0.02 I 1 1 1 1 1 r^- i 1 r — I • 1 1 1 1 1—i 1 1 — | 1 1— Only the z=0.00cm line can be said to exhibit a concave shape, with positive variance. The further out axially the source is moved the less behaved is the distribution -vacillating wildly about zero. Constraining the parameter to a value below zero has the affect of truncating the above graph to only the lower plane (upper ceiling at 0). It was decided to allow for a complete variation in the parameter value at this stage. The restriction of the parameter to a value below zero will only become necessary (if at all) during the proceeding interpolative fitting stage (see section 7.3.5). 7.2.5. A . R E S U L T S The final model parameter a5 has a parallel function to that of a3 except along the axial ly (z) coordinate. Since minimal axial scatter asymmetry is observed, it 78 Chapter 7: Execution may be expected that the parameter remains near zero for all points in the planagram. Figure 7-10 below reveals this to be the case: Figure 7-10 a5 Distribution a Variation The parameter straddles the zero mark. Unlike a3 which expressed an odd symmetry in r, a s is even in r (although odd in z by definition). The parameter is taken to be quadratic in r (its variation is taken to be cubic in z; see section 73.1). 7.3. G E N E R A T I N G T H E C O M P L E T E K E R N E L In order to characterize the scatter response across the whole of the planagram domain the defining parameter values must be interpolated (plus the scatter fraction thus delineated) over the unsampled planagram domain. This is made possible by constraining a model for the parameter distributions to pass through the sampled parameter values. Since the parameters are only defined on the region supporting a 79 Chapter 7: Execution scattering material (delimited by the phantom boundaries in the projection) the modeled parameter values are restricted to a transverse 20cm diameter. The program thus envisaged for constructing the kernel h requires the assembly of a position variable kernel by way of referencing to position dependent parameter distributions via a look-up table implementation i.e. where the look-up table assembly is summarized in the following diagram: Figure 7-11 Look-Up Table Implementation Of Assigning h : — : _ • — lx The kernel assigned to the point lx',l' is generated from the subordinate parameters a ^ assigned to the point / / , / ' . The process of indexing the kernel is thus reduced to the indexing of its argument parameters. The proceeding section therefore accomplishes the task of generating the parameter distributions across the (reduced 20cm diameter) planagram domain. 80 Chapter 7: Execution 7.3.1. F O R M I N G T H E P A R A M E T E R D I S T R I B U T I O N S Based on the sampled parameter values of section 7.2, the complete parameter spectrum is determinable for each of. aj-a5. This is achieved by fitting a separate polynomial function through each of the sampled parameter spectra. The functional model describing the parameters is assumed separable in the IJ. directions: Based on the symmetry considerations detailed in section 7.2, the models are either even or odd along lx,ly.about the planagram center - herein taken to be bin Z/=64,-/y'=16. The algorithm used to perform the fit is the same already elaborated in section 6.3. In the axial'/ direction, the.inherent tomograph symmetry imposed the constraint that the leading coefficients for all the odd polynomial powers were identically zero except for the ease of a5. Moreover, powers above two were too sensitive to fluctuations in the parameter samples and tend to admit multiple inflections (multiple miii-max regions) in the axial fit which is not born out by physical argument. By construction the axial reflective symmetry of a5 is odd not even (recall, the coefficient multiplier changes sign upon translation through z=0.00cm). This means that the functionality of the axial variation must be similarly odd. This fact predisposes the function to a polynomial up to the cubic power in this case only. Along the transverse planagram direction, separate treatments were required for the disparate parameters in order to achieve successful fits. Inspection of the sampled parameter values suggests an even function for all the parameters except a3 which is clearly odd in symmetry. Due to the possibilities of introducing unphysical multiple inflections for all the even parameters (with the exception of al), only quadratic 81 Chapter 7: Execution transverse models were used. The subsequent functional forms used in modeling the parameter distributions are summarized below along with the commensurate distributions. 7.3.2. A, The proposed model for the projection distribution of a1 is as follows: 7.2 at(lx,ly) = {bx+b2*(lx-lx')2+b3*(l x-lx')>{b4+b5*(ly-ly')2} with respect to the planagram bin number and where lx',l' are (64,16) as mentioned earlier. In order that ax be translated to a physically meaningful quantity, the parameters are scaled by either (.3125cm/bin)n or (.3375cm/bin)n depending on whether the parameter multiplies lx or I and where h is equal to the exponent of the coefficient multiplier. b1 and b 1 2 are not scaled. For all the ensuing model specifications the former scaling will be implicit. Based on the above model the corresponding distribution appears as in Figure 7-12: 82 Chapter 7: Execution Figure 7-12 Planagram Distribution For a t 6.5 6.0 5.5 7.3.3. A2 The model function used to describe parameter a2 is: 7.3 a2(/x,/y) = ^ ! + fo2*(/,-4')2}*{^ + &5*(/,-/;)2} The parameter a2 refers to the inverse of the width of the scatter distribution along the transverse direction through the planagram (and inverse square of the deviation). As the peak is displaced axially, the parameter decreases, meaning that the width increases, all be it slightly. When the peak is displaced in the radial coordinate the reverse occurs. The joint effects are explainable if one considers what happens when the source is displaced first radially for a given axial position and then axially for a given radial position. When the peak is shifted radially, the scatter becomes compressed and 83 Chapter 7: Execution more localized in one region of planagram. This same effect is not evidenced from the transverse perspective when the shift is in z. In this latter case the overall scatter fraction decreases relative to that found at z=0cm. This effect implies a flattening of the scatter distribution (less relative counts at the peak position). This flattening or broadening will appear as an increased width. This same affect should be evidenced for the width in the axial distribution of the scatter response as reflected by a4 although to a less pronounced extent (since the perspective is reduced to the 31 versus 128 sampling bins; see section 7.3.5 for more details). The subsequent distribution is given in Figure 7-13: Figure 7-13 Planagram Distribution For a 2 84 Chapter 7: Execution 73A. Recalling the earlier discussion of a3, its.model description deviates from the previous two. ' ' 7.4 a3(lx,ly) = {bl+b2Hlx-lx') +b3Hlx-lx')3}*{b4+b5*(ly-ly')2} netting the-fit of Figure 7-14: Figure 7-14 Planagram Distribution For a3 As expected the function is odd in r and relatively flat. Its value is zero in around r=0cm conforming our expectation that in this domain the scatter kernel should be radially symmetric implying the coefficient of the asymmetric term should fall to zero. .85 Chapter 7: Execution . 7.3.5. A4 Second order polynomials were found to describe the data adequately. The model reads: 7-5 . aA(ls,ly) = {k+h^d^-l^^+h, *(ly-i;f] • with the matching distribution of Figure 7-15: . ; Figure 7-15 Planagram Distribution For a4 o.oo° r -0.002 p-_o.oot U -0.006 -0.008 U -0.01 °J ; 7.3.6. A. Lastly a5 affords special Care due to its poor parameter estimation and odd axial symmetry. The corresponding model;is given by: 7.6 a5(ix,iy) = {b{ +b2*(ix-KfMh +.K*(iy-i;)} with the matching fit given by Figure 7-16: 86 Chapter 7: Execution Figure 7-16 Planagram Distribution for a 5 . 0 . 0 0 1 0 i.oooo _ o . o o i ° P -0.0020 p -0.0030 p _o.oo*°J 7.4. I N T E R P R E T A T I O N The sundry parameter distributions do not conform to physically meaningful values everywhere on the planagram domain. This is clear since by definition they have meaning only on the region supported by a scattering material. Based on the scatter model elicited from each sample planagram, a scatter fraction was determined according to equation 3.10. Like the heretofore expounded parameters, each angular viewof the point source effectively contributed an SF sample at a point on the planagram domain indexed by the projection peak, position lx',l' of the 87 Chapter 7: Execution underlyingly sampled planagram. Fitting the resultant.samples of the SF with equation 7.2 gives the following: Figure 7-17 Scatter Fraction Distribution on the Planagram Domain 0.35 0.30 0.20Y-0.10k 0.05J. Outside the scattering object, no material exists to assume the perfunctory role of scatterer. Extrapolating into this region the sampled SF values or in equal measure the sampled parameter values constitutes a violation of the maxim conjecturing h as a material response to an embedded point source. If one wants to define the scatter kernel parameters in a region larger than 20cm in diameter, then a wider phantom must be used. This fact will assist in the execution of the correction since it reduces the physical space over which the numerically intensive convolution takes place. 88 Chapter 8:Test Of Correction , - •.' 8. T E S T O F C O R R E C T I O N 8.1. A L G O R I T H M The method of scatter, correction using h has already been elaborated: the planagrams are extracted from the sinograms, and then iteratively convolved with the kernel h indexed through the parameter arguments to the projected source coordinates •//,/'. The convolution regiment requires the preassembly of each lx';l' specific kernel array. For a convolution domain* given by -10cm</v'<10cm (64 transverse projection bins) -5.3cm</y'<5.3cm (31 axial projection bins) this requires 64x31 separate kernels. Each kernel extends 127 transverse by 61 axial bins. The large size insures that the kernel is able to cover the whole of the convolution domain even when centered on the first element of the convolution domain i.e. lx'=bin 1, /y'=bin 1. The subsequent convolution reduces to a program of scaling every kernel.element in a kernel array by the image pixel to which the kernel is indexed. This process is repeated for each of. the 64x31 convolution domain elements and their representative kernels. The result of each scaling is summed to give an overall scatter characteristic to the" input image planagram. Due to the size of the kernel and nature of the convolution process, the For the purposes of a clinical implimentatipn, the convolution domain is chosen to. coincide the projected portion"pf the subject being investigated. Chapter• 8-Test Of Correction scatter array is 190+ transverse by 61 axial bins large. However, only the central 64x31 of this last summed series correspond to the actual scatter for the input image planagram. The full algorithm is diagramatically represented in Figure 8-1. The 190 is the sum of 64 direct domain elements plus 63 padded transverse elements to either side of the regular domain. • • . • 90 ;'; Chapter 8:Test Of Correction Figure 8-1 Convolution Algorithm Fox? Planagr etm. 1^ 31. C o n v o l n t l D n : , - " " Each P ixe l Element i ' j in P]anagram Multipl ies Elements o f Kernel i , j For the purpose of final implementation, this algorithm is optimized by prudently selecting out for scaling only the 64x31 kernel sub-array that directly contributes to the Chapter 8:Test Of Correction ultimate 64x31 scatter array. This improvement obviates.the need to construct the full 127x61 kernel as or the 190x61 holding matrix for the output. The scatter correction is performed iteratively. The measured image planagrams are deemed an initial estimate of true primary count distribution (see section Since the kernel h is defined complimentary to the primary distribution P, the subsequent convolution on the whole of M will return an exaggerated scatter contribution: 8.1 S\lx,ly) = P\lx,ly)®h{lx,ly); {P0(lx,ly) = M(lxJy)}>P(lx,ly) Consequently the second iterative estimate of the primary distribution P 3 =M-S° is underestimated. Because of the underestimation of P in the second iteration, equation 8.1 yields an underestimated value of S, S1 . The subsequent estimate of P given by P 2=M-S 3 will therefore be greater than the correct value i.e. P2>P: This cycle repeats, oscillating around the true primary count value with decreasing amplitude after each iteration (see Figure 8-6). With each iteration the convolution corrected estimate of P converges to the optimal (in the sense of a %2 test) primary count value. The stopping criterion for the correction is designated as the iteration which . leaves the primary count distribution estimate statistically unchanged from the previous estimate. This measure is calculated using the familiar %2 statistic calculated from the difference between the original measured distribution and the most, recent primary count estimate (normalized to the original variance in each data point). When 92 Chapter 8:Test Of Correction consecutive enumeration's of the %2 statistic return values whose difference A%2 is less than 1, the primary count estimate remains effectively unchanged between iterations. This implies that the convolution correction may be terminated. For eventual clinical application however, the series of iterative calculations in the algorithm may be terminated early. Little improvement in the correction is gained by proceeding beyond this point. This fact is demonstrated in Figure 8-6 where the amplitude in the oscillations about the optimal %2 dampens to negligible amounts following the fourth iteration. 8.1.1. N O R M A L I Z A T I O N The normalization of the kernel is a critical component of the correction. Equation 5.4, establishes the normalization criterion: 8.2 N = S F J/iCfl^ ,03,a4,as.) projection whereby N corresponds to a multiplicative normalization factor calculated at each point lx',ly' in the planagram on which the convolution is centered. Since the implementation of the convolution occurs on the reduced 64 by 31 bin planagram region, only the kernel integrand overlapping this integration domain contributes to the new normalization factor N. -However, applying the new normalization factor using the scatter fraction SF distribution calculated from the offset source (see Figure 8-2 below) will underestimate the overall normalization. 93 Chapter 8:Test Of Correction Figure 8-2 SF Measurements From The Offset Source SF Variation wifh Transverse Source Position in the Planagram s The reason lies in the way the SF was calculated across the planagram domain. Consider a cross-section through a cylindrical phantom of 20cm diameter supporting a point source offset along the central chord from the cylinder center by 7.34cm. The resultant planagram projection at each view is a superposition of the primary activity and scatter due to the source emitted y-ray's interaction with the surrounding scattering material. Both the scatter and primary distributions will be peaked at the planagram coordinate corresponding to the projected source position. Assuming a method exists to extricate the scatter and primary counts from the total measured count, then a scatter fraction SF value in the sense of equation 3.10 may be assigned to the planagram bin underneath the projected source position. 94 Chapter 8:Test Of Correction When this same source is viewed at a different angle inside the tomograph the projected source position shifts accordingly, The matching scatter and primary counts will also adjust due to the different material distribution now intersecting the source and detector array relative to the previous view. The result is a new SF determination corresponding to a new projected source position. Each succeeding view therefore generates a sampling of the SF at disparate points in the planagram projection. This sampling process may equivalently be conceived through the rotation of the source supporting phantom through angular increments equal to the angular variation between adjacent views. When viewed as a transverse slice through the planagram projection (/y=constant), the change in perspective corresponds to a transformation between a source stationary frame-of-reference to a planagram stationary frame-of-reference delineated by the rectangular coordinate system (lx,r). The lx coordinate references the usual transverse planagram coordinate while the r coordinate refers to a L O R ray running perpendicular to it. The process is summarized in Figure 8-3 with the SF sampling result (for each of three scans of section 5.3.2) given in Figure .8-2. 95 Chapter 8:Test Of Correction Figure 8-3 SF Sampling Procedure The Sampling Procedure Viewed From the Perspective of a Point Source Positioned inside a Phantom at (lx,r): 10cm —i LOR2 (r) o -10cm ((7.34cm)A 2-(lx')A 2) "1/2 Point Source -10cm 10cm Corresponding Scatter and Primary Distributions on the resulting Projections, Determines an SF Evaluation For Each source Position (lx,r): SF(lx=0cm;r =7.34cm ) -10cm 0cm 10cm SF(lx=8cm;r =((7.34cm)A 2 -(lx')A 2) «i/2 ) Primary Scatter -10cm ." • 0cm 10cm If one is to occasion a second such SF measurement but now based on a source displaced along a LOR from the original offset position, the corresponding distribution is changed in scale from the original. On first approximation this effect along with other differences between disparate source based scatter profiles may be overlooked. This was in fact done in section However, the scatter correction is for the purposes of higher precision sensitive to variations in scatter fraction (as opposed to the scatter shape which is manifest in the form of the kernel h). The implication of this sensitivity 96 ' . Chapter 8:Test Of Correction is that the ultimate scatter fraction SF used in the correction must take into account the fact that different points (lx,r) inside an arbitrary distributed activity source' will A n arbitrary source distribution may be resolved conceptually into a superposition of spatially distinct (although contiguous) radiating point sources. The corresponding SF at any point in the planagram may therefore be expressed as the superposition of the SF contribution from the resolved element point sources along the L O R specific to each planagram bin. Put another way, the effective SF that is assigned to each projection bin must be considered an ensemble average of SF's generated by a hypothetical string of activity sources stretched across the whole LOR. The normalization factor N must account for this superposition. This is achieved through averaging of the SF contributions along the whole extent of the LOR. The SF distribution of Figure 7-17 only corresponds to SF samples taken at a single point in each LOR: r= J\l34cm)2 - (lx f for -7.34cm</i<7.34cm -5.3cm</y<5.3cm. generate different SF values. A separate determination of an effective SF is required for every point in the projection over which the kernel is defined, namely -10cm<^<10cm, The methodology j of calculating an SF for a given'//,/' in the above domain is developed, for the sake of illustration, on a ID projection consisting of a single lx coordinate. The source position within the phantom is designated" by the rectangular coordinate pair (lx,r), where - . redg<r<redge where redg=^(l0cm}2 -(lxf for the particular case of a 20cm diameter (10cm radii) phantom. Hence each source placement at a position {l=lx',r=r') determines an 97 Chapter 8:Test Of Correction independent SF=SF(/x',r')- Herein, the scatter fraction is regarded as a function of lx,r, i.e. SF=SF(Zx,r). The determination of an effective scatter fraction at. a projection position l=V reduces to that of calculating (SF(/x=/x',r))r, where the <»)r operator corresponds to an analytical averaging of the argument over the range of r. The averaging operation must be done with a limited sampling of the scatter fraction given by the combined values in Figure 8-2 (r=7.34cm) and SF(/x=0cm,r=0cm) measured when the source was placed at the center of the cylindrical phantom (see Figure 3-4). It is therefore necessary to elaborate a functional form for SF(/x,r). The SF(/x,r) variation may at first be assumed to separate out into a variation along lx and r: SF(/X,r) =L(/x)R(r). L(ZJ is already given by the scatter fraction distribution of Figure 7-17. R(r) must be deduced. The first inference on the form of R(r) is based on the coincident nature of acquisition in PET (and the subsequent reflective symmetry about r=0cm) R(r)=R( I r\). Secondly, the nature of R's variation with r is inferred from the fact that the SF is highest at r=0 and lowest at the phantom edge r=10cm. In between it diminishes quickly near r=0cm, with the SF fall off down towards its low value as the source is displaced out to r=10cm (at Zx=0cm).46 This last fact suggest an R bearing an exponential decay fromr=0cm to r=rcis=^(\0cm)2 ~(lxf .* A n improved evaluation of R is outside the scope of the present research. Any future research into the form of R may amend the calculation of (SF(/x=Zx',r))r/ accordingly. * Based on a cylindrical phantom wi th a circular profile at /y=constant 98 Chapter 8:Test Of Correction Thus R( I r I) is taken to have an exponential form given by 8.3 R(\r\) = ae~m Since SF(/x,r) takes its maximum value at r=0cm and its minimum at the phantom edge r-redge the boundary conditions determine a-R(0) (which may be expressed p m , SF(lx,0), . . An(SF(lx,Ocm))-ln(SF(lx,red.e)) • ;; as R(0) ) and b = j—-j - — . These boundary conditions fully characterize R( I r I) at a given lx, namely: 8.4 R{\r\) = LEMle-( W ™ L(4) The average scatter fraction measured for a projected source position at l=lx' (SF(Zx=/x',r))r, is therefore the ray length (2x\rat \) normalized integral of the S¥(lx,r) distribution along r. r«dge redge \drSF(lx=lx',r) L(lx=lx') \drR{r) 8.5 {SF(lx =lx', r)r = ~*- ^ = — edge [df jdr Substituting equation 8.4 into equation 8.5 and evaluating gives 8.6 (SF(lx = K ,r)r SF(lx=lt\0) C • l n ( i f ( I l = I x ' , 0 ) ) - l n ( J F ( I t = ( t > ^ „ ) ) I ' l \ ln(SF(lx =/x • ,0))-ln(SF(lx =lx' ,r^,)) Equation 8.6 prescribes how to evaluate the effective scatter fraction for each 1-e ) point lx on the projection irregardless of L(Zx). However, it does assumes a knowledge of the projection distributions SF(/x/0) and SF(/ x ,r % ) which according, to the earlier ln(*(0)) - HR(redge)) = ln(*(0)) - HR(redge)) + ln(L(/x)) - ln(L(/x)) = l n ( ^ - ) - l n ( ^ ^ ) L(L) L(lx) 99 Chapter 8:Test Of Correction formalism on the variable separation of SF(/x,r) should just be L(QR(0) and L(/JR( I redge I) respectively. Describing SF(/x,0) and.SF(Zx/r ) though proves more convenient in terms of the sampled scatter fraction values directly (non-analytically). SF(lx,redge) is taken: from the samples for the three scans at ^(lOcm)2 + (lxf =7.34cm (Figure 8-2). SF(Zx,0) may be recovered from SF(lx,redge) based on two observations. First, the scatter fraction at /i=0cm,r=0cm is .08 higher than it is at /I=0cm,r=7.34cm (previously measured). This suggests that SF(lx,0)=SF(lx,re((ge) +.08. However, this last ansatz must be corrected for the following: when the source at r=0cm is displaced along the lx coordinate to the phantom boundary at I lx\ =10cm the resulting scatter fraction SF(/;c=10cm,0) is identical to SF(Zx=10cm,r ) with redg=0 (as IlxI =>10cm, redge=>0cm). This fact mitigates against a direct adoption of SF(/x,0)=SF(Z;c,r e^) +.08. Instead a variation of this last form is adopted which accounts for the former extenuating consideration: OR 8.7 SF(lx,0) = SF(lx,reJgJ + -^-H\Ocm-\lx\) 10cm 08 where the correction factor '.——-*(l0cm -\lx\) insures SF(/r,r)'s compliance with its 10cm expected value at lx=0cm and I lx I =10cm. The origin of equation 8.7 is wholly heuristic. Future research efforts into this area may improve upon this result. 8.2. A L G O R I T H M D E M O N S T R A T I O N In order to demonstrate the validity of an iterative convolution subtraction using a position dependent kernel the scatter correction was performed on three test scans. The first of the three scans was a. 3D scan performed on a phantom consisting of a 100 Chapter 8-.Test OfCorrection central 2.3cm diameter activity solution filled cylinder centered within a bath of lower "background" activity enveloped inside an elliptical shell. The Plexiglas shell had a major axis of 19cm, a minor axis of 14.5cm and an axial length of 14cm. The contrast in activity between the two chambers was measured to be 3.8. The second scan was a 2D scan of the same two chamber source. The 2D scan and subsequent scatter correction is used as the "gold standard". The third scan, a variation on the. former idea of contrasting, activities, was performed on a three chamber phantom. In this last case two 2.3cm diameter cylinders, one of high activity the other of no activity were diametrically offset within a bath of intermediate background activity (high-to-background ratio of 4.0). The reconstruction of a transverse slice's activity distribution (non-scatter corrected) in Figure 8-4 demonstrates the geometry of each phantom. 101 Chapter 8:Test Of Correction Figure 8-4 Test Phantoms 2 Chamber Phantom 3 Chamber Phantom For the sake of clar i ty the greyscale has been reversed so that the range of whi te to black corresponds to a scale of l o w to h i g h act ivi ty. A l s o for ease of reference the two source objects w i l l here in be denoted as the hot spot source (2 chamber phantom) and the hot and c o l d spot source (3 chamber phantom). The appellat ions hot and cold a l l u d i n g to the presence of higher and lower activity chambers relative to the " b a c k g r o u n d " act ivi ty of the outer chamber. A n iterative correction of the type ou t l ined i n section was executed on each of two scans. The test is per formed us ing the already de te rmined scatter kernel h i ndexed b y the pos i t ion var iable parameters a2(lx',ly), a 3 ( / / / / y ' ) / a 4(Z/,/ y '), as{lx',ly) and norma l i za t ion factor N(lx',l') (adjusted for the e l l ip t ica l construct ion of the scanned objects). 102 Chapter 8:Test Of Correction For the purpose of illustration the procedure is detailed in the case of the central hot spot scan. Although in the final implementation the scan is neither mashed nor normalized to plane efficiencies, for the sole purpose of efficacy in demonstration, here the planagrams are plane-wise normalized and mashed by a factor of 12. The result is a series of 32 angular projection planagrams of the kind in Figure 8-5. Figure 8-5 Planagram View of Scanned Central Hot Spot Object 4 0 0 ° r I 2000 k iooop-The central plateau extending along the axial coordinate concedes the higher activity in this region. The scatter component of the projection is clearly visible outside the projected source in the transverse direction. A n iterative convolution subtraction of the kind described in the last chapter is then applied to this projection data on a planagram by planagram basis. The result is a series of improved estimates of the primary count 103 Chapter 8:Test Of Correction , ,' - - • distribution. This improvement is witnessed in the following plot of the %2 and change in x2:A%2 between iterations: '." .. • Figure 8-6 Variation in %2 and A%2 with Successive Iterations * — - 2 X — E 3 — 'rj "H i—* i—t 1 — £ 1—f J—( 1—E 1 — ' Number-of Iterations. The result after 17 iterations on the first of 32 mashed planograms is summarized : in the following pictorial (where the upper left projection corresponds to that of the raw measured data i.e. Figure 8-5) ' , 104 Chapter 8:Test Of Correction Figure 8-7 Evolution in Scatter and Primary Count Planagram Distributions P C 1 6 > 1 - 5 0 0 ooo 5 O 0 O J . Following 17 iteration the difference in X? (between successive iterations) falls below 1. Depending on the'view, the number of iterations necessary for convergence ranged between 14 and 18 on the remaining 31 planagrams. The seeming success of the correction is evidenced by considering, the transverse and axial perspectives. (summed along the respective conjugate coordinate) of the original and scatter corrected planagram. For the transverse plane in Figure 8-8, where the dotted line corresponds to the corrected planagram value, it is clear that the correction is capable of eliminating the obvious scatter fail outside the activity source: , 105 Chapter 8:Test Of Correction Figure 8-8 Transverse View Through Planagram: Original (solid upper), Corrected (dotted) and Scatter (solid lower) • - 2 0 ' , - 1 0 0 10 ?0 • • . r (cm) . Similarly the axial view through the planagram appears in Figure 8-9:' 106 Chapter 8:Test Of Correction •. Figure 8-9 Axial View Through Planagram: Original (solid upper), Corrected (dotted), Scatter (solid lower) 2.0x10" 1.5x10" cn •£ <= 5 1.0x10" o o 5.0x10^r-- 6 1 r r-_ J i L _ 0 z (cm) The net affect of the correction on the first sample planagram is to generate a scatter fraction (measured within the 20cm diameter central region) of 32%. 8.3. C O M P A R I S O N OF T H E C O R R E C T I O N E FF ICACY For the final implementation the former algorithm (in its optimized form) is applied to both scans. O n this occasion the acquired sinograms are not modified prior to correction. Only following correction is a normalization and attenuation correction applied. For the purpose of comparison, a position invariant convolution correction detailed in Bailey et al. 3 6 was simultaneously evaluated. Both corrections are then juxtaposed to a scatter uncorrected reconstruction. In addition the central hot spot 2D 107 Chapter 8:Test Of.Correction scan is reconstructed in 2D mode using a linear scatter background subtraction technique." The form of the comparison is to calculate the hot to background activity ratios in the reconstructed images. Only the 31 direct planes are used in this comparison. The activity level in each chamber is determined by delineating a region of interest (ROI) within each imaged chamber. The effective generating activity (inside the source) is thus determined as the averaged pixelated intensity inside each image ROI. The construction of the ROI is demonstrated in Figure 8-10 for each of the two scans. A scatter estimate based on the activity from outside the F O V is subtracted from the measured activity inside the F O V . 108 Chapter 8:Test Of Correction Figure 8-10 ROI Positioning Inside the Test Phantoms Hot Spot Source Hot and Cold Spot Source Based o n the R O I of F igure 8-10 three ratios are constructed. The first, part icular to the hot spot source, is the ratio of act ivi ty i n the hot spot to that i n the background . The second is ident ica l to the first, but evaluated for the hot and co ld spot source. The th i rd is the act ivi ty ratio of the co ld spot to the backg round i n the hot and co ld spot source. The result of the calcula t ion for the hot spot source is p lo t ted i n F igure 8-11. 109 Chapter 8:Test Of Correction > .-. - ••' . Figure 8-11 Activity Contrast in the Hot Spot Source Position Variable Convolution Correction A., Position Invariant Convolution Correction - B - - Uncorrected • * •] • - Q - - 2D Scatter Corrected 'eS • « T3 S S 2 ' (5X1 CS CQ •is © a, . © ' 15 . 20 Direct Plane The prescan contrast measurement of 3.8 is reproduced to varying success by all three scatter correction schemes. The position variable convolution correction at first inspection appears to most closely coincide with the true ratio value on average. Plainly though the variance of this correction is significantly larger than that in either of the other two competing correction schemes. This is unlikely however to originate from correction introduced noise. In fact, the position variable convolution correction's ratio, as a function of transverse plane, closely follows the axial profile of the ratio from the uncorrected image. This suggests that the alternate two corrections may inadvertently have suppressed the underlying variation in ..the ratio. Although further study is required to explain this preliminary result it would seem that position variable, scatter correction reproduces axial contrast variations. One explanation for the position invariant convolution correction's failure to reproduce this same observation may Chapter 8:Test Of Correction ' reside with its use of a spherically symmetric kernel rather than art elliptically symmetric kernel. This may impact on its ability to accurately characterize the axial scatter variation and so mime sharp variations in the underlying activity ratios along this direction. For the case of the hot and cold spot source the two ratios of interest are plotted in Figure 8-12 and Figure 8-13. . •. - " Figure 8-12 Hot Spot to Background Activity Contrast in Hot and Cold Spot Source • Comparison of Correction Performance "« x) s . 3 , •©• . u. two •a , % « . o . o, © a. 4.5 4 3.5 3 .2.5 2" 1.5 | 1 I m J.-rfa ± 1 1^ j a : fa Position Variable Convolution Correction •_— Position Invariant Convolution Correction - B - - Non-Scatter Corrected 10 . 1 5 20 Direct Plane 25 30 111 Chapter 8:Test Of Correction Figure 8-13 Cold Spot to Background Activity Contrast in Hot and Cold Spot Source , Comparison of Correction Performance 10 e3 S a © « CQ © c 2 © -2 | 1 j c > - -R- 1 UOIUWU T Ui ItlMlV V - U H I U I U U U U V-W11^V-LHJ1I 1(111 IU111 .MZA^ V l H l U J C I i a < U G U ; 1 - Position Invariant Convolution Correction 1 J 1 _ , J 1 1 1 t 1 t —< 1 I > < M < \ X 1 >-( > L t \ I 1 1— > C > > ) t <t • • 1 c f - -f ) ' t ( I » < > N N M N ft! 10 25 ' 30 15 20 Direct Plane In both comparisons, the performance of the position variable convolution correction parrots that of the position invariant convolution correction. They both return the expected measured value of 4.0 for the hot to background ratio and 0 for the cold to background ratio. Any Value departing from zero in this last measurement would signify overcorrection by the algorithm. The single key difference that does arise between the performance of the two convolution corrections is the scatter overestimation by the position variable convolution correction in the outer axial planes. The source of this difference should be better elucidated in future investigations into the convolution correction. 112 Chapter 9: Conclusion 9. C O N C L U S I O N The nature of the imaging process in PET results in a significant portion of the collected coincident y-rays to arise from in situ Compton Scattering. The L O R through the scanned object corresponding to these collected events does not reflect the location of a decayed radio-tracer. Including these events in the final reconstruction results in a blurring of the activity contrast levels inside the reconstructed subject. This blurring proves anathema to the program of activity quantitation inside the source. This scatter based error is exacerbated in 3D volumetric imaging. The increased problem motivates an improved scatter correction scheme which takes into account the variation in scatter response from different points inside the object. The program for scatter correction that was therefore developed in this research is an evolution of the 2D integral transformation correction first posited by Bergstrom.42 In the revised variation, a, scatter response is characterized in a 2D projection (reordered into a planagram array) for various projected point activity source positions inside the FOV. The point source constitutes a position localized point source of activity. When it is imaged inside, a 20cm diameter cylindrical water flood, the impulse activity determines a point source scatter response (due to the emitted radiation's interaction with the surrounding water medium). The scatter distribution for each projected point source position is found through a five parameter gaussian fit to the measured activity distributions lying outside the projected point source (which arises entirely from scatter). The fitted parameter values for each projected point source position thus defines a scatter kernel (indexed to that point source projection position). 113 Chapter 9: Conclusion ,\ Interpolation of the kernel determining parameter values across a 20cm transverse and 10.6cm axial square region in the center of the projection planagram establishes a kernel for each sampling element within this last domain. The scatter, correction is based on an iterative position variable convolution of a measured scan with a position variable kernel. The position variability is implemented through the pre-assembly of a distinct kernel for each projection element in the convolution domain. Following acquisition, the raw scan data is reformatted into a planagram array. The subsequent convolution takes the form of multiplying each of the kernel arrays by the projection element pixel (inside a 20cmxl0.6cm central region within the image planagram) value to which the kernel is indexed. The sum of arrays generated from each of these multiplication operations determines a scatter estimate for the input image planagram. For an iterative implementation, this last scatter estimate is subtracted from the initial planagram and the result is used as input to the convolution algorithm again. This iterative cycle then repeats* until the statistical benefit of a succeeding iteration becomes insignificant. When this algorithm was executed on two test activity sources the result of the scatter correction was shown to be of equal or better success to established scatter correction mechanisms. Future and more exhaustive comparisons into the comparative merit of the position variable convolution algorithm can elaborate on these initial findings. In particular effort should be expended on determining the exact nature of the SF sensitivity along the LOR, since the accuracy of such a characterization appears to place an upper limit on the performance of the correction. Such sensitivity will be best 114 Chapter 9: Conclusion elucidated by actively measuring out the SF for a more complete set of point source positions inside the FOV. Furthermore, the performance of the position dependent integral transformation may possibly be improved by the adoption of a deconvolution type correction. This requires the calculation of a kernel h' from h, which acts on the measured projection data directly (rather than on succeeding estimates of the primary activity distribution). The software necessary to generate h' from h has already been developed in the course of this research. However the level of computation as is proves prohibitive without significant optimization and parallelisation in the algorithm. Ultimately the benefit of calculating h' will derive from the elimination of an iterative program for scatter estimations (in the case of K acting as the convolution kernel, a single iteration should reproduce the whole of the scatter component on each planagram). As is, the iterative generation of S by-way of M's convolution with h, is already practical for clinical purposes. O n a Sun Sparc 10 computer, the correction takes only five minutes. This time may be cut even shorter since the calculation may be accurately truncated to only a few iterations (where the number is still to be determined). There appears to be little relative difference in the evaluation of S (with increasing iterations), beyond the first few iterations. 115 • Appendix : ..• • ; Appendix! Comparison of Model Fittings on Central Source Planagram -n=i • view -> (1) Ix = 66 Iz = 16 ' scat_Ix += 66 scat_Iz= 16 scat fraction* = 0.509228 scat_peak_ratio§: = 2.597480e-03 ; ; . X 2 = 1.256502e+00 prob" = 2.623144e-01 data_num=1770 the parameters" are as follows: Specifies the projection peak position along the lx coordinate. ^Specifies the peak position of theScatter distribution in the lx direction. * Specifies the Scatter Fraction ' ." • s Specifies the ratio of the planagram peak height to the scatter peak height. • The probability is given.as cumulative probability for discovering a %2 greater than the value determined. It is determined as the compliment of the incomplete gamma function (which describes the probability distribution of the random variable % v 2, where v=data_num - 5). t + The uncertainty in the parameter values corresponds to the diagnal elements of the %2 hessian (second derivative Taylor expansion abouit the %2 minimum). This value gives the autoyariances (diagnal elemhts of the covariance matrix) at the minimum. The errors therefore constitute a systematic Appendix a[l] = 7.355251e+01 +/-1.071809e+09 a[2] = -8.492778e-02.+/- 1.051736e+09 a[3] = -1.064817e-03 +/.- 1.047024e+09 a[4] = -2.338183e-02 +/- 1.055705e+09 a[5] = -1.752937e-03 +/- 1.053271e+09 n=2 view -> (1) Ix = 66 Tz = 16 scat_Ix = 66 scat_Iz= 16 scat fraction = 0.501616 scat_peak_ratio = 2.420076e-03 X 2= 1.183910e+00 prob = 2.765616e-01 data_num=1770 the parameters are as follows: a[l] = 6.852900e+01+/-• 1.071240e+09 a[2] = -4.954871e-02 +/- 1.050053e+09 a{3] = -1.116408e-03 +/- 1.047234e+09 a[4] =. -1.566652e-02 +/- l.D54900e+09 a[5] = -1.85G733e-03 +/- 1.053305e+09 uncertainty in the determined parameter settings and not an experimental uncertainty of the determined setting to the true minimizing values. 117 Appendix n=3 view -> (1) Ix = 66 Iz = 16 scat_Ix = 66 scat_Iz= 16 scat fraction = 0.494461 scat_peak_ratio = 2.293940e-03 V - 1.143834e+00 prob = 2.848436e-01 data_num=1770 the parameters are as follows: ' a[l] •= 6.495722e+01 +/-1.070824e+09 a[2] = -2.918766e-02 +/- 1.048577e+09 a[3] •= -1.156430e-03 +/- 1.047413e+09 a[4] = -1.072931e-02 •+/- 1.054039e+09 a[5] = -1.952813e-03 +/- 1.053332e+09 n=4 view -> (1) . Ix = 66 Iz = 16 scat_Ix = 66 scat_Iz= 16 scat fraction = 0.487710 scat_peak_ratio - 2.l99530e-03 X 2= 1.126725e+00 prob = 2.884749e-01 ; ' data_num=1770 118 Appendix * the parameters are as follows: a[l] = 6.228382e+01 +/- 1.070534e+09 a[2] = -1.730073e-02 +/- 1.046945e+09 a[3] = -1.188603e-03 +/- 1.047547e+09 a[4] = -7.513261e-03 +/- 1.053235e+09 a[5] = -2.054310e-03 +/- 1.053353e+09 n=5 view -> (1) Ix = 66 Iz = 16 scat_Ix - 66 scat_Iz= 16 scat fraction = 0.481332 scat_peak_ratio = 2.125978e-03 X 2= 1.126529e+00 . prob = 2.885169e-01 data_num=1770 the parameters are as follows: a[l] = 6.020105e+01 +/- 1.070187e+09 a[2] = -1.029644e-02 +/- 1.045485e+09 a[3] = -1.215254e-03 +/- 1.047615e+09 a[4] = -5.375972e-03 +/- 1.052499e+09 119 Appendix a[5] = -2.152969e-03 +/- 1.053369e+09 n=6 view -> (1) Ix = 66 Iz - 16 scat_Ix = 66 scat_Iz= 15 scat fraction = 0.475323 scat_peak_ratio = 2.067315e4)3 %2= 1.139088e+00 prob = 2.858450e-01 data_num=1770 the parameters are as follows: a[l] = 5.852163e+01 +/- 1.069943e+09 a[2] = -6.142816e-03 +/- 1.043900e+09 a[3] = -1.238236e-03 +/- 1.047676e+09 a[4] = -3.922392e-03 +/- 1.051784e+09 a[5] .= -2.248059e-03 +/- 1.053382e+09 120 Appendix Appendix 2 Comparison of Model Fittings on Offset Source Planagram n=l view ->. (4) • rx = 41 Iz = 16 scat_Ix = 41, scat_Iz= 16 ; -scat fraction = 0.245255 scat_peak_ratio = 1.193243e-03 cisquare = 1.130331 e+00 prob = 2.877047e-01 data_num=1611 the parameters are as follows: a[l] = 6.556608e+00 +/- 1.067988e+09 , ' a[2] .=;-5,468943e-02 +/- 1.055301e+09 ' . a[3] =-2.792275e-02+/-1.052379e+09 a[4] = -2.107656e-02 +/- 1.05?.192e+09 ; • . a[5] = 1.035348e-02 +•/- 1.056835e+0?' n=2 view -> (4) ix = 41 . iz = i6 - - .-. 121 Appendix ;.. ". ;•' • scat_lx = 40 scat_Iz= 16 scat fraction = 0.238190 scat_pcak_ratio = 1.1247.00e-03 cisquare = 1.137718e+00 prob - 2.861350e-01 data_num=1611 the parameters are as follows: ' a[l] = 6.170343c+00 +/- 1.067528e+09 a[2] = -3.091754e-02 +/- 1.053743e+09 a[3] =-2.940650e-02+/-1.052432e+09 a[4] = -1.417771e-02 +/-1.058322e+09 . : a [51 = 1.060790c-02 +/- 1.056880e+09 n=3 : ... •• view -> (4) . • " Ix = 41 Iz = 16, ; s.cat_Ix = 35 scat_Iz= 18 scat fraction = 0.231539 scat_peak_ratio = 1.091260e-03 • cisquare = 1.145295e+0'0 prob = 2.845362e-01 data_num=1611 the parameters are as follows: a[l] = 5.893847e+00 +/- 1.067072e+09 122 Appendix ... a[2] = -1.768274e-02 +/-1.052083e+09 a[3] = -3.062286e-02 +/- 1.052476e+0? a[4] = -9.692262e-03 +/-1.057514e+09 a[5] .= 1.072645e-02 +/- 1.056915e+09 n=4 view -> (4) Ix = 41 Iz = 16 scat_Ix = 30' scat_Iz= 19 scat fraction = 0.225293 scat_peak_ratio = 1.078337e-03 cisquare = 1.152806e+00 prob = 2.829626e-01 data_num=1611 the parameters are as follows: a[l] = 5.685357e+00 +/- 1.0667356+09^  a[2] =-1.018446e-02 +/-1.050571e+09 a[3] =-3.166093e-02+/-1.052514e+09 ' ' ' " a[4] = -6.742479e-03 +/- 1.056777e+09 a[5] = 1.073121e-02 +/- 1.056942e+09 123 Appendix n=5 view -> (4) ' Ix = 41 Iz = 16 scatjbc = 26 scat_Iz= 20 scat fraction = 0.219441 scat_peak_ratio = 1.073977e-03 cisquare = 1.160124e+00 prob = 2.814398e-01 data_inum=16.il the parameters are as follows: a[l] = 5.522734e+00 +/- 1.066506e+09 a[2] = -5.895704e-03 +/- 1.048952e+09 a[3] =-3.256214e-02 +/-1.052547e+09. • a[4] = -4.770741e-03 +/- 1.056022e+09 a[5] = 1.066849e-02 +/- 1.056964e+09 n=6 view -> (4) Ix = 41 Iz = 16 scat_Ix = 24 scat_Iz= 21 scat fraction = 0.214010 scat_peak_ratio = 1.073409e-03 cisquare = 1.167185e+00 prob = 2.799803e-01 data_num=1611 124 Appendix the parameters are as follows: a[l] = 5.391014e+00 +/- 1.066286e+09 a[2] = -3.421610e.-03 +'/- 1.047474e+0? a[3] = -3.337300e-02 +/- 1.052577e+09 a[4] = -3.421162e-03 +/- 1.055224e+09 a[5] = 1.056553e-02 +/- 1.056973e+09 [1] Townsend D.W. , "Optimization of Signal in Positron Emission Tomography Scans: Present and Future Developments", C I B A Foundation Symposium 163:57-69 discussion pp. 69-75,1991 [2] Saha, G . B., "Physics and Biology of radionucleisdes",Springer-Verlag,1993 [3] Bailey, D . L . , "Aspects of Quantitation in 3D Positron Emission Tomography",Project Report-March 1992;p.l(ppl-9) [4] Hoffman E. J., Phelps M . E-., "Positron Emission Tomography: and Autoradiography: Principles and Applications for the Brain and Heart" N e w York: Raven Press ©1986 [5] Liang, Z./ 'Detector Response Restoration in Image reconstruction of H i g h Resolution Positron Emission Tomography", IEEE Transactions on Medica l Imaging, June(2) 1994 13(2); pp.314 -319 [6] Bailey D . , Defrise M. ,Kinahan P.E..,Townsend D .W. , "Theoretical and Practical Aspects of Image reconstruction for Volume PET Scanners", IEEE N u c l Science Symposium & Medica l Imaging Conference 1994 [7] Bailey D . L . , Townsend D. W., Jones T., Geissbuhler A: , . Defrise M . , "Implimentation of 125 References Quantitative 3D PET", preprint . [8] Tornai M. P., Germane- G., Hoffman E. J./'Positioning and Energy Response of PET Block Detectors with Different Light Sharing Schemes", IEEE Trans Nucl Med 41(4) 1994;pp.l458-1463 [9] Anger H. O., "Scintillation Camera", Rev. Sci. Instr.29 1958; p.27 [10] Spinks T.}., Jones T., Bailey D. L., Townsend D. W., Grootoonk S., Bloomfield P. M.,Gilardi H. C , Casey H. E., Sipe B., Reed J., "physical Performance, of a Positron Tomograph for Brain Imaging with Retractable Septa", Phys. Med. Biol., vol. 37,1992; pp. 1637-1655 [11] Mazoyer B., Trebossen R., Deutch R., Casey M., Blohm K.,"Physical Characteristics of the ECAT 953B/31: A New High Resolution Brain Positron Tomograph", IEEE Transactions on Medical Imaging 10(4) 1991;pp.499-304 [12] Bailey D, Jones T, Spinks T, "A Method for. Measuring the Absolute Sensitivity of Positron Emission Tomographic Scanners", Eur J Nucl Med 18 1991; pp.374-379 [13] Strother S. C , Casey M., Hoffman E. J., "Measuring PET Scanner Sensitivity: Relating Countrates to Image Signal-to-Noise Using Noise Equivalent Counts", IEEE Trans Nucl Sci 37(2) 1990:; pp.783-788. . . . [14] Bailey D. L., Jones T., Spinks T. J., Gilardi M. C. Townsend D. W., "Noise Equivalent Counts in a Neuro-PET Scanner with Retractable Septa", IEEE Trans. Med. Imag. [15] Knoll G. F., "Radiation Detection and Press", Wiley; p.692 [16] Radon J, (translation by Parks P. C), "On the Determination of Functions from Their Integral Values Along Certain Manifolds", IEEE Trans, on Med. Imag. 5(5) 1986; pp.170-176 [17] Shao L., Karp J. S., Countryman P.,"practical Considerations Of the Wiener Filtering Technique on Projection Data for PET", IEEE Trans, on Nucl Sc. 41(4) 1994;pp.l560-1565 [18] Stark H v "Image Theory and Applications", Academic Press, Inc. 1987 126 References [19] Colsher J G , " " Ful ly Three-Dimensional Positron Emission Tomography", Phys. M e d . Biol . 25(1): 1980; pp. 103-115 [20] Kinahan P E , Rogers J G , "Analyt ic 3D Image Reconstruction Us ing All Detected Events", IEEE Trans, on N u c l . Sci. 36(1) 1989; pp.964-968 ,-• ; [21] Hi l l ie r , R.,"gamma Ray Astronomy",Claredoh Press 1984 [22] Mazoyer B. , trebossen R., Deutch R., Casey M . , Blohm K.,"Physical. Characteristics of the E C A T 953B/31: A N e w H i g h Resolution Brain Positron Tomograph", IEEE Transaction son Medica l Imaging 10(4) 1991;pp.499-504 [23] Casey M., Hoffman E. J., "Quantitation in Positron Emission Tomography: 7. A Technique To Reduce Noise In Accidental Coincidence Measurements and Coincidence Efficiency Calibration", Journal of Computer Assited Tomography.,7 1983;pp .845-850 [24] Harr ison R. L . , Haynor D . R., lewellen T. K. / 'Limita t ions of Energy-based Scatter Correction for Quantitative PET" , Proceedings of the Nuclear Science and Medica l Imaging Conference, 1992;pp.862-864 ' [25] Liang, Z . , "Detector Response Restoration in Image Reconstruction of H i g h Resolution Positron Emission Tomography", IEEE Trani M e d . Im., 2 June 1994; pp.862-864 [26] Thompson. C. J./ 'The Problem Of Scatter Correction i n Positron .Volume imaging", IEEE Transactions O n Medica l Imaging 12(1) 1993;pp. 124-132 [27] Tornai M . P., Germane G. , Hoffman E., "Positioning and Energy Response of PET Block Detectors wi th Different Light Sharing Schemes" Nuclear Science Symposium and Medica l Imaging Conference, 1993 [28] Grootnoonk S., SpinksT.J . , Kennedy A . M . , Bloomfield P. M . , Sashin D. , Jones T. ,"Dual Energy Window Scatter Correction For Positron Emission", Nuclear Science Symposium and Medical Imaging Conference, 1992;pp.942-944 References [29] Spinks T. J., Grootoorik S., Bailey D . L . , Schnorr L . , " A Comparison of 3D Scatter Correction Methods in a neuroPET Tomograph", to appear in the Proceedings of the I E E E / M I C 1994, N o v . . 1994, Norfolk Virg in ia U S A [30] Harr ison R. L . , Haynor D . R„ Lewellen T. K.,"Limitations of Energy-Based Scatter Correction for Quantitative P E T " , Proceedings of the Nuclear Science Symposium 1992;pp.862-864 [31] Shao L . , Freifeldder R., Karp J. S."Triple Energy, Window Scatter Correction Technique in PET",pp.924-926 [32] Bendrium B./ Trebossen R., Frouin, Syrota A . , " A PET Scatter Correction Us ing Simultaneous Acquisitions w i th Lower Energy Thresholds", in press ppl -5 . [33] Barney S., Harrop R., Atkins M . S., " A d d i t i o n of Noise by Scatter Correction Methods in P V I " , IEEE Trans'. N u c l . Sc. 41 1994; pp.1551-1555 : , . [34] Sossi V . , Barney S., Harr ison R. L . , Ruth T:,"Effect of Scatter from Radioactivity Outside-the F O V " , IEEE Trn . N u c l . Sc. 42,4 1995; i n press : : [35] Hoverath H , Kuebler W K , Ostertag H J, D o l l J, Zeigler S I , Knopp M V a n d L o r e n z W J Scatter' Correction In Transaxial Slices of a Whole Body Positron Emission Tomograph", Phys M e d Biol 38 1993; pp.717-28 , ' •  / : ' ' ' ' [36] Bailey D . L . , Meikle S. R . , . " A Convolution-Subtraction Scatter Correction For 3D PET" , Phys: Med. 'Bio l . (39)1994; 411-424 : - V " [37] Mckee B. T., Gurvey A . T., Harvey P. J., Howse D . C , "A Deconvolution Scatter Correction for a -3D Pet Sytem", IEEE Trans. M e d . Im. 11 1992; pp.56G-569 ' <• [38] Shao L . , karp S.,"Cross Plane,Scattering Correction-Point Source Deconvolution in PET" , IEEE Transactions on MedicalTmaging 10(3) 1991;pp.234-239 - : ,, [39] Barney J. S:, Harrop R., Dykstra C. D.,"Source Distribution Dependent Scatter Correction for PVI" , IEEE trans N u c l Sc. 40(4) 1993;pp.l001-1007 , . ' . 128 References [40] Msak i P., Axelsson B., Dahl C. M . , Larrson S. A . , "Generalized Scatter Correction in SPECT using Point Scatter Distribution Functions", J. N u c l . Med . , 1987(28);pp.l861-1869 [41] Press W. H . , Flannery B. P., Teukolsky S. A . , Vetterling W. T. ,"Numerical Recipes in C " , © Cambridge University Press 1988;pp.1-732 - chapter 12 [42] Bergstrom M . . , "Correction For Scattered Radiation", J Comput Assist Tomography 7(1) 1983 [43] Msaki P, Erlandsson K , Svensson L , Nolstedt L , " The Convolution Subtraction Hypothesis and its Val id i ty Domain in Radioisotope imaging" Phys M e d Biol 38(1993) 1359-1370 [44] Links J. M . , Wagner H . N . , "Specification of Performance of Positron Emission Tmography Scanners", (Letter) J N u c l M e d 23:p. 82 1982 [45] Strother S. C. , Casey M . E. , Hoffman E. J.,"Measuring PET Scanner Sensitivity: relating. Countrates to Image Signal-to-Noise Ratios Using Noise Equivalent Counts",IEEE Transactions on Nuclear Science 37(2) A p r i l 1990:pp.783-788 (p.783) [46] Barney S.,"Scatter Correction in Positron Volume Imaging", Ph .D thesis 1993 [47] W u C , Caesar E. C*.,Chen C , "Characterization and Correction for Scatter in 3D PET Using Rebinned Plane Integrals", IEEE Trans. N u c l . Sc. 41(6) 1994; pp .2758-2764 129 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items