UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 ANALYSIS OF TITAN HYBRID ALGORITHMS FOR SPECT SIMULATION By KATHERINE ELIZABETH KELLER ROYSTON A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEG REE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2010 PAGE 2 2 2010 Katherine Elizabeth Keller Royston PAGE 3 3 To my parents for all their encouragement and my husband for his support PAGE 4 4 ACKNOWLEDGMENTS I would like to ex press my gratitude to Dr. Alireza Haghighat, my advisor, for his patience and guidance in this work. His advice enabled me to successfully overcome the many obstacles I encountered. I thank committee member Dr. Ce Yi for taking the time to work with me a nd provide vital assistance in this research. I also thank Dr. David Gilland for his time as a committee member and his suggestions for the completion of this thesis. I would like to thank my fellow graduate students, Mike Wenner and Will Walters, for pr oviding discussion of problems as they arose. Mike especially provided support on many occasions for which I am very grateful. My family deserves acknowledgement for their endless support and encouragement of my academic studies. My parents have taught m e to always do my best and provided me with the skills to succeed. I would like to thank my husband Sean for always having in faith in me. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................. 4 LIST OF TABLES ............................................................................................................ 7 LIST OF FIGURES .......................................................................................................... 8 ABSTRACT ................................................................................................................... 10 CHAPTER 1 INTRODUCTIO N .................................................................................................... 12 Motivation ............................................................................................................... 12 SPECT Background ................................................................................................ 12 NCAT Code ............................................................................................................ 13 2 THEORY ................................................................................................................. 17 TITAN ..................................................................................................................... 17 Deterministic Transport Code ........................................................................... 17 Differencing Schemes ...................................................................................... 19 Adjoint Transport Equation ............................................................................... 20 Monte Carlo ............................................................................................................ 21 Multigroup Approximation ....................................................................................... 23 3 METHODOLOGY ................................................................................................... 25 General Simulation Methods ................................................................................... 25 TITAN Methods ....................................................................................................... 26 SIMIND Methods .................................................................................................... 27 MCNP5 Methods .................................................................................................... 28 Adjoint Methods ...................................................................................................... 30 4 RESULTS ............................................................................................................... 37 TITAN Analysis ....................................................................................................... 38 SIMIND Analysis and Comparison .......................................................................... 39 Investigation into SIMIND and TITAN Differences .................................................. 42 Adjoint Analysis and MCNP5 Comparison .............................................................. 43 5 CONCLUSIONS AND FUTURE WORK ................................................................. 60 PAGE 6 6 APPENDIX A INPUT FILES .......................................................................................................... 62 B DATA PROCESSING CODES ................................................................................ 77 LIST OF REFERENCES ............................................................................................... 93 BIOGRAPHICAL SKETCH ............................................................................................ 95 PAGE 7 7 LIST OF TABLES Table page 4 1 Comparison between myocardium pixels in group 1 anterior projection images with increasing order of quadrature for different number of voxels (meshes) ............................................................................................................ 48 4 2 Maximum absolute relative diff erence between TITAN flux distributions in the myocardium of the 64 x 64 x 64 model relative to the 128 x 128 x 128 model ... 49 4 3 Infinity norm and 2norm for TITAN S6DDZ and S6DTW compared to SIMIND for 131 million particles per projection ................................................... 50 4 4 Three group computation time comparison for four projection angles ................ 50 4 5 Comparison of relative differences of TITAN S6DDZ and S6DTW compared to SIMIND for 131 million particles per projection including the infinity norm, the next largest difference after the infinity norm, and the percentage of pixels with absolut e relative differences greater than 10% ................................. 51 4 6 Single hole collimator responses for TITAN compared with MCNP5 .................. 56 4 7 Single hole collimator responses for TITAN with PNTN ordinate splitting compared with MCNP5 ....................................................................................... 56 4 8 Four hole collimator responses for TITAN compared with MCNP5 .................... 57 4 9 Nine hole collimator responses for TITAN compared with MCNP5 .................... 57 4 10 Ideal and fitted direction weights for each collimator case and the resulting recommended order of quadrature ..................................................................... 59 4 11 Calculation times for 16 processors in MCNP5 and TITAN for each collimator .. 59 PAGE 8 8 LIST O F FIGURES Figure page 1 1 Geometry of a SPECT study .............................................................................. 15 1 2 Cut away view of NCAT material distribution ...................................................... 15 1 3 Cut away view of NCAT source distribution ........................................................ 16 3 1 Example of pixel numbering in projection images ............................................... 33 3 2 Schematic of a fictitious quadrature set with circular ordinate splitting ............... 33 3 3 YZ planar view of MCNP5 simulation geometry with a single collimated detector voxel of air located over the heart ......................................................... 34 3 4 MCNP5 simulation collimators ............................................................................ 34 3 5 Depiction of a photon passing through the minimum path length w from one collimator hole to another and being detected as a function of hole length l hole diameter d and septal thickness t .............................................................. 35 3 6 TITAN geometry for adjoint simulation using SN and CM regions ....................... 35 3 7 Source directions chosen in four octants for S62 used in adjoint calculations to approximate the acceptance angle shown by the circle ................................. 36 4 1 Normalized TITAN projection images at four angles around the patient ............ 48 4 2 Infinity norm and 2norm relative to 131 million photons per projection in SIMIND ............................................................................................................... 49 4 3 SIMIND three energy group anterior projection images ..................................... 49 4 4 TITAN three energy group anterior projection images ........................................ 50 4 5 SIMIND and TITAN normalized projection images of the heart ......................... 50 4 6 Plot of normalized flux along a single profile of anterior projection images created by the SIMIND and TITAN code along wit h the relative difference through A) column 44 and B) row 38 .................................................................. 52 4 7 Plot of normalized flux along a single profile throug h row 38 of SIMIND and TITAN projection images and the relative difference ......................................... 53 4 8 Group 1 adjoint f lux through the heart from the TITAN model of a singlehole collimator ............................................................................................................ 55 PAGE 9 9 4 9 Single hole collimator relative difference in TITAN response compared to MCNP5 as a function of direction weight (order of quadrature) .......................... 56 4 10 Four hole collimator relative error in TITAN response as a function of direction weight (order of quadrature) ................................................................. 57 4 11 Nine hole collimator relative error in TITAN response as a function of direction weight (order of quadrature) ................................................................. 58 4 12 Optimal TITAN direction weight versus the cosine of the collimator acceptance angle in MCNP5 .............................................................................. 58 A 1 NCAT input file ................................................................................................... 63 A 2 CEPXS input file ................................................................................................. 65 A 3 TITAN input file for SPECT simulation ................................................................ 66 A 4 SIMIND parameter file ........................................................................................ 70 A 5 MCNP5 input file for the four hole case .............................................................. 71 A 6 TITAN input file for adjoint methodology ............................................................. 74 B 1 Code to compare binary projection images from the TITAN and SIMIND codes .................................................................................................................. 78 B 2 Program to calculate the detector response using the adjoint flux produced from an adjoint calculation in the TITAN code .................................................... 81 B 3 Program to convert NCAT material distribution into material numbers for MCNP5 input ...................................................................................................... 86 B 4 Linear regression program ................................................................................. 89 B 5 Program to compare flux distributions of two different mesh sizes ..................... 91 PAGE 10 10 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science ANALYSIS OF TITAN HYBRID ALGORITHMS FOR SPECT SIMULATION By Katherine Elizabeth Keller Royston May 2010 Chair: Alireza Haghighat Major: Nuclear Engineering Sciences Traditionally, Monte Carlo methods have been used in the simulation of Single Photon Emission Computed Tomography (SPECT), however computer advances have lead deterministic methodologies to become more feasible than previously. The TIT AN code is a hybrid deterministic transport code that solves the linear Boltzmann equation using several approaches: the discrete ordinates method (SN), the characteristics method (CM), and ray tracing. While the SN method is the most commonly used and ge nerally more efficient, the MOC is best in low scatter regions. By allowing the use of these two methods in different regions of a problem, the TITAN code is more efficiently able to solve problems like SPECT simulation than a code allowing only one of th e methods. In this research, benchmarking of the TITAN code was performed on its application to a SPECT simulation of a myocardial perfusion study. The TITAN results were compared with two Monte Carlo codes: SIMIND and MCNP5. The SIMIND code was used as an established Monte Carlo SPECT simulation code for comparison with projection images created using TITAN. In the heart region, general agreement was found between the SIMIND and TITAN projection images with differences attributed PAGE 11 11 largely to the differences in cross sections used in the two codes. Examination of calculation times showed that for this model the TITAN code is 1.5 to 48 times faster than the SIMIND code, depending on the particular simulation parameters. The TITAN projection images were al so compared for two types of differencing schemes and mesh sizes and several orders of quadrature. Neither of the two differencing schemes stood out as better than the other and the larger 0.625 x 0.625 x 0.625 cm3 mesh size was adequate. The adjoint methodology of TITAN was tested for its ability to represent a collimator by comparing detector response with the MCNP5 code. Several collimator hole sizes were used to examine potential limitations of the TITAN code. In each case, a quadrature that gave les s than a 1% relative difference was found. Also, a linear relationship was demonstrated between the TITAN adjoint source direction weight and the response. A limitation of this adjoint methodology is the increased calculation time for the high orders of quadrature needed to represent increasingly small acceptance angles. However, the TITAN adjoint calculation was still demonstrated to be more than a factor of 3 faster than the MCNP5 response calculation time. PAGE 12 12 CHAPTER 1 INTRODUCTION Motivation Traditional ly, Single Photon Emission Computed Tomography (SPECT) simulations have been conducted using Monte Carlo methods.1 However Monte Carlo methods have the drawback of requiring large amounts of calculation time in general compared to deterministic methods.2 Recently, the hybrid deterministic transport code TITAN has been applied to the simulation of SPECT3 with the hope that it can be used to develop methods of improving image quality and reducing patient dose. The immediate goal of this research is the benchmarking of the TITAN code for the simulation of a myocardial perfusion study. The results of the TITAN code will be compared with those of the SIMIND and MCNP5 codes. The effects of variations in mesh size, quadrature order, and differencing scheme on the results of the TITAN code will also be investigated. SPECT Background Nuclear medicine consists of two main imaging modalities: Single Photon Emission Computed Tomography (SPECT) and Positron Emission Tomography (PET). The SPECT images use radionuclides that decay by gamma ray emission, whereas PET uses positron emitting radionuclides. A general nuclear medicine study involves injecting a radiopharmaceutical into the patient. The radiopharmaceutical consists of a radionuclide tagged to a pharmaceutic al that is designed to be preferentially absorbed by the region of interest. The radiation emitted from the patient can then be used to provide a functional image. Projection images of the patient are obtained using gamma cameras positioned at various angles around the patient, as shown in Figure 1 1. If PAGE 13 13 they have been obtained at enough angles, these projection images can then be used to reconstruct a 3dimensional image of the patient. The acquisition of a projection image in SPECT requires the use o f a collimator in front of the gamma camera to provide spatial resolution. A collimator consists of a slab of a high z material, typically lead, with holes that are usually parallel. The collimator restricts photons that are counted to photons from the phantom along the detector normal from each collimator hole, within an acceptance angle. Rejection of photos outside of an energy window then reduces scatter counts. Unfortunately, these techniques also greatly reduce the number of counts at the detector and therefore increase the length of time needed to acquire an image of adequate quality. Nuclear medicine images are obtained for many diagnostic studies including the determination of coronary artery disease, neurological disorders, and cancers. The ra dionuclide most commonly used in SPECT is 99mTc, which is frequently applied to measuring myocardial perfusion through the radiopharmaceutical Sestamibi. Myocardial perfusion studies determine if the myocardium, or heart wall, is receiving enough blood su pply and can be used to diagnose coronary artery disease based upon how much of the radionuclide is absorbed into the myocardium. In this research, a myocardial perfusion study is simulated by placing a uniform source in the heart wall of the phantom. NCA T Code To provide a phantom of a patient for the SPECT simulations, the NURBS based cardiac torso (NCAT) code4 was used. The NCAT code uses nonuniform rational B splines (NURBS) to generate the phantom organs based upon data from the Visible Human Projec t. An input file is used to specify all the details of the desired phantom PAGE 14 14 and the NCAT code then outputs two distributions: a material distribution and a source distribution. The material distribution contains the linear attenuation coefficient of each of the materials in the phantom and the source distribution contains the relative intensity of each voxel. For the simulations completed in this work, a 40 x 40 x 40 cm3 phantom simulating an upper torso was generated using the NCAT code with default input s. Two voxel sizes were used: 0.625 x 0.625 x 0.625 cm3, which results in 64 x 64 x 64 voxels in the phantom, and 0.3125 x 0.3125 x 0.3125 cm3 for a 128 x 128 x 128 voxel phantom. Thirteen different materials were represented including: air, water (body) intestine, lung, rib bone, liver, kidney, spine bone, muscle, cartilage, spleen, heart, and blood. Figure 1 2 is a visual representation of the phantom with the top cut away so that the internal organs are visible. In the source distribution, a 140.5 k eV source was placed in the heart wall to simulate uptake of 99mTc with an intensity of 75 specified in the NCAT input file. A far less intense source with a specified intensity of 2 was placed in the body, muscle, and blood to account for some radiopharm aceutical uptake of these materials as shown in Figure 13. PAGE 15 15 Figure 1 1. Geometry of a SPECT study Figure 12 Cut away view of NCAT material distribution PAGE 16 16 Figure 13. Cut away view of N CAT source distribution. The red is the intense source in the myocardium and the light blue is a less intense source. PAGE 17 17 CHAPTER 2 THEORY TITAN The TITAN code is a hybrid deterministic transport code that solves the linear Boltzmann equation (LBE).5 The cod e is a hybrid because it allows the use of different solvers of the LBE in different regions (or coarse meshes) of a problem: a Discrete Ordinates (SN)2 solver, a Characteristics Method (CM)6 solver, and a ray tracing solver.6 The TITAN code is especially applicable in problems with low scattering regions, as is found in the SPECT problem, where the CM and ray tracing solver can be used to speed up the calculations. Deterministic Transport Code Deterministic codes seek to numerically solve the LBE for neut ral particles, shown in time independent integrodifferential form in Equation 21 for a non multiplying medium.7 The LBE is derived from a balance of particles in a phase space as given by d3rdEd H S in V (2 1) where, H t( r E ) d E d 40 s( r E E ) In Equation 21, is the angular flux, t is the total macroscopic cross section, s( r E E ) is the differential macroscopic scattering cross section and S is the independent source. The first two terms on the left side of Equation 21 are the loss from the phase space and called the collision and streaming terms, respectively. The third term on the left side and the right hand side of Equation 21 are the gains to the PAGE 18 18 phase space: the scattering and independent source terms, respectively. T o solve this equation, the phase space can be discretized. Energy is generally discretized using the multigroup approximation7 in deterministic methods and the angular domain is discretized into discrete ordinates or SN.2 In the TITAN code, the spatial d omain is discretized using either meshes for the Discrete Ordinates method (SN), or arbitrarily shaped regions for the Characteristics Method (CM) and ray tracing. The differential scattering crosssection is also expanded in terms of Legendre polynomials Pl( o), where o = as shown in Equation 22. s( r E E ,o) 2 l 1 4 sl( r E E ) Pl(o)l 0 (2 2) Equation 22 is then truncated at some order L according to the degree of anisotropy in the particular problem of interest. In regions defined by the user the TITAN code enables the use of the SN method to solve the LBE. The SN method solves Equation 21 for the angular flux in a discrete set of directions.2 N is called the quadrature order and in a 3D geometry there are N ( N 2 ) / 8 directio ns per octant. Each direction is weighted to preserve the moments of angular flux (e.g. flux and current), and the combination of the directions and their weights is called the quadrature set. Differencing schemes must then be used to deal with the strea ming term and solve the LBE for the quadrature set. The CM is especially useful in low scattering regions. The same angular quadrature method used to discretize the angle in the SN method is usually used, however instead of using spatial meshes, the LBE i s solved along parallel directions called characteristic rays.5 The integral transport equation can then be used to calculate the angular flux along a characteristic ray. A uniform fine mesh is built on the PAGE 19 19 coarse mesh boundaries and characteristic rays are drawn along the quadrature directions from the centers of the fine meshes.6 The CM does not require a differencing scheme to solve the transport equation, but can be very memory intensive. Differencing Schemes If the SN method is being used to solve t he LBE, a spatial approximation is required and the problem geometry is divided up into cell volume elements. The TITAN code allows the use of two spatial differencing schemes: diamond differencing with zero fix up (DDZ)8 and directional thetaweighted (D TW).9 The diamond differencing method (DD) simply uses a linear average of two opposite boundary angular fluxes to determine the cell average angular flux. Equation 23 gives an example of this formulation for 3D Cartesian coordinates along the x axis wh ere gnijk is the cell average angular flux in group g direction n and in spatial mesh ijk gnijk1 2 gn i 1 2 jkgn i 1 2 jk (2 3) In Equation 23, gn,i 1/2,jk and gn,i+1/2,jk are the boundary angular fluxes (or incoming and outgoing fluxes) on either side o f the cell average angular flux. Similar equations are used along the y axis and z axis. In a typical problem, the incoming flux will be known from boundary conditions so that the outgoing and cell average angular flux can be solved for. As can be seen f rom Equation 23, it is possible to obtain negative outgoing angular flux values from the DD scheme, especially in regions with steep flux gradients. One method of dealing with this is DDZ and entails setting the negative flux value to zero and recalculat ing the cell average angular flux before continuing with the calculation. This method is expensive PAGE 20 20 and leads to numerical difficulties.8 Hence, Petrovic and Haghighat developed a positive method called DTW (directional theta weighted). The DTW scheme doe s not equally weight the incoming and outgoing boundary angular fluxes, as the DDZ scheme does. Instead the fluxes are weighted as a function of space, energy, and angle. Equation 24 gives the weighting scheme for 3D Cartesian coordinates along the x ax is for a cell, where the weighting coefficient (for the DTW scheme), a is calculated according to Equation 25. aout x ( 1 a )in x (2 4) a 1 n x in x(n)n y in yn z in z qn( 2n y 2n z )in x (2 5) (n) n 2 (2 6) Similar weights of b and c are calculated along t he y axis and z axis and all three weights depend on a parameter, given in Equation 26 for the x axis. In Equation 25, n, n and n are the direction cosines for direction index n and qn contains the scattering and independent source terms. Adjoint Transport Equation The adjoint equation is derived using the adjoint property7 given in Equation 27, where symbolizes adjoint and symbolizes integration over all independent variables. H H (2 7) In Equation 27, H is the transport operator, H is the adjoint transport operator, and is the adjoint flux. The adjoint identity requires the boundary condition that the flux of incoming particles and the adjoint flux of outgoing particles are both zero (vacuum PAGE 21 21 bou ndary condition). The adjoint form of the transport equation for the timeindependent and nonmultiplying case is given by Equation 28.7 H S in V (2 8) where, Ht( r E ) d E d 40 s( r E E ) In the adjoint equation s( r E E ) is the scatt ering cross section from energy E to E and direction to S( r E ) is the adjoint source, and the rest of the terms are defined as in the LBE. If the adjoint source is defined as a detector cr osssection, the adjoint flux will be a measure of the contribution, or importance, to that detector response. Monte Carlo Monte Carlo is a method that involves simulating each particle history separately and calculating some average behavior usually using tallies, in contrast with deterministic methods that seek to solve the transport equation for an average value. The method is highly dependent on having a good random number generator to sample the probability density functions (PDFs) of the particles i nteraction types, path length between interactions, scattering angle, etc. In general, Monte Carlo codes may require less computer memory than deterministic codes but are far more computationally expensive to complete calculations with an acceptable relat ive error.2 Monte Carlo methods allow for continuous treatment of the position, energy, and angle variables and are generally more accurate. As a result of having continuous energy, continuous (or point wise) energy cross sections can be used, rather than the approximated multigroup cross sections used in deterministic codes. Monte Carlo methods also have the PAGE 22 22 advantage of being naturally parallel since each particle transport is an independent simulation. The term analog Monte Carlo refers to methods that do not use any variance reduction techniques. In nonanalog Monte Carlo, variance reduction is used to speed up the problem by biasing the PDF to increase the probability that a particle contributes to the desired quantity. The outcome is then prese rved by applying a weight (w) to the biased result as shown in Equation 29. wbiased pdfbiased wunbiased pdfunbiased (2 9) Non analog Monte Carlo is frequently necessary to obtain a precise result in a reasonable amount of time. As long as the particle weights are conserved, as in Equation 29, nonanalog Monte Carlo can produce accurate results, however it must be used with caution to avoid unintentional biasing of results. In Monte Carlo methods, the statistical error corresponding to the final result is of the utmost importance to know how precise the result is. S21 N 1 ( xi x )2 i 1 N (2 10) The sample variance is calculated according to Equation 210, where x is the sample mean and N is the number of histories.1 The sample variance of the mean is then given by the Central Limit Theorem1 to be Equation 211. S x 2 S2N (2 11) The relative error can then be estimated using Equation 212. R S x x (2 12) PAGE 23 23 Multigroup Approximation In order to solve the LBE, deterministi c methods use the multigroup approximation to discretize the energy variable into energy groups ( g =1, G ). Within an energy group, the multigroup approximation assumes that all particles interact with the same probability. Equation 2 13 gives the multigroup transport equation.7 g( r )g( r )g( r ) d s g g( r ) g ( r )4 g 1 g Sg( r ) (2 13) In Equation 213, the subscript g denotes the energy group and the angular flux for group g is given by Equation 214. g( r ) dE( r E )EgEg 1 (2 14) The multigroup cross sections are defined in Equations 215 and 216. g( r ) 1 g( r ) dEt( r E )( r ,E )EgEg 1 (2 15) s g g( r ) 1 g ( r ) dE d E s( r E E )( r E )E g E g 1EgEg 1 (2 16) Finally, the multigroup source is defined in Equation 217. Sg( r ) dE S ( r E )EgEg 1 (2 17) By convention, the highest energy group is group 1 and the lower groups follow down to group G The differential scattering cross section s g g( r ) is expanded in the orthogonal Legendre polynomials Pl(o)2, where o in Equation 2 18. PAGE 24 24 s g g( r ,o) 2 l 1 4 s l g g( r ) Pl(o) l 0 (2 18) Using the orthogonal properties of the Legendre polynomials, Equation 219 can then be written. s lg' g( r ) dEEgEg 1 d E sl( E E )( E )E g E g 1 d E ( E )E g E g 1 (2 19) In Equation 219, is the scalar flux and weights the differential scattering cross section. To generate the multigroup cross sections, the CEPXS code10 was us ed. CEPXS was designed to create multigroupLegendre cross sections for use in coupled electronphoton transport calculations with the ONEDANT code. CEPXS calculates the cross sections of arbitrary materials using the cross sections of the constituent el ements according to Equation 220. NiAiA i i wii i (2 20) In Equation 220, Ni is the number fraction of element i Ai is the atomic weight of element i wi is the weight fraction of element i and I is the macroscopic cross section of element i The CEPXS code uses a zero order weighting function in Equation 219, because it does not have any knowledge of the scalar flux. Equation 219 is then written as Equation 221. s lg' g( r ) dEEgEg 1d E sl( E E )E g E g 1E g 1 E g (2 21) PAGE 25 25 CHAPTER 3 METHODOLOGY In order to simulate a myocardial perfusion study in each of the codes, some code specific modifications had to be made. The following sections describe the methods used to perform a variety of SPECT simulations in the TITAN, SIMIND, and MCNP5 codes. General Simulation Methods Each single photon emission computed tomography (SPECT) simulation produces projection images like those that would be obtained by a gamma camera placed at several angles around a patient injected with 99mTc. Four projection angles were simulated in bo th the SIMIND and the TITAN codes: anterior, left lateral, posterior, and right lateral. In MCNP5, only the anterior view was simulated due to time limitations. Projection images plot pixel (1,1) in the upper left corner as depicted in Figure 31. The m ultigroup calculations used three energy bins and cross sections generated using the CEPXS code with ICRP Publication 8911 data on elemental composition of body tissues. Technetium 99m emits a 140.5 keV gamma ray and SPECT systems typically have an energy resolution of approximately 10%12, so the group 1 bounds were chosen to be 10% above and below the source energy (126.45154.55 keV). The second group range was chosen to have the same width as the first (98.35126.45 keV) and the third group then extended to the cutoff energy of 10 keV. Projection images were normalized to the highest pixel value. In the projection images, it is clear that the peripheral pixels will receive few counts. These low count pixels are in unimportant regions far from the sour ce and will necessarily have large errors. Also, in a myocardial perfusion study only accuracy in the heart is sought and PAGE 26 26 other regions of the body can be disregarded. For these reasons, numerical comparison between projections was confined to pixels wit h normalized values greater than 0.65 in the heart region. This cutoff value was chosen to restrict the analysis to the high count pixels in the myocardium. Quantitative comparisons between different projection images at the same angle were made using the infinity norm and 2 norm given by Equations 31 and 3 2. Infinity norm = max[( P ij Pij) / Pij] (3 1) 2 norm = P ijPijPij 2 i, j 1 / 2 (3 2) In Equations 31 and 3 2, Pij is the pixel value in position ( i,j) of the projection image. TITAN Methods The TITAN code allowed for di rect use of the NCAT codes material and source distributions as well as the multigroup cross sections generated using the CEPXS code. The flux distributions of all energy groups are solved for in one calculation and a projection image is output for each energy group at each specified projection angle. In TITAN SPECT simulations, the SN method was used to solve the LBE within the NCAT phantom and the ray tracing method was used to transport photons from the phantom surface to the detector. For each projection angle, TITAN uses a fictitious quadrature technique to calculate the boundary angular flux along the projection angle.3 The detector is not directly simulated in the TITAN code, however the effects of a collimator are incorporated using an ordinate s plitting (OS) technique.13 The OS technique involves further splitting a chosen direction in a quadrature set. Ordinate splitting can be applied when a highly peaked angular dependent flux or source is PAGE 27 27 present in a problem. To simulate a SPECT collimator, the TITAN code implements the circular ordinate splitting technique3 in which the split directions are located on the edge of a circle centered on the original projection direction as depicted in Figure 32. The SPECT collimator acceptance angle determ ines the radius of the circle. The angular flux reaching the detector location is then averaged over the angular fluxes of the original and split directions. In this way, collimator blur is accounted for in the simulated projection image. A limitation o f this method is that interactions with the collimator septa are not simulated. The sensitivity of the results of the TITAN code to quadrature order, differencing scheme, and mesh size was investigated with the 64 x 64 x 64 and 128 x 128 x 128 voxel NCAT p hantoms. The TITAN code allows the user to choose either the DDZ or DTW differencing scheme. These differencing schemes were used in conjunction with different quadrature orders of S4, S6, S8, S10, and S12, which correspond to 3, 6, 10, 15, and 21 direct ions per octant, respectively. A brief comparison of the flux distributions in the 64 x 64 x 64 and 128 x 128 x 128 phantoms is given in Chapter 4. SIMIND Methods The SIMIND code is a Monte Carlo code designed specifically for SPECT simulation.14 The cod e is able to directly use the NCAT material and source distributions for the phantom, however it does not use any sort of provided cross section file. An input file allows the specification of many parameters, including collimator details. To match the T ITAN simulation as closely as possible, the modeling of a collimator was turned off in the SIMIND code and an acceptance angle equal to that used in TITAN (0.66) was specified. To obtain results for each energy group, the code was executed separately wit h an energy window set to contain the energy group PAGE 28 28 desired. As output, the SIMIND code produces a projection image for each angle that has been simulated. A major drawback of the SIMIND code was the fact that it does not provide any statistical output. B ecause of this, the first step with using this code was to determine the appropriate number of simulated source particles needed to converge on a solution. To achieve this, calculations with increasing numbers of particles were performed and compared using an infinity norm and a 2norm. Once it was determined that the code was converging on a solution, the number of particles needed to obtain an appropriate level of precision could be found. The projection images produced by the SIMIND code were then com pared with those created in TITAN as described at the beginning of this Chapter. MCNP5 Methods The Monte Carlo N Particle, Version 5 (MCNP5) code is a general purpose Monte Carlo code for neutron, photon, and electron transport developed by Los Alamos National Laboratory.15 The MCNP5 model used a lattice of universes containing different materials, ultimately resulting in the same 40 x 40 x 40 cm3 phantom that the SIMIND and TITAN codes read in. Only the 64 x 64 x 64 voxel phantom was used in the MCNP5 si mulations. A single 0.625 x 0.625 x 0.625 cm3 detector voxel of air was positioned over the heart and 10 cm from the phantom surface with a collimator placed in front of it. Figure 33 shows the MCNP5 geometry. In the MCNP5 code, the collimator was model ed as a lead, parallel hole collimator with circular holes and a depth of 5.7 cm. Three collimator hole configurations were used and these are displayed in Figure 34. There were singlehole, four hole, and ninehole collimators with acceptance angles of 3.14, 1.42 and 0.98, respectively, in PAGE 29 29 the 0.625 x 0.625 cm2 face of the lead. The acceptance angle was defined as tan 1( r / l ) in which r is the collimator hole radius and l is the collimator hole length. To effectively eliminate any stray photons from reaching the single voxel detector by traveling through the collimator walls or outside of the collimator, a 0.625 cm thickness of lead was added around the collimator and detector. The minimum thickness between adjacent holes was calcu lated to limit septal penetration to less than 5%.12 Figure 35 depicts a photon passing through a collimator septum with the shortest path length possible to reach the detector, w Using this figure, an approximate relationship can be made between the s eptal thickness t the collimator hole diameter d and the collimator hole length l t w 2 d l w (3 3) To limit the septal penetration to 5%, exponential attenuation was used to obtain Equation 34. w 3 / (3 4) In Equation 3, is the linear attenuation coefficient of the septal material. Placing Equation 34 into Equation 33 and rearranging then gives Equation 35. t 6 d /(l 3 ) (3 5) The detector response R (counts/s) of a detector with cross section d was cal culated from the flux distribution according to Equation 36. R dEdV( r E ) S( r E ) g id g, i Vi g 1 G Vi Vd (3 6) In Equation 36, is the scalar flux, S is the adjoint source density, Vd is the detector volume, and Vi is the volume of cell i The adjoint source density was chosen to be 1 PAGE 30 30 in each energy group for simplicity, because the purpose of this study is to compare the responses not obtain the exact value. A flux tally was used in the detector volume to obtain the scalar flux for each energy group. Adjoint Methods A djoint calculations were used to compare the detector response from the phantom in the TITAN code simulations to the response determined from the MCNP5 code described in the previous section. These comparisons were made for a single voxel detector, which was the adjoint source. Figure 36 displays this configuration with the SN and CM regions used in the TITAN code labeled. The purpose of this comparison was to develop and validate an adjoint method of representing the collimator in the TITAN code. The deterministic TITAN code only solves the LBE for a limited number of directions. However, since the detector response is limited to particles approaching within the acceptance angle of the simulated collimator, the adjoint source needed to be restricted to within this same angle. To achieve this, angles much smaller than those found in S12 were needed and therefore high orders of quadrature were used. Some experimentation was done to determine the best quadrature choice. In each case, a quadrature was us ed that contained only one direction within the acceptance angle along the +y axis. The adjoint source was then defined along that direction in the four octants that face the +y axis. Figure 3 7 demonstrates how these four directions approximate the col limator acceptance angle. The ordinate splitting technique was also used to determine if increasing the number of directions within the acceptance angle improved the accuracy of the response. PAGE 31 31 The TITAN adjoint model of collimation is fundamentally different from the MCNP5 collimator model. However, it is interesting to consider the area subtended by the solid angle for comparison with the area of the opening in the MCNP5 collimator holes. The quadrature direction weight, w is directly related to a nd the area, A, as shown in Equations 37 through 3 9. w 2 (3 7) A r2 (3 8) w 2 A r2 (3 9) In Equations 37 through 39, r is the collimator thickness and the /2 is a normalization factor for the s olid angle of an octant. Since the source will be specified in four octants, the total subtended area can finally be found according to Equation 310. A 4 A 42 r2w (3 10) By using this linear relationship between area and direction weight, the relative difference in the response of the TITAN for different orders of quadrature can be viewed as a function of area. The adjoint calculation in the TITAN code results in one adjoint flux distribution for each energy group. The response of a detector placed at the location of the adjoint source can then be calculated according to Equation 311. R dEdV( r E ) S ( r E ) d g, i Sg i Vi g 1 G V i V s (3 11) In Equation 311, is the adjoint flux and S is the source density. The MCNP5 response with collimation will have an entrance surface area into the collimator that is PAGE 32 32 determined by the number and size of the collimator holes. However, the TITAN codes adjoint source will be emitted from the entire surface of the detector voxel. To account for this difference, a ratio of surfac e areas was multiplied by the TITAN response. PAGE 33 33 Figure 31. Example of pixel numbering in projection images Figure 32. Schematic of a fictitious quadrature set with circular ordinate splitting PAGE 34 34 Figure 33. YZ planar view of MCNP5 simulation geometry with a single collimated detector voxel of air located over the heart A B C Figure 34. MCNP5 simulation collimators A) Singlehole collimator with 3.14 acceptance angle B) Four hole collimator with 1.42 acceptance angle C) Nine hole collimator w ith 0.98 acceptance angle PAGE 35 35 Figure 35. Depiction of a photon passing through the minimum path length w from one collimator hole to another and being detected as a function of hole length l hole diameter d and septal thickness t Figure 36. TITAN geometry for adjoint simulation using SN and CM regions PAGE 36 36 Figure 37. Source directions chosen in four octants for S62 used in adjoint calculations to approximate the acceptance angle shown by the circle PAGE 37 37 CHAPTER 4 RESULTS In this section, the previously described methodologies are put into practice. The projection images produced by TITAN in a SPECT simulation are compared for increasing orders of quadrature for two differencing schemes and mesh cases. The SIMIND code is analyzed to verify that enough particles are being simulated and that the solution is converging. The projection images from the SIMIND and TITAN codes are compared both visually and quantitatively for four projection angles. Differences between images produced for each energy group are also examined. For several chosen cases, the computation times of the TITAN and SIMIND codes are compared. All of the above calculations were completed in serial using the PC cluster designated Einstein. Einstein contains eight nodes with two process ors per node. Each processor is a 2.4 GHz AMD Dual Opteron processor and each node contains 4 GB of DDR SDRAM. For the adjoint methodology in TITAN, the detector responses from several different orders of quadrature are compared to the forward calculati ons performed by the MCNP5 code for each collimator case. This allows an optimal order of quadrature for a particular collimator simulation to be found for each case. An attempt is then made to develop a methodology to recommend the best order of quadrat ure to use in the TITAN code, for a given collimator acceptance angle from the MCNP5 collimator models. Finally, the calculation times for each collimator case are compared. These calculations were performed on the PC cluster designated Dirac. Dirac con tains 16 nodes with eight processors per node. Each processor is a 2 GHz Intel Xeon Quad Core processor with 4 GB of DDR2 SDRAM. PAGE 38 38 TITAN Analysis Normalized group 1 projection images created by TITAN at four different angles around the patient are displayed in Figure 41. The sensitivity of TITANs results to varying orders of quadrature for two differencing schemes and mesh sizes was examined. Anterior projection images created using the TITAN code were compared for the diamond differencing with zero fix up (DDZ) and directional theta weighted (DTW) differencing schemes and 64 x 64 x 64 and 128 x 128 x 128 mesh cases. Results are given in Table 41 for each projection image compared to the next higher order of quadrature. The infinity norm and 2norm val ues demonstrate that for either mesh size, the DTW differencing scheme converges to its solution at a lower order of quadrature than the DDZ method. It can be seen that both differencing schemes maintained similar infinity norms between the 64 x 64 x 64 a nd 128 x 128 x 128 mesh cases. The infinity norm values, regardless of the particular case, indicate that a low order of quadrature is sufficient for this simulation. The differences in the flux distributions for the 64 x 64 x 64 (largemesh) and 128 x 128 x 128 (small mesh) phantoms in TITAN were examined. Since each voxel in the large mesh phantom will correspond to 8 voxels in the small mesh phantom, two steps were taken to compare these cases. First, only voxels with the heart wall material in them w ere compared. This means that for a single largemesh voxel there could be anywhere from 0 to 8 small mesh voxels compared to it. Second, the relative difference for a largemesh voxel was defined as the average of the absolute relative difference between the largemesh voxel and each corresponding small mesh voxel meeting the first step. The maximum absolute differences are given in Table 42 for each differencing scheme and quadrature order used. No consistent pattern of increasing or PAGE 39 39 decreasing absol ute relative differences is seen with increasing order of quadrature. The DDZ and DTW schemes also do not show any pattern to indicate that one of them has higher differences between the two voxel sizes. The differences are all small (<0.2%) and the same order of magnitude, which indicates good agreement between the two models regardless of differencing scheme and quadrature order. SIMIND Analysis and Comparison In order to determine the number of particles needed to converge on a solution with no statist ical data from SIMIND, several calculations were performed with an increasing number of particles simulated. Figure 42 shows the infinity norm and 2norm relative to a run of 131 million photons per projection as a function of increasing number of photons per projection image. It is clear from the trend of the plot that the SIMIND codes solution is converging. The appropriate number of photons per projection needed will depend on the level of precision desired Visual comparison of projection images pro duced by the SIMIND and TITAN codes revealed some notable differences in the lower energy groups. Figures 43 and 4 4 display the anterior projection images produced by the SIMIND and TITAN codes respectively for each of the three energy groups simulated. Due to photoelectric absorption, dense and high Z materials, like bone, will greatly reduce the number of photons reaching the detector, especially at low energies. In the TITAN projection images this effect is clearly visible, however it is not so apparent in the SIMIND code images. Figures 43B and 43C show the SIMIND projection images for energy groups 2 and 3 with the expected attenuation in bone not visible and the photons appearing to have spread out in the phantom. The group 3 SIMIND image shows a faint shadowing of the rib cage. The TITAN energy group 2 projection image, Figure 44B, shows a light PAGE 40 40 outline of the rib cage and Figure 44C shows the TITAN group 3 result with the clavicle, sternum and ribs clearly visible. This matches well with t he drastic increase in attenuation that is expected with decreasing photon energy due to the photoelectric effect. As mentioned previously, the SIMIND input does not include any cross sections and so its calculations are based on its knowledge of the material linear attenuation coefficients. The limited information SIMIND works with appears to leave the code unable to accurately represent energy groups lower than the first group containing the source (i.e. groups made up of scattered particles). From these results it was determined that further analysis with the SIMIND data shall be performed for the first energy group only. Since a SPECT study desires images without scatter, the first energy group is the one of primary interest and this limitation is ac ceptable for this study. The results of the SIMIND and TITAN codes were compared using the infinity norm and 2norm of the normalized projection images at four difference angles around the phantom. The most accurate SIMIND result (131 million particles per projection) and the TITAN code results for S6DTW and S6DDZ were chosen for comparison. The choice of S6 for order of quadrature in the TITAN code had an infinity norm of about 1%, as found in Table 41. Figure 4 5 displays the normalized projection i mages of the heart produced by each code for a visual comparison and the images are in good agreement. The numerical comparison can be found in Table 43 and shows that comparisons with S6DDZ and S6DTW were similar. The anterior projection images show the smallest infinity norm of 15% and the right lateral images have the largest infinity norm of 49.2%. The heart is closest to the anterior surface of the body and so this projection is expected PAGE 41 41 to be more accurate than other angles in SIMIND. Comparing the DDZ and DTW results reveals that neither method is better than the other for this model. A likely contributor to the significant differences between the TITAN and SIMIND images is the fact that SIMIND uses attenuation values provided in the NCAT phantom, while cross sections were generated separately for TITAN. Further, the large difference in the right lateral image could be attributed to a higher expected uncertainty for the SIMIND calculations at that particular angle. The computation times for each code to create four projection images are compared in Table 44 for SIMIND with 131 million and 1 million photons per projection and for TITAN S6 using DDZ and DTW. All calculations were performed in serial using the Einstein PC cluster. A speedup f actor is defined relative to the SIMIND case with 131 million photons per projection. Table 44 clearly shows that the TITAN code was able to complete the calculations in significantly less time than SIMIND; however, the time is highly dependent on the number of photons per projection simulated in SIMIND. Even for the relatively fast 1 million photons per projection SIMIND case, the TITAN code took approximately a third less time. The speedup of the TITAN code over the SIMIND code will further improve if the number of projection angles simulated is increased.3 This is because TITAN has the advantage of finding a solution for the angular flux distribution in the phantom once and then creating however many projection images are required. In contrast, the SIMIND code must repeat the simulation for each projection image. The DDZ scheme was 30 s faster than the DTW method, which could become a reason to use DDZ for larger problems since their accuracy appears to be equal for this case. However, this is for comparing the methods PAGE 42 42 using the same mesh size and, in principle, the DTW scheme yields accurate solutions for larger mesh sizes, for which the DDZ scheme is not accurate, and therefore requires a larger number of meshes and more computation time. For the present study, very small mesh sizes were needed to properly represent the heart wall in the phantom. For these mesh sizes the DDZ scheme is as accurate as DTW, and somewhat more efficient. Investigation into SIMIND and TITAN Differences The differences seen between the TITAN and SIMIND code results in the previous section were significant. Further analysis was completed to better understand these differences. First relative differences between normalized projection image pixels greater than 0.65 were c onsidered in greater detail. The fraction of these pixels with absolute relative differences greater than 10% is given in Table 45. The second largest relative difference is also found in Table 45. The results show drops of 10% and 20.2% to the next l argest relative difference in the right lateral projection images for the S6DDZ and S6DTW cases, respectively. This indicates that the differences between the two images are overstated by the infinity norm. The majority of pixels have a relative differ ence of less than 10% for the anterior, left lateral, and right lateral projection images. The posterior projection images have differences greater than 10% for most pixels. Single profiles of pixels, either columns or rows, were plotted from the SIMIND a nd TITAN S6DDZ projection images. These plots display the behavior of the flux through the projection images. Profiles were chosen that passed through the heart and are plotted in Figures 46 and 47 along with the relative difference for each pixel. F igure 46 plots a column and a row through the heart for the anterior projection angle. It is clear PAGE 43 43 that the SIMIND and TITAN plots have the same shape; however there appears to be a shift between the two data sets in the column plot (Figure 46A). Attem pts to discover a cause for this shift have not yet been successful, but it should be noted that the edges of the NCAT heart phantom line up exactly with the position of the edges in the TITAN anterior projection image. This would indicate that the SIMIND data is shifted. Further studies using the MCNP5 code are planned and should help determine if there is actually a shift. Profile plots for a row through the heart of the left lateral, posterior, and right lateral projection images are found in Figure 47. Figure 4 7A shows good general agreement with a difference in plot shape seen around pixel 32. While the overall behavior is similar in the posterior image, Figure 47B, there appear to be differences in magnitude. The right lateral image, Figure 47 C, again shows that the behavior of the two plots matches very well; however a shift between the data points is visible again. Adjoint Analysis and MCNP5 Comparison The single hole collimator with a 0.625 x 0.625 x 0.625 cm3 voxel detector located 10 cm from the phantom surface and centered over the heart had a response of 0.07493 2.13% in MCNP5. Relative differences for the TITAN response calculations were then found relative to this value. MCNP5 stated uncertainties are all one standard deviation. The singlehole MCNP5 collimator had an acceptance angle of 3.14 and several orders of quadrature were used to calculate the response in TITAN to determine the best match to this acceptance angle. In each case the direction closest to the +y axis is the only direction from the quadrature set within the acceptance angle, and the adjoint source is defined along that direction in each of the four octants along the +y  PAGE 44 44 axis. As was seen previously, a quadrature order of 6 is perfectly adequate for this simulat ion; however, in the following simulations much higher orders of quadrature will be used for the sole purpose of obtaining a particular direction in the quadrature set. The lowest quadrature order used here, S42, has the source direction approximately equal to the acceptance angle, while the higher quadrature orders will have direction cosines smaller than the acceptance angle cosine. The group 1 adjoint flux obtained from the TITAN model of a singlehole collimator using S62 is shown in Figure 48 for an axial slice and a sagittal slice through the adjoint source location. The adjoint flux is shown on an exponential scale due to its significant dropoff away from the adjoint source. At the far left of Figure 48 the peak adjoint flux value can be seen i n the voxel containing the adjoint source. To the right of the figures and far from the adjoint source ray effects2 (i.e. the effects of not having a quadrature direction actually along the y axis) are visible. However this is seen in the air behind the torso phantom where there is no source and should therefore not contribute to the calculated response. Also, the majority of the detector response will come from the intense source in the heart, which is in the front of the torso where this effect is not visible. Table 46 gives the responses for several different quadrature orders of the angular adjoint source and the percent difference relative to the MCNP5 solution. The direction weight and the area represented at the collimator entrance by this weight are also tabulated. The results in Table 46 show that the best choice is S62 and this order of quadrature places the adjoint source direction well within the acceptance angle. The solid angle range that this direction represents ( ) then covers the di stribution of PAGE 45 45 acceptable angles. The solid angle range represented is not a direct match to the circular acceptance angle of the MCNP5 collimator; however, it seems that the TITAN methodology makes a good approximation. The relative difference in the res ponse for S62 is 0.23% and this demonstrates that the adjoint methodology of simulating a collimator in TITAN can produce accurate results when the appropriate adjoint source direction can be obtained. In Figure 49 the relative difference in each of the four orders of quadrature is plotted versus the adjoint source direction weight, which is proportional to the area subtended by the adjoint source direction. A linear regression was used to fit a line to these points and is shown in Figure 49 with a cor relation coefficient of 0.99998. This linear relationship means that the appropriate quadrature order to use could be determined from two test points. In the singlehole case above, the adjoint source was represented using only four directions in TITAN. To investigate this as a potential source of error, the ordinate splitting technique was used to split the source direction in each octant into several directions whose weights summed to weight of the original direction. These cases are presented in Table 4 7 along with their nonsplit counterparts for comparison. The S60, S62, and S68 cases used PNTN splitting of various orders. It is clear from the table that while there is a measurable change in the response with ordinate splitting, the difference is slight and one direction per octant is adequate. The remaining collimator cases will use just one direction per octant. The four and ninehole cases were considered next with one adjoint source direction per octant. Tables 48 and 49 give the response and relative differences for various orders of quadrature compared to the MCNP5 response for the four and nine PAGE 46 46 hole cases, respectively. The most accurate quadrature is S134 for the four hole collimator with a relative difference of 0.94% and S188 for th e ninehole collimator with a relative error of 0.12%. As with the singlehole case, it is clear that the adjoint methodology in TITAN can give accurate calculations of response with collimation. The area corresponding to the direction weight that is closest to the MCNP5 response is also consistently larger than the collimator hole area in MCNP5. The relative differences as a function of direction weight for the four and ninehole cases are plotted in Figure 410 and Figure 411, respectively. These pl ots also have lines from linear regression fits that show, as seen in the singlehole case, that the relationship between the relative difference in response and the direction weight is linear. This relationship makes sense when one considers that the dir ection weight is directly proportional to the solid angle represented by the direction. Given the above described linear regression fits for the single, four, and ninehole collimator cases, an attempt was made to create a model to predict the best order o f quadrature for a given acceptance angle in MCNP5. The fitted lines for the collimator cases were used to find direction weights corresponding to relative differences of zero. These ideal direction weights were then plotted versus the cosines of the c ollimator acceptance angles. There appeared to be a linear relationship between the weight/area and the direction cosine, therefore, a line was fitted to this data using a linear regression and each of the acceptance angle cosines for the collimators was then substituted into the fitted equation to get a predicted direction weight. The predicted direction weights were then compared with the direction weights of different quadrature orders to obtain the closest match. PAGE 47 47 Initially, the best choice of quadrature order was correctly predicted for the singlehole case; however, wrong quadrature orders were obtained for the other two collimator cases. The predictions were close to the best answer so it was hypothesized that if the MCNP5 data had too large of a variance the fitted lines would not be accurate enough to predict the best direction weights. To this effect, the four hole collimator MCNP5 response had a standard deviation of 4.3% and was noticeably larger than the other two cases, which had standard deviations of ~3%. Hence, for the four hole collimator case, the number of histories was increased until a standard deviation of ~3% was achieved. Figure 412 gives the plot and the fitted line from a linear regression performed using the new ideal direc tion weight for the four hole collimator. The data showed linear behavior with a correlation coefficient of 0.999998, resulting in accurate prediction of the quadrature orders that had given the smallest relative differences. These findings are summarized in Table 410. More analysis is needed with additional acceptance angles to verify that this method of predicting the best order of quadrature can work consistently, however these results are promising. The time for calculation of the single, four, and nine hole collimator cases for each code are tabulated in Table 411. All calculations were performed in parallel on the Dirac PCcluster using 16 processors. For the MCNP5 code, these results represent the time taken to achieve the uncertainty given wit h the response values in previous tables. For the TITAN code, it is clear that time becomes an important consideration as the acceptance angles are reduced and higher orders of quadrature are required. For these cases, TITAN shows a significant speedup over the MCNP5 code. PAGE 48 48 A B C D Figure 41. Normalized TITAN projection images at four angles around the patient A) Anterior B) Left Lateral C) Posterior D) Right Lateral Table 41. Comparison between myocardium pixels in group 1 anterior projectio n images with increasing order of quadrature for different number of voxels (meshes) Differencing scheme Quadrature order 64 x 64 x 64 128 x 128 x 128 Infinity norm 2 norm Infinity norm 2 norm DDZ S 4 2.36% 0.0536 2.24% 0.1410 S 6 1.12% 0.0272 1.15% 0.0621 S 8 0.48% 0.015 0.57% 0.0389 S 10 0.34% 0.0104 0.38% 0.0322 S 12 DTW S 4 1.58% 0.0404 1.71% 0.1038 S 6 0.82% 0.0209 0.87% 0.0439 S 8 0.39% 0.0121 0.38% 0.0268 S 10 0.20% 0.0063 0.21% 0.0156 S 12 PAGE 49 49 Table 42. Maximum absolute relative differ ence between TITAN flux distributions in the myocardium of the 64 x 64 x 64 model relative to the 128 x 128 x 128 model S4 S6 S8 S10 S12 DDZ 0.127% 0.168% 0.159% 0.188% 0.187% DTW 0.154% 0.191% 0.168% 0.173% 0.175% Figure 42. Infinitynorm and 2norm relati ve to 131 million photons per projection in SIMIND A B C Figure 43. SIMIND three ener gy group anterior projection images A) Group 1 B) Group 2 C) Group 3 PAGE 50 50 A B C Figure 44. TITAN three energy group anter ior projection images A) Group 1 B) Group 2 C) Group 3 A B C D E F G H Figure 45. SIMIND and TITAN normalized projection images of the heart. A) SIMIND anterior. B) SIMIND left lateral. C) SIMIND posterior. D) SIMIND right lateral. E) T ITAN anterior. F) TITAN left lateral. G) TITAN posterior. H) TITAN right lateral. Table 43. Infinity norm and 2norm for TITAN S6DDZ and S6DTW compared to SIMIND for 131 million particles per projection Projection angle Anterior Left lateral Posterior Right lateral DDZ Infinity norm 15.0% 29.0% 31.8% 42.7% 2 norm 0.445 0.674 0.893 1.06 DTW Infinity norm 19.9% 30.1% 27.6% 49.2% 2 norm 0.545 0.679 0.868 1.08 Table 44. Three group computation time comparison for four projection angles Code Time (minutes) Speedup factor SIMIND (131 million particles) 377.0 1.0 SIMIND (1 million particles) 11.9 31.8 TITAN (S 6 DDZ) 7.8 48.3 TITAN (S 6 DTW) 8.3 45.4 PAGE 51 51 Table 45. Comparison of relative differences of TITAN S6DDZ and S6DTW compared to SIMIND for 131 million particles per projection including the infinity norm, the next largest difference after the infinity norm, and the percentage of pixels with absolute relative differences greater than 10% Projection angle Anterior Left lateral Posterior Right lateral DDZ Infinity norm 15.0% 29.0% 31.8% 42.7% Next largest difference 14.6% 28.7% 27.7% 32.7% Fraction > 10% 24.2% 38.2% 59.5% 36.3% DTW Infinity norm 19.9% 30.1% 27.6% 49.2% Next largest difference 17.1% 27.4% 26.3% 29.0% Fraction > 10% 27 .3% 34.3% 60.6% 45.1% PAGE 52 52 A B Figure 46. Plot of normalized flux along a single profile of anter ior projection images created by the SIMIND and TITAN code along with the relative difference through A) column 44 and B) row 38 PAGE 53 53 A B Figure 47. Plot of normalized flux along a single profile through row 38 of SIMIND and TITAN projection images and the relative difference for A) left lateral, B) posterior, and C) right lateral angles PAGE 54 54 C Figure 47. Continued PAGE 55 55 A B Figure 48. Group 1 adjoint f lux through the heart from the TITAN model of a singlehole collimator for A) and axial slice and B) a sagittal slice PAGE 56 56 Table 46. Single hole collimator responses for TITAN compared with MCNP5 Response (counts/s) Relative difference Area (cm 2 ) Direction weight MCNP5 0.0749 2.13% 0.3068 S 42 0.1676 123.71% 0.7180 3.517E 3 S 60 0.0800 6.75% 0.3532 1.730E 3 S 62 0.0748 0.23% 0.3313 1.621E 3 S 68 0.0618 17.48% 0.2753 1.349E 3 Figure 49. Single hole collimator relative difference in TITAN respons e compared to MCNP5 as a function of direction weight (order of quadrature) Table 47. Single hole collimator responses for TITAN with PNTN ordinate splitting compared with MCNP5 Code Source directions per octant Response (counts/s) Relative difference MCNP5 0.0763 2.13% TITAN S60 No splitting 1 0.0800 6.75% P N T N Order 3 6 0.0800 6.78% P N T N Order 6 21 0.0801 6.88% TITAN S62 No splitting 1 0.0748 0.23% P N T N Order 3 6 0.0748 0.15% P N T N Order 5 15 0.0749 0.08% TITAN S 68 No splitting 1 0.0618 17.48% P N T N Order 6 21 0.0620 17.24% PAGE 57 57 Table 48. Four hole collimator responses for TITAN compared with MCNP5 Response (counts/s) Relative difference Area (cm 2 ) Direction weight MCNP5 0.0125 2.94% 0.0631 S 90 0.0287 123.83% 0.1575 7.713E 4 S 130 0.0138 7.20% 0.0756 3.703E 4 S 134 0.0130 0.94% 0.0712 3.486E 4 S 136 0.0126 2.99% 0.0691 3.384E 4 S 138 0.0122 4.78% 0.0671 3.287E 4 Table 49. Nine hole collimator responses for TITAN compared with MCNP5 Response (counts/s) Relative d ifference Area (cm 2 ) Direction weight MCNP5 0.00704 3.58% 0.0297 S 130 0.01459 107.38% 0.0756 3.703E 4 S 180 0.00767 9.05% 0.0395 1.730E 4 S 188 0.00704 0.12% 0.0362 1.621E 4 S 190 0.00690 1.94% 0.0354 1.349E 4 Figure 410. Four hole collimator relative error in TITAN response as a function of direction weight (order of quadrature) PAGE 58 58 Figure 411. Ninehole collimator relative error in TITAN response as a function of direction weight (order of quadrature) Figure 412. Optimal TITAN direction weight versus the cosine of the collimator acceptance angle in MCNP5 PAGE 59 59 Table 410. Ideal and fitted direction weights for each collimator case and the resulting recommended order of quadrature Collimator Acceptance angle cosine Ideal direction weight Fitt ed direction weight Recommended order of quadrature Single hole 0.99850 1.6222E 3 1.6221E 3 S 62 Four hole 0.99969 3.4534E 4 3.4646E 4 S 134 Nine hole 0.99985 1.7709E 4 1.7610E 4 S 188 Table 411. Calculation times for 16 processors in MCNP5 and TITAN for each collimator Collimator Code Time (min) Speedup factor Single hole MCNP5 366.5 TITAN 39.3 9.3 Four hole MCNP5 1066.8 TITAN 180.7 5.9 Nine hole MCNP5 1348.8 TITAN 390.5 3.5 PAGE 60 60 CHAPTER 5 CONCLUSIONS AND FUTURE WORK Several aspects of s ingle photon emission computed tomography (SPECT) simulation using the TITAN code were investigated. For this simulated myocardial perfusion study, low orders of quadrature were found to be adequate with increasing the order resulting in only about 2% improvement or less in projection images. The results using the diamond differencing with zero fix up (DDZ) and directional thetaweighted (DTW) differencing schemes showed no significant advantages for either method. This is attributed to the fact that mes h sizes are small enough for DDZ to be accurate. In comparing with SIMIND, the projection images created using the TITAN code were found to be in good visual agreement for energy group 1. However, significant differences ranging from 15% to 49% were found when the projection images were compared numerically. These differences are believed to be primarily the result of TITAN using cross sections generated using the CEPXS code, whereas the SIMIND code does not allow the input of any cross sections. Compar ison of the lower energy group projection images revealed discrepancies in attenuation in bone that are also thought to be the result of SIMIND not having cross section information. It was determined that the SIMIND methodology could not accurately repres ent lower energy groups that do not contain the source. A SPECT model was created in MCNP5 to provide a Monte Carlo comparison to the adjoint methodology of the TITAN code The adjoint methodology used in TITAN to simulate collimation was found to provide highly accurate results for detector response with three different collimator angles obtaining responses within 1% of the MCNP5 result. A linear relationship was seen between the direction weight used to define the PAGE 61 61 adjoint source in TITAN and the relativ e difference of the response with the MCNP5 result. A line obtained via linear regression was able to predict the best order of quadrature for the TITAN code for simulating a given MCNP5 collimator acceptance angle. The TITAN code also proved to be at least three times as fast as the MCNP5 calculation for each acceptance angle simulated. However, TITAN calculation times did increase significantly as extremely high orders of quadrature were needed to represent smaller acceptance angles. With future work, it is hoped that the TITAN code can be used to develop algorithms resulting in improved image quality. By studying and applying the results of SPECT simulations in the TITAN code, the inevitable scatter that reduces image quality in patient SPECT images c an be characterized and predicted. If a strong understanding of the scatter in an image can be obtained, this source of image degradation can be reduced. With improved image quality, disease diagnosis can be more accurate, treatment can be more effective, and doses delivered to patients can be reduced without sacrificing image quality. PAGE 62 62 APPENDIX A INPUT FILES The general input files used in this research are found below. The parameters adjusted for different cases are described. To conserve space most files had to be abbreviated however any important parameters should still be present. The NCAT input file consists mostly of default parameters. The parameters of note are the pixel and slice width, array size, Start_slice, and end_slice. These are given in Figure A 1 for the 64 x 64 x 64 voxel phantom and were adjusted to create the 128 x 128 x128 voxel phantom. PAGE 63 63 Figure A 1. NCAT input file PAGE 64 64 Figure A 1. Continued PAGE 65 65 The CEPXS input file to generate the multigroup cross sections is found in Figure A 2. Figure A 2. CEPXS input file PAGE 66 66 The TITAN input file for SPECT simulation includes a section, Section 10, where the SPECT specific parameters are specified. Some parameters that were varied for different cases include: oquad, cmdi ff, cmxfin, cmyfin, and cmzfin. Figure A 3 is the input file for the 64 x 64 x 64 voxel phantom with S4 quadrature and the DDZ differencing scheme. Figure A 3. TITAN input file for SPECT simulation PAGE 67 67 Figure A 3. Continued PAGE 68 68 Figure A 3. Continued PAGE 69 69 Figure A 3. Continued PAGE 70 70 The SIMIND code provides a program for setting up parameters and outputs a parameter file. One of these is given in Figure A 4. Figure A 4. SIMIND parameter file PAGE 71 71 The MCNP5 input file is in Figure A 5 for the four hole case To create the single and ninehole cases the number, position, and size of the collimator hole surfaces have to be changed. Figure A 5. MCNP5 input file for the four hole case PAGE 72 72 Figure A 5. Continued PAGE 73 73 Figure A 5. Continued PAGE 74 74 The TITAN input file f or the adjoint methodology follows the same form as the SPECT input file in Figure A 3 except that Section 10 is no longer used. In Figure A 6, the S62 quadrature is used with the DTW differencing scheme. The srcang parameter specifies the adjoint source in four octants along the quadrature direction closest to the +y axis. The srcdis parameter indicates the single voxel that the adjoint source is in. The splitq parameter is commented out in this figure, but would be used for ordinate splitting. Fig ure A 6. TITAN input file for adjoint methodology PAGE 75 75 Figure A 6. Continued PAGE 76 76 Figure A 6. Continued PAGE 77 77 APPENDIX B DATA PROCESSING CODE S The codes in this section were written to prepare input files, format data and analyze outputs. The imageCompare prog ram was developed to compare the binary projection images created by the SIMIND and TITAN codes. The program prompts the user for two input files and gives the infinity norm and 2norm of the second file relative to the first. PAGE 78 78 Figure B 1. Code to compare binary projection images from the TITAN and SIMIND codes PAGE 79 79 Figure B 1. Continued PAGE 80 80 Figure B 1. Continued PAGE 81 81 The response program in Figure B 2 reads in the flux distribution from an adjoint calculation in TITAN and returns the detector response. It r equires the problem source distribution. Figure B 2. Program to calculate the detector response using the adjoint flux produced from an adjoint calculation in the TITAN code PAGE 82 82 Figure B 2. Continued PAGE 83 83 Figure B 2. Continued PAGE 84 84 Figure B 2. Continued PAGE 85 85 T he phanMat program in Figure B 3 turns the NCAT material distribution into an array of material numbers that can be directly placed into the MCNP5 input. In the MCNP5 input, these material numbers represent universes in a lattice. The program creates a s imple distribution of material numbers and a more compact distribution using a shorthand notation readable by MCNP5. PAGE 86 86 Figure B 3. Program to convert NCAT material distribution into material numbers for MCNP5 input PAGE 87 87 Figure B 3. Continued PAGE 88 88 Figure B 3. Continued PAGE 89 89 The regression program in Figure B 4 was written to perform a linear regression on an input file containing the number of points followed by the x and y values. Figure B 4. Linear regression program PAGE 90 90 Figure B 4. Continued PAGE 91 91 The program in F igure B 5 reads in a 64 x 64 x 64 flux distribution and a 128 x 128 x 128 flux distribution and calculates a maximum absolute relative difference. Figure B 5. Program to compare flux distributions of two different mesh sizes PAGE 92 92 Figure B 5. Continued PAGE 93 93 L IST OF REFERENCES 1. M. H. K ALOS and P. A. WHITLOCK, Monte Carlo Methods Volume 1: Basics John Wiley and Sons, New York (1986). 2. E. E. LEWIS and W. F. MILLER, Computational Methods of Neutron Transport John Wiley & Sons, New York (1984). 3. C. YI and A. HAGHIGHAT Hybrid SN and Ray Tracing with Fictitious Quadrature for Simulation of SPECT, Proc. Int. Mtg. Mathematical Methods for Nuclear Applications (M&C 2009) Saratoga Springs, New York, May 2009, American Nuclear Society (2009) (CD ROM). 4. W. P SEGARS Dev elopment and application of the new dynamic NURBS based cardiatorso (NCAT) phantom, PhD Thesis, University of North Carolina (2001). 5. C. YI and A. HAGHIGHAT, A Three Dimensional Block Oriented Hybrid Discrete Ordinates and Characteristics Method, Nucl. Sci. Eng. 164 221 (2010). 6. C. YI, Hybrid Discrete Ordinates and Characteristics Method for Solving the Linear Boltzmann Equation, PhD Thesis, University of Florida (2007). 7. G. I. BELL and S. GLASSTONE, Nuclear Reactor Theory Robert E. Krieger Publishing Malabar, FL (1985). 8. B. G. PETROVIC and A. H AGHIGHAT Analysis of Inherent Oscillations in Multidimensional SN Solutions of the Neutron Transport Equation, Nucl. Sci. Eng. 124 31 (1996). 9. B. G. P ETROVIC and A. HAGHIGHAT A New Directional Weighted SN Differencing Scheme, Trans. Am. Nucl. Soc. 73, 195 (1995). 10. L. J. LORENCE J. E. MOREL, and G. D. VALDEZ Users Guide to CEPXS/ONELD: A One Dimensional Coupled ElectronPhoton Discrete Ordinates Code Package, Sandia National Laboratory (1989). 11. ICRP, Basic anatomical and physiological data for use in radiological protection: reference values. ICRP Publication 89, Annals of the ICRP (2001). 12. S. R. CHERRY, J. A. SORENSON, and M. E. PHELPS, P hysics in Nuclear Medicine, Saunders, Philadelphia, PA (2003). 13. G. LONGONI and A. HAGHIGHAT, Development of New Quadrature Sets with the Ordinate Splitting Technique, Proc. Int. Mtg. Mathematical Methods for Nuclear Applications (M&C 2001), Salt Lake City, Utah, September 2001, American Nuclear Society (2001) (CD ROM ). PAGE 94 94 14. M. LJUNGBERG, S. STRAND, and M. A. KING, The SIMIND Monte Carlo program: A Monte Carlo Calculation in Nuclear Medicine, Applications in Diagnostic Imaging 11, 145 (1998). 15. X 5 MONTE CARLO TEAM MCNP A General Monte Carlo for Neutron and Photon Transp ort, Version 5, Los Alamos National Laboratory (2003). PAGE 95 95 BIOGRAPHICAL SKETCH Katherine Royston was born in West Palm Beach, Florida in 1984. She graduated from Wellington High School in 2002 and started her undergraduate career at the University of Florid a in the fall of that same year. In the summer of 2005, she completed the Research Experiences for Undergraduates Program at the University of Notre Dame. She completed her Bachelor of Science in physics with a minor in business administration at the Uni versity of Florida in May of 2006. She then worked for a University of Florida startup company, Viewray, until August of 2007 when she was accepted to graduate school at the University of Florida. 