Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The development of multiple line transmission sources for SPECT Sitek, Arkadiusz 1998

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

Item Metadata

Download

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

Full Text

THE DEVELOPMENT OF MULTIPLE LINE TRANSMISSION SOURCES FOR SPECT By Arkadiusz Sitek M.Sc , Warsaw University, 1994 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF PHYSICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March 1998 © Arkadiusz Sitek, 1998 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of Br i t i sh Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Physics The University of Br i t i sh Columbia 6224 Agricul tural Road Vancouver, B . C . , Canada V 6 T 1Z1 Date: Abstract In this work a new kind of transmission source for S P E C T imaging, based on multiple line sources, is presented. The source may be used for either simultaneous or sequential S P E C T transmission imaging to measure patient specific attenuation maps. Two different, new transmission sources based on the multiple line concept, the col-limated line sources and the multiple line array, were proposed and tested using computer simulations. Based on simulations, prototypes of the transmission sources were built, and series of experiments were performed. Computer simulations and experiments showed that multiple line transmission sources provided accurate attenuation maps. Compared with other transmission sources, the multiple line transmission sources do not require complicated hardware or expensive maintenance, but st i l l give good quality attenuation maps. The multiple line array transmission source has been adopted for commercial use by Siemens Medical Systems with their ecam S P E C T camera. 11 Table of Contents A b s t r a c t i i Tab l e o f C o n t e n t s i i i L i s t o f F i g u r e s v i i i L i s t o f Tab l e s x i x A c k n o w l e d g m e n t s x x i 1 I n t r o d u c t i o n 1 1.1 A i m of This Work 1 1.2 Structure of the Thesis 1 1.3 Diagnostic Medical Imaging 3 1.4 Nuclear Medicine Imaging 3 1.4.1 Planar and S P E C T Imaging 4 1.5 Anger Camera 6 1.5.1 Collimator 6 1.5.2 Scintillation Crystal 7 1.5.3 Photomultiplier Tubes and Electronics 7 1.5.4 Computer 8 1.6 Data Reconstruction in S P E C T 12 1.6.1 Fil tered Backprojection 12 1.6.2 Iterative Methods 14 1.7 Basic Physics of Photon Attenuation 16 i i i 1.7.1 Linear Attenuation Coefficient 17 1.8 Effect of Attenuation on Reconstructed S P E C T Images 18 1.9 Correction for Attenuation 20 1.10 Importance of Attenuation Correction 21 2 Di f f e ren t M e t h o d s of M e a s u r i n g A t t e n u a t i o n M a p s 23 2.1 Computed Tomography ( C T ) - X-rays 24 2.2 Sheet Source 25 2.3 Scanning Line Transmission Source 26 2.4 Fan and Cone Beam Geometry Transmission Source 28 2.5 Line Sources perpendicular to the axis of rotation 29 2.6 Conclusion 30 3 M u l t i p l e L i n e T r a n s m i s s i o n S o u r c e G e o m e t r y 31 3.1 Collimated Line Sources (CLS) 33 3.2 Mul t ip le Line Array 35 3.3 Summary 37 4 M e t h o d s o f R e c o n s t r u c t i o n o f the A t t e n u a t i o n M a p s 38 4.1 Collimated Line Sources - Attenuation Map Reconstruction 38 4.1.1 Data Preprocessing - Rebinning 39 4.1.2 Strip Geometry - Fil tered Backprojection 40 4.1.3 Exact Geometry - Iterative Approach 43 4.2 Mult ip le Line Ar ray -Map Reconstruction 50 4.2.1 Data Preprocessing - Missing Counts Correction Algor i thm . . . . 50 4.2.2 Fil tered Backprojection 51 4.2.3 Iterative Methods 51 iv 4.3 Summary and Conclusions 54 5 S i m u l a t i o n E x p e r i m e n t s 56 5.1 Goals of the Simulations 56 5.2 Methods 57 5.2.1 Analyt ic Phantom 57 5.2.2 Data Generation 58 5.2.3 Attenuation M a p Reconstruction 59 5.2.4 Error Analysis 59 5.3 Results - Simulations of C L S 60 5.3.1 Geometry of the Simulated C L S 60 5.3.2 Reconstruction Methods for C L S 62 5.3.3 Number of Divisions 73 5.3.4 Beam Profile, Thickness of the Lines 74 5.3.5 Spatial Resolution - C L S 77 5.4 Results - Simulations of M L A 79 5.4.1 Geometry of the Simulated M L A Source 79 5.4.2 Reconstruction Methods for M L A 81 5.4.3 M L A Beam Profile 91 5.4.4 Spatial Resolution - M L A 92 5.5 Summary and Conclusions 93 6 E x p e r i m e n t s - G e n e r a t i o n o f A t t e n u a t i o n M a p s 96 6.1 Goal of the Experiments 96 6.2 Methods 96 6.2.1 Prototype of the C L S 96 6.2.2 Design of the M L A 100 v 6.2.3 Phantoms Used in Experiments 103 6.2.4 Data Acquisi t ion and Reconstruction 103 6.3 Results 104 6.3.1 Experimental Results - C L S Transmission Scans 104 6.3.2 Experimental Results - M L A Transmission Scans 107 7 E x a m p l e o f the A p p l i c a t i o n o f the A t t e n u a t i o n M a p s 116 7.1 Methods 117 7.1.1 Simultaneous Transmission/Emission Data Acquis i t ion 117 7.1.2 Cross Talk Correction 117 7.1.3 Emission Reconstruction 119 7.2 Results 120 7.2.1 Estimation of the Ratio Parameter k 120 7.2.2 Effect of Attenuation Correction 120 7.3 Discussion and Conclusion 121 8 C o n c l u s i o n s 125 8.1 Summary of the Work 125 8.2 Advantages of Mul t ip le Lines Transmission Sources 127 8.3 Cl in ica l Application 128 8.4 F ina l Word 132 B i b l i o g r a p h y 133 A G e n e r a t i o n of the S y s t e m M a t r i x for the C o l l i m a t e d L i n e S o u r c e 141 A . l Calculation of aij 142 A.2 Compression of Sparse Ma t r i x w 142 A.3 Code for Generation of the System Mat r i x for C L S 144 v i B G e n e r a t i o n o f the S y s t e m M a t r i x for the M u l t i p l e L i n e A r r a y 152 B . l Calculation of the Ma t r i x Elements 152 B.2 Ordering of the Elements in the Series 153 B. 3 Code for the Generation of the System Mat r ix for M L A 154 C A n a l y t i c a l T r a n s m i s s i o n P r o j e c t i o n G e n e r a t o r ( A T P G ) 160 C. l Collimator Acceptance Function, C A F 161 C.2 Source Coll imation Function, S C F 163 C.3 Attenuation Length, AttLength 163 C.4 Normalization constant, A . Poisson Noise 164 D V e r i f i c a t i o n o f E q . 4.10 166 E I m p l e m e n t a t i o n o f A t t e n u a t e d P r o j e c t o r / B a c k p r o j e c t o r 170 v i i List of Figures 1.1 Planar and S P E C T imaging, (a) Information from objects lying on a line parallel to the z-axis is overlapped in the final image, (b) S P E C T imaging can distinguish objects overlapping in the planar study by reconstruction of a 3D image from views (projections) taken at different angles 5 1.2 Parallel hole collimator. View from the top and from the side, (a) photon coming at 0° is absorbed by the septum (b) photon coming 0° is detected, (c) photon is absorbed in the collimator; its incident angle was more the acceptance angle, (d) photon is detected after septal penetration, (e) photon incoming with an angle smaller than the acceptance angle is detected, (f) photon with an angle smaller than the acceptance angle is absorbed in the septum. A typical acceptance angle is about 3-4° 10 1.3 Schematic drawing of the Anger camera and of its main components. . . . 11 1.4 Schematically described backprojection. (a) slice of the object and 3 sample projections, (b) simple backprojection results in star artifacts, (c) filtered backprojection removes most of the star artifacts around the object. 13 1.5 Definition of system matrix w. Element W{j is the area of overlap of pixel i with the parallel strip coming from bin j 15 1.6 (a) Schematically described narrow beam geometry, (b) definition of the linear attenuation coefficient 18 1.7 Fraction of radiation escaping from a water cylinder 30 cm in diameter containing a point source vs. position of the point source on the diameter. The fraction is an average calculated over 360° 19 v i i i 1.8 Reconstructed heart phantom with its profile. The bright area is the uni-form heart. The reconstructed heart has darker and brighter areas due to attenuation, the part of the heart deeper in the phantom is darker. The reconstructed warm background also appears to be not uniform due to attenuation 20 1.9 The effect of attenuation on the number of counts in projections for (a) source without attenuation and (b) with attenuation 21 1.10 Attenuation artifacts in a reconstructed Jaszczak phantom. The Jaszczak phantom is a water cylinder, 22 cm in diameter filled with uniform back-ground activity. There are 1 cold and 2 hot spheres inside the phantom. The cold spot is an air bubble; it is not seen at al l due to severe attenu-ation artifacts. The shapes of the hot spheres are deformed, and the warm background of the phantom is not uniform due to streaks going from the hot spheres 21 1.11 Attenuation correction reconstruction of the Jaszczak phantom, (a) no correction (b) with correction (c) attenuation map used for the correction of this reconstruction 22 2.1 Configuration of the sheet source (a) without collimator on the source and (b) with collimator. For the sheet source without collimator (a) a lot of scatter radiation can be detected changing the resulting attenuation coeffi-cient (broad beam effect). This scatter is mostly eliminated by placing a collimator in front of the source (b) 25 ix 2.2 The scanning line source. During each S P E C T projection the whole area of the detector is scanned. A narrow area on the detector opposite the line source (shaded on the figure) is electronically windowed to detect only the transmission data; the rest of the detector acquires the emission data. . . 27 2.3 The fan beam geometry of the transmission source. A single line source can be positioned in the centre (left), or to reduce truncation artifacts, off centre (right). The line source has to be located at the focal position of the converging collimator 29 2.4 The line sources perpendicular to the axis of rotation 30 3.1 Geometry of multiple line sources. Collimated multiple line sources are placed parallel to the axis of rotation of a S P E C T camera, (a) view per-pendicular to the axis of rotation (b) view along the axis of rotation. . . . 32 3.2 (a) The geometry of the collimated line sources (CLS) transmission source. The simulated (b) and experimental (c) profiles of the blank scan photon beam 34 3.3 (a) The geometry of the multiple line array ( M L A ) transmission source (only 8 line sources are shown). The simulated (b) and experimental (c) profiles of the blank scan photon beam 36 4.1 Pre-reconstruction rebinning. In this example each fan beam was divided into 3 smaller fan beams (divisions) 40 x 4.2 Strip geometry in F B P . (a) The first two images show the two opposite projections which wi l l be combined into one data set. The shaded strips represent our approximation of the fan beams (solid lines). The third image is the combination of these two opposite projections. The dot in the middle represent the axis of rotation. This procedure can be performed on data where the ini t ia l beam is divided into any number of subdivisions, (b) Rebinning procedure where the beam is divided into 2 divisions 42 4.3 Method of calculating system matrix w for multiple fan beams 44 4.4 A schematic presentation of the matrix Qj,fc. Horizontal and vertical neigh-bours contribute to the penalty with weight 1, diagonal with 0.71, and 0 otherwise 50 4.5 The geometry of the M L A source. Schematically shown are one pixel of the reconstruction area, i, and one projection bin, j 53 5.1 Central slice of the analytic phantom. The two square areas are the regions of interest (ROIs) used later in the analysis 57 5.2 Simulated profiles of a one line source, (a) no collimation on the source (b) low collimation (c) high collimation 61 5.3 Accuracy and precision of the ROIs plotted vs. cutoff frequency. The top row corresponds to low activity in the C L S source, the middle, and the bottom rows to middle and high activities respectively 63 5.4 Central slices of the attenuation maps from the medium C L S source activity reconstructed using F B P with the Butterworth filter (order 5), with the cutoff frequencies of 0.45, 0.65, 0.85, and 1.00 respectively for images (a), (b), (c) and (d) 64 x i 5.5 The local accuracy and % uncertainty, plotted vs. number of iterations for the maximum likelihood algorithm using medium activity in the source. The 4 lines correspond to M L and P M L with different penalty parameter reconstructions. The value of penalty parameter is given (in the brackets) 5.6 The value of log-likelihood plotted vs. number of iteration for maximum likelihood reconstruction. The figure corresponds to a total activity in the C L S source equal to 4.8 G B q . The 4 lines correspond to M L and P M L with different penalty parameter reconstructions. NB: The plotted function is not an M L objective function but only the value of the log-likelihood. The value of penalty parameter is given (in the brackets) in the legend 66 5.7 Central slices of the attenuation maps reconstructed by the maximum like-lihood method with (a) 20 iterations of non-penalized maximum likelihood, (b) 20 iterations of P M L with penalty parameter a = 0.4, (c) 20 iterations of P M L with a = 1, and (d) 20 iterations of P M L with a = 2. The total activity in the C L S source was 4.8 G B q 67 5.8 Local accuracy and % uncertainty of ROI-1 plotted vs. number of itera-tions. Rows a), b), and c) correspond to the total activities in the C L S source equal to 1.4, 4.8, and 11.8 G B q respectively. The value of penalty parameter is given (in the brackets) in the legend 69 5.9 Accuracy and % uncertainty of ROI-2 plotted vs. iteration number. Rows a), b), and c) correspond to different activities in the C L S source equal to 1.4, 4.8, and 11.8 G B q respectively. The value of penalty parameter is given (in the brackets) in the legend 70 in the legend 65 xn 5.10 Central slices of the attenuation maps reconstructed by the least squares method with a total activity in the C L S source of 4.8 G B q . (a) 10 iterations of L S , (b) 10 iterations of W L S , (c) 20 iterations of P W L S with the penalty parameter a = 5, and (d) 20 iterations of P W L S with penalty parameter a = 10 71 5.11 The attenuation maps obtained from data rebinned into different number of divisions and reconstructed by F B P for C L S medium activity. The number at each figure indicates the number of divisions used in data processing. . 74 5.12 Profiles through the source activity on the detector, (a) 1 m m line sources with uniform activity distribution across the sources, (b) 1 m m and (c) 3 m m diameter line sources with non-uniform activity distribution 75 5.13 Attenuation maps acquired using C L S with different beam profiles: (a) exponential with 1 m m , (b) uniform with 1 m m , and (c) exponential with 3 m m in diameter line sources 76 5.14 The reconstructed attenuation map of the set of 27 bars. The image is 256x256 pixels with pixel size of 0.24 cm (the area shown is 61.4x61.4 cm). 77 5.15 Ac t iv i ty distribution for the M L A . (a) the distances between the lines are shown in cm. Comparison of the beam profiles for a large 60 cm water cylinder obtained by a simple calculation (b) and using simulation of a blank projection by A T P G is presented in (c). The calculated line is produced by; Profile = Norm • exp(0.15 • d) where Norm is a factor normalizing Profile to the A T P G simulation 80 5.16 Accuracy and precision of ROIs in the reconstructed image plotted vs.-cutoff frequency. The top row corresponds to low activity in the source, the middle one to medium, and the bottom one to high activity in the M L A system 82 x i i i 5.17 The central slices of the attenuation maps of the phantom reconstructed by F B P with the Butterworth filter (order 5) and different cutoff frequencies equal to 0.45, 0.65, 0.85, and 1.00 for medium activity of M L A system. . 83 5.18 Accuracy and precision plotted vs. iteration number for the maximum like-lihood algorithm, for a total activity in the M L A source equal to 2.8 G B q . The value of penalty parameter is given in the legend 84 5.19 The value of log-likelihood vs. the number of iterations for maximum like-lihood reconstruction. The graph is for a total activity in the M L A source equal to 2.8 G B q . The lines correspond to different penalty parameters given in the legend. NB: The plotted function is not an M L objective function but the value of the log-likelihood 85 5.20 Central slices of the attenuation maps of the phantom reconstructed by the maximum likelihood method with (a) 20 iterations of non-penalized ( M L ) (b) 20 iterations of P M L with penalty parameter a = 0.4, (c) 20 iterations of P M L with a = 1, and (d) 20 iterations of P M L with a = 2 for medium activity of M L A source 86 5.21 Accuracy and precision for ROI-1 vs. the number of iterations. Rows (a), (b), and (c) correspond to activities in the M L A source of 1.4 G B q , 2.8 G B q , and 9.3 G B q , respectively. The value of penalty parameter is given (in the brackets) in the legend 88 5.22 Accuracy and precision for ROI-2 vs. the number of iterations. Rows a), b), and c) correspond to activities in the M L A source of 1.4 G B q , 2.8 G B q , and 9.3 G B q , respectively. The value of penalty parameter is given (in the brackets) in the legend 89 x i v 5.23 Central slices of the attenuation maps of the phantom reconstructed by the least squares method with a) 10 iterations of least squares, b) 10 iterations of W L S , c) 20 iterations of P W L S with penalty parameter a = 5, and d) 20 iterations of P W L S with a — 10 for medium activity of the M L A source. 90 5.24 Reconstructed image of 27 bars, each 1 cm in diameter. The image contains 128x128 pixels of 0.48 cm. Acquis i t ion was over 180°, from 0° to 180°. . 93 6.1 Picture of the prototype of C L S . The Sophy detector with 2 heads is shown. The transmission source containing 6 line sources is mounted on one of the heads (the upper one); the other head is used for data acquisition. Between the heads, the Jaszczak Phantom is being scanned 98 6.2 The design of line source for the C L S . A detailed scheme of the design is shown in (a) and a photograph of the line source prototype in (b) 99 6.3 Experimental setup for the M L A . The single head Siemens Diacam camera is used. The picture shows data acquisition for a patient heart perfusion scan. The M L A transmission source is attached to the gantry by a steel bar. Opposite the M L A , behind the patient, the detector head can be seen. 101 6.4 Coll imation of the M L A source. A block of a standard high sensitivity parallel hole collimator (a) is placed on top of the aluminum, lead shielded box containing line sources (b) (only 5 are drawn) 101 6.5 Prototype of the M L A . The line sources are shown on the bottom of the aluminum box. The source collimator used is shown on top of the box. . . 102 xv 6.6 Slices from the reconstructed attenuation maps of the lung phantom ob-tained by the CLS prototype. The rows correspond to the bolts on the lid of the phantom, a slice from the lung region, and a slice from the uniform re-gion, respectively. In each row the columns correspond to the 20 iterations of PWLS reconstructions with different values of the penalty parameter (a = 2, 5, and 10). The ROIs from which the values of [i were calculated are marked in the first column 106 6.7 Position of metal bars in experiment with M L A system 107 6.8 The four sinograms of log-ratio for 3 aluminum bars. Each sinogram shows 64 projections of 180° data starting at 0°,90°,180° and 270° (see text). . . 108 6.9 Reconstructed metal bars with PWLS. The positions 0° and 180° corres-pond to the first and the last projection 109 6.10 The relation of the maximum [i reconstructed for a given bar vs. the average distance of that bar from the detector for the M L A . The actual value of fj, was different for experiment and simulation I l l 6.11 Slices from the reconstructed attenuation maps of the phantoms. The rows correspond to the experiments with the Jaszczak, lung phantom, and lung phantom with bags. In each row the columns are the PWLS reconstructions with different values for the penalty parameter, a = 0.5, 2.0, and 5.0. The white square in column 1 designates the region in which n is determined. 113 6.12 Slices of the attenuation maps of patients taken from the lung region. Rows correspond to different patients: small, medium, and large. Columns cor-respond to different values of the penalty parameter in PWLS reconstruc-tion, a=0.5, 2.0, and 5.0 115 xvi 7.1 Schematic of the energy spectrum of the simultaneous transmission/emission acquisition. Spectrum of photons scattered in N a l is not shown. The num-ber of these photons is small compared to number photons scattered in the object 118 7.2 Energy window positions and width for the simultaneous transmission/emission ( G d / T c ) M L A transmission source 119 7.3 Example slices of the lung phantom schematicaly presented in (a) recon-structed with no attenuation correction (a) and with attenuation correction (b). Reconstruction was performed by O S E M 122 7.4 Example slices of emission reconstruction with O S E M of the lung phantom, reoriented along the short axis of the heart, and a profile, (a) with no attenuation correction and (b) with attenuation correction 123 7.5 Emission reconstructions with M L E M of the lung phantom with bags, (a) no attenuation correction and (b) attenuation correction. The contrast of both images was adjusted to show more clearly the streak artifacts. . . . 124 8.1 Slice of emission image of the heart. Reconstruction was performed by O S E M (a) without and (b) with attenuation correction. The presented image is reoriented along short axis of the heart 130 8.2 Prototype for the new Siemens 2-head ecam camera with 2 "wings" con-taining two M L A transmission sources attached to the gantry by a single mounting bar 131 A.3 Mathematical formulation of the calculation of a;j 142 A.4 A l l the possibilities of the relative position of two parallel lines to the pixel . Each has to be treated differently because the area of the pixel between them has to be calculated from different formulas 143 xv i i A . 5 Scheme of matrix WitJ; compression 143 B . 6 The ordering of the elements in the series in the system matrix for M L A . The closer the element is to the detector, the earlier in the series it is positioned after the ordering 153 C. 7 Definition of the collimator acceptance function ( C A F ) . Shaded is the area "seen" by photons incoming with angle 9 161 C.8 The collimator acceptance function (probability that photon wi l l be detec-ted) vs. the incidence angle 6 163 C.9 The S C F for the C L S 164 D.10 The approximations made during the binning of I D projection data, (a) Ex -act projection uses continuous values acquired in transmission and blank scans (b) Data is acquired in bins; the value of the bin is calculated as an average over the bin from the continuous distribution, (c) further rebinning of (b) into new bins. Method (c) was used in this thesis for C L S system. The M L A data was calculated as described in (b) 168 D . l l Error in least squares reconstruction of the simulated projection data binned with different methods 169 xvin List of Tables 1.1 Some examples of applications of nuclear medicine imaging 4 5.1 Comparison of the different methods of the reconstruction for C L S . . . . 72 5.2 Accuracy in the ROIs for F B P reconstructed maps for low and medium activity in C L S source for different number of divisions 73 5.3 The values of accuracy and precision (expressed in % uncertainty) in the ROIs and the global accuracy for different beam profiles, reconstructed by F B P with total medium activity in the C L S system 76 5.4 Comparison of the different methods of reconstruction for M L A with 2.8 GBq total activity in the source 91 5.5 The value of accuracy and precision in the ROIs for different beam profiles, reconstructed by F B P with the Butterworth filter (order 5, cutoff 0.4) for the M L A system. 92 6.1 Values of \x [ cm - 1 ] and their standard deviations found in R O I for the lung phantom scanned using transmission C L S system. Data was reconstructed with 20 iterations of P W L S with different values of the penalty parameter a. 105 6.2 Values of p in [ cm - 1 ] and their standard deviations found in ROIs for phantoms scanned using the transmission M L A system. The reconstruc-tion was done by 20 iterations of P W L S with different values of the penalty parameter a 112 6.3 Values of n [ cm - 1 ] and their standard deviations found in ROIs for patients scanned using transmission M L A system. Data was reconstructed with P W L S with different values for the smoothing parameter a 114 xix 7.1 Values of the ratio k of scatter counts between the transmission and scatter energy windows. The standard deviation, a, of k between the projections is also presented 121 C . l Parameters of the collimators. Data for the Sophy D S T Camera 162 xx A cknowle dgment s Wdziecznosc Sa. ludzie ktorzy maj§, serca czyste. S§ wspomnienia, perelki czasu. Rzeczy czasami, te zbyt oczywiste, Umkn§, naszym ubogim slowom. I nie mowie, nie p isze . . . pamietam. W i t h this short poem, I would like to thank my wonderful supervisor, Dr . A n n a Celler, for her help and guidance throughout this work. I have enjoyed working with her from the first to the last day. It was not only because she was an excellent scientist, and her help was tremendous, but also because of her interpersonal skills; she was and she is a great person. I also thank Dr . Ronald Harrop for great help with putting my work together. Dr . Harrop, Dr . Christine Dykstra and Dr . Alex M a c K a y put huge effort correcting my English, my special thanks to them. I must also express my gratitude to all who help me with many different aspects of my work: Dr . Dave Axen for help with monte carlo, Siemens Medica l Systems and T R I -U M F for building transmission devices, technologists in Nuclear Medicine Department for preparing radioactive sources, and all the doctors at V H H S C for support. I would like to thank all my friends from M I R G for great work environment. There are people who are my friends and who were al l the time with me, from them I learned a lot while my time in school, my great thanks to them: to my brother Robert, A T , and Artur . M y final thanks are for my parents. I do not think any words express my gratitude. I just say: thank you. xx i Chapter 1 Introduction 1.1 Aim of This Work Single photon emission computed tomography ( S P E C T ) is a diagnostic medical ima-ging technique based on the detection of photons emitted by radioactive tracers from inside the body of the patient. Photon attenuation within the patient is one of the major effects which reduces the quantitative accuracy of an image acquired by S P E C T . Cor-rection for photon attenuation can greatly improve the image accuracy, but it requires knowledge of the patient specific distribution of attenuating medium represented by an attenuation map. Attenuation maps can also be used for accurate corrections for scatter, another factor degrading image quality in S P E C T . The objective of this work is to develop a method for obtaining attenuation maps using a new type of transmission source and to demonstrate that this method gives accurate, good resolution attenuation maps which are free of artifacts. 1.2 Structure of the Thesis This thesis consists of eight chapters. The first chapter presents the basic ideas of 1 Chapter 1. Introduction 2 nuclear imaging; it explains the phenomenon of photon attenuation and shows how this affects a S P E C T image. There are several types of transmission sources which have been developed and used. They are described in Chapter 2 where their advantages and disadvantages are discussed. Chapter 3 describes our new type of transmission source. There are two transmission source systems presented there, both based on a usage of multiple line sources. The rationale behind the new transmission source design along with a presentation of the advantages of our transmission sources over other solutions is presented there. Methods for the reconstruction of attenuation maps from the acquired information by transmission source are essential for this work to be complete. Chapter 4 presents the methods used later in the thesis. In Chapter 5, data acquisition and reconstruction methods are tested using computer simulations for the new transmission sources. The geometry of the sources is optimized. Knowledge of the optimized geometry was then used for designing the prototypes of the transmission sources. This design is described in Chapter 6 where experiments with the prototypes are also presented, and experimentally measured attenuation maps are shown. A n example of the usage of attenuation maps for attenuation correction is described in Chapter 7. Attenuation maps were acquired in simultaneous transmission/emission scans with our multiple line source system. This chapter shows that the attenuation correction based on our attenuation maps removes artifacts caused by attenuation. Finally, a summary is presented in Chapter 8 in which some future cl inical application are also discussed. In a series of Appendices at the end of the thesis, some of the computationally chal-lenging tasks of the thesis are described, and validation is provided for a mathematical approximation made in Chapter 4. Chapter 1. Introduction 3 1.3 Diagnostic Medical Imaging In 1895 W i l l i a m Roentgen discovered X-rays, and this year may also be considered as the beginning of diagnostic medical imaging. The ability to look at human organs without the use of a scalpel has become a very useful tool of modern medical diagnosis. A t the present time almost all the hospitals around the world are equipped with medical imaging systems. There are different techniques used for medical imaging, and they use different phys-ical phenomena to obtain information about patient anatomy and/or physiology. These techniques are ultrasonography ( U S G ) , computed tomography ( C T ) , plain film x-ray, magneto-encephalography ( M E G ) , magnetic resonance imaging ( M R I ) , and the tech-niques of nuclear medicine: positron emission tomography ( P E T ) , planar gamma ima-ging, and single photon emission computed tomography ( S P E C T ) . 1.4 Nuclear Medicine Imaging Nuclear medicine imaging investigates the localization of radiopharmaceuticals inside the body of a patient by the detection of photon radiation emitted by a radionuclide from inside of the patient. The radionuclide is attached to a radiopharmaceutical which serves as a carrier for the radioactivity. Radiopharmaceuticals are administered to the patient by injection, orally, or by inhalation. The important characteristic of nuclear medicine imaging is that physiological and metabolic aspects of a patient can be studied by imaging the distribution of the "carried" radionuclides (tracers). There is a vast variety of radiopharmaceuticals which can be used in nuclear medicine, some of which are listed in Table 1.1. There are two types of nuclear medicine imaging: planar, where only a 2D image Chapter 1. Introduction 4 of radioactivity distribution inside the patient is acquired, and single photon emission computed tomography ( S P E C T ) which can create 3D images. Table 1.1: Some examples of applications of nuclear medicine imaging. Organ Radi ophar maceuti cal Energy [keV] Some Medica l Applications Heart 9 9 m T c - S e s t a m i b i 143 perfusion studies Skeleton 9 9 m T c - M D P 143 metastatic lesions, bone tumor, fracture, osteomyelitis Kidney 1 2 3 I-Orthoiodohippurate 159 functional status, arterial stenosis Lung 1 3 3 X e gas 81 regional ventilation Thyroid 1 3 1 I - S o d i u m Iodine 346 thyroid function 1.4.1 Planar and SPECT Imaging A planar image is acquired by a detector positioned close to the patient. During the scan, the detector does not move (Fig. 1.1(a)). In a planar study, al l photons originating in organs lying on the same line perpendicular to the detector contribute to the same area in the final image since the acquisition of the data sums all the photons coming from the sources along this line. It makes the diagnosis process more difficult because of uncertainty of the depth localization of the tracer (Fig. 1.1(a)). S P E C T acquires information about the tracer distribution from multiple views by rotating the detector around the patient(Fig. 1.1(b)). Using this information in the re-construction process, the 3D tracer distribution can be found. Usually in S P E C T , the detector stops at 60 to 120 positions around patient and at each position takes a planar image. Such a planar image wi l l throughout this thesis be referred to as a projection. Chapter 1. Introduction 5 Figure 1.1: Planar and S P E C T imaging, (a) Information from objects lying on a line parallel to the z-axis is overlapped in the final image, (b) S P E C T imaging can distinguish objects overlapping in the planar study by reconstruction of a 3D image from views (projections) taken at different angles. Chapter 1. Introduction 6 1.5 Anger Camera The Anger camera is the basic device used for detection of radiation in nuclear medicine [1, 2]. The important part of the camera is the scintillation crystal. A photon striking the crystal creates light (scintillation light) which is measured by photomultiplier tubes. The electronics behind the tubes determines the position of each event on the surface of the crystal, as well as the amount of scintillation light. The latter is proportional to the energy of the photon. The main components of the Anger camera are: collimator, crystal, photomultiplier tubes, electronics, and computer. 1.5.1 Collimator A collimator placed in front of the crystal ensures that radiation detected in the camera is coming from a known direction. This knowledge is crucial for the image reconstruction process in S P E C T and for the interpretation of the image in planar imaging. The collimator is basically a thick lead sheet with many small hexagonal holes in it. This device cuts the efficiency of photon detection by several orders of magnitude (approximately 104) leading to poor statistics in acquired images. There are several types of collimators. The most commonly used is a parallel hole collimator. It allows detection of photons coming only from the direction perpendicular to the detector surface. Other types of collimators include converging, diverging, and pin-hole. The parallel hole collimator which wi l l be used throughout this thesis ideally accepts only photons coming with an incident angle of 0°. The incident angle is the angle between the photon direction and line normal to the detector. The finite size of the holes also allows photons coming at some small angle, less than the acceptance angle of a collimator, to be detected (Fig . 1.2). The acceptance angle of a collimator is defined as the maximum Chapter 1. Introduction 7 incident angle for a detected photon passing through a collimator hole. A photon with incident angle larger than the acceptance angle can be detected, but only if it penetrates the lead wall of the hole (septum); this effect is known as septal penetration (Fig.1.2). 1.5.2 Scintillation Crystal The scintillation crystal is a single piece of sodium iodine N a l ( T l ) , activated with thallium (TI). Pure N a l is a poor scintillator at room temperature, but the addition of up to 1% of impuri ty (TI) substantially improves its performance. A gamma photon interacting with atoms of the crystal loses its energy, producing a large number of low energy electrons; these then recombine with the atoms giving off photons (visible light) which are detected by the photomultiplier tubes. The advantage of using N a l as a scintillator in nuclear medicine is its relatively low production cost. The size of the crystal is generally about 40x50 cm and its thickness usually about 0.95 cm. The thickness of the crystal depends on the energy of the photons for which the camera is designed. For the detection of low energy photons, the crystal can be thinner (about 9 mm) without loss of efficiency. A thinner crystal has better spatial resolution, but for higher energy radiation, it wi l l be less efficient because the probability that the photon wi l l go through the crystal undetected increases. Choice of crystal thickness is a trade off between spatial resolution and efficiency. 1.5.3 Photomultiplier Tubes and Electronics Photomultiplier tubes ( P M T s ) convert the light created in the scintillation crystal into an electric signal and amplify it so that it can be further processed by the electronic circuitry of the camera. There are usually around 50-100 P M T s in the cameras used in nuclear medicine. The position of the event on the detector surface is determined Chapter 1. Introduction 8 by the position logic circuits. The pulse-height analyzer determines the energy of the ini t ial photon by measuring the light output induced by the photon in the crystal. The information about the position of the event ( X and Y coordinates) and the energy of the photon is passed on to the computer. Because the energy of the gamma-ray is measured, it is possible to use two radioactive nuclides simultaneously and distinguish between them. Thus, 9 9 m T c (140 keV) can be used for emission studies while using 1 5 3 G d (100 keV) to measure attenuation maps of the patient. A n important feature of the camera is that it can process only one event at a time because the position logic circuits do not work correctly if two or more events occur con-currently. Because the information about energy passed to the computer is proportional to the light output of the scintillator, if two or more photons hit the detector at the same time the Anger camera determines them to be one event, but with energy equal to the sum of their energies (pile-up effect). Modern digital Anger cameras are able to process more then one event at a time using multi-zone triggering and local centroid weighting, which divide electronically the detector into several independent detectors. This reduces pile-up effects. Each camera has a certain dead time, which is the time needed for processing an event. For the Anger camera this time is on the order of microseconds. 1.5.4 Computer The information from the detector is collected, digitized and passed to a computer. The computer is needed for processing the information, so that it can be displayed on the screen. Usually the surface of the detector is divided into a square matrix from 64x64 to 1024x1024 pixels and each event is assigned to the appropriate pixel according to its x and y positions. Such digitized information can then be processed (image reconstruction, filtering, reorientation statistical analysis etc.) and displayed. Computers are also used Chapter 1. Introduction for transmission of the information, and for storage. Chapter 1. Introduction 10 A Acceptance angle i (a) (b) (c) (d) (e)(f) Figure 1.2: Parallel hole collimator. View from the top and from the side, (a) photon coming at 0° is absorbed by the septum (b) photon coming 0° is detected, (c) photon is absorbed in the collimator; its incident angle was more the acceptance angle, (d) photon is detected after septal penetration, (e) photon incoming with an angle smaller than the acceptance angle is detected, (f) photon with an angle smaller than the acceptance angle is absorbed in the septum. A typical acceptance angle is about 3-4°. Chapter 1. Introduction 1 1 Computer Electronics Detector cover \ A 1 PM Tubes ] - Light guide _- Nal crystal ^~ Collimator Radioactive source Figure 1.3: Schematic drawing of the Anger camera and of its main components. Chapter 1. Introduction 12 1.6 Data Reconstruction in SPECT The data acquired during a SPECT study constitutes a set of projections taken at different angles around the object. The goal of the reconstruction process is to recover a 3D image of the source distribution from acquired 2D projections. Reconstruction of a 3D image can be performed by either a 2D or a 3D reconstruction. In a 2D reconstruction the projections are divided into sets of slices, and each slice is reconstructed independently. In a 3D reconstruction the whole image of an object is directly reconstructed from the 2D projections. In this thesis only 2D reconstructions are considered. The most common reconstruction approaches are filtered backprojection (FBP) and iter-ative reconstructions. 1.6.1 Filtered Backprojection Filtered backprojection is the most common method used for clinical reconstructions. The major advantage of this method is its speed; reconstruction of a 2D slice by FBP takes only a fraction of a second. The idea behind backprojection is very simple. The value of each bin in the projection is back-projected onto the image along a line perpendicular to the projection. Superposition of these lines from different projections give the location of the source. Mathematically in a continuous, space such a operation can be described by: where /(#, y) is the reconstructed source distribution, and p(t, 9) is the value of projection at angle 6 and position t = x cos 6 + y sin 9. Simple backprojection is not enough for useful reconstruction, because it blurs the activity across the image, producing star artifacts (Fig. 1.4). To avoid that, projection data p(t,9) have to be convoluted with some filter function iu(r) as in Eq. 1.2. Filtered (1.1) Chapter 1. Introduction 13 (a) (b) (c) Figure 1.4: Schematically described backprojection. (a) slice of the object and 3 sample projections, (b) simple backprojection results in star artifacts, (c) filtered backprojection removes most of the star artifacts around the object. backprojection is then the operation of convolution of the projections with a filter function W(T) and then backprojection according E q . 1.1. The procedure of filtered backprojection as described here is done in image space; it could also be done in frequency space. For further information see [3, 4, 5, 6]. Mathematically, filtered backprojection ( F B P ) can be described by: p(T,6)-w(t-T)dTd0 -oo (1.2) involution Chapter 1. Introduction 14 F B P is usually done with a 180° data set; a 360° set can be reduced to a 180° set by taking a mean from opposite projections. Iterative methods of reconstruction have become more and more popular because they produce images which are superior to those reconstructed with F B P [7]. Although they require much more computational time, the low cost of high power computers wi l l eventually make the iterative methods available not only for research but also for clinical use. Iterative reconstruction starts with some ini t ia l estimate of the source distribution. A set of projections P is calculated from this estimate by E q . 1.3. The vector A is the current estimate of the source. The size of this vector is equal to the number of pixels in the reconstruction space (e.g. for a 64x64 reconstruction, this size is 4096). The projection elements Pj and Pj are, respectively, the values measured during the scan and those calculated from the current estimate of source distribution. The number of projection elements is equal to the number of projections times number of bins in the projection. The matrix w is a system matrix. The element Wjti of w is the probability that a photon emitted in the i-ih pixel wi l l be detected in the j-th projection bin. In 2D reconstruction, the element of the system matrix is proportional to the common area of pixel i with a parallel strip coming from the projection bin j (Fig. 1.5). The number of elements in such a matrix is very high; for a reconstruction of 128x128 pixels, 64 projections, and 128 bins per projection, the number of elements is 134,217,728, but fortunately the system matrix w is sparse (most elements are equal to 0), so either w 1.6.2 Iterative Methods Chapter 1. Introduction 15 Reconstruction area divided into pixels Figure 1.5: Definition of system matrix w. Element u > , j is the area of overlap of pixel i with the parallel strip coming from bin j. can be compressed and stored, or non zero elements of w can be calculated (with some approximations) on the fly during reconstruction. The calculated projections P are compared with the measured projections P, and then the image estimate is updated according to some recursive formula. This process is iterated until it meets some preselected criterion of optimal solution. The criterion could for example be that of least square error [8, 9, 10], where the function of differences between measured and calculated projections is minimized. Another approach would be to treat Eq. 1.3 as a simple inverse problem of the set of linear equations. Many different algorithms can be used to solve this problem, for example singular value decomposition (SVD) [11, 12] and algebraic reconstruction technique (ART) [13, 14, 15] or its variations: multiplicative-ART, Dykstra-ART [14, 16]. Chapter 1. Introduction 16 The method best suited for the reconstruction of S P E C T images is maximum likeli-hood [17, 18] because of its statistical nature. This method models Poisson noise which is always present in the experimental projection data P . The drawback of maximum likelihood is that this method has very slow convergence when used with the expectation maximization (MLEM) algorithm. A recently modified form of M L E M , the maximum likelihood Ordered subset [19], has become very popular because its convergence is about an order of magnitude faster than M L E M . This algorithm, O S E M , is starting to be used in some commercially available SPECT systems. 1.7 Basic Physics of Photon Attenuation Radioactive isotopes used in nuclear imaging produce gamma ray photons. The photons emitted from inside the body of a patient travel throughout the tissue before they are detected in the camera, and on their paths to the detector they may interact with the tissue. There are 3 main types of photon interaction which are important at the energies used in nuclear medicine (from 60 keV to 500 keV). (a) P h o t o e l e c t r i c Ef fec t . The photoelectric effect takes place when a photon is com-pletely absorbed by the electron. It usually produces a vacancy in the atomic shell because the electron gains enough energy to be ejected from the atom. The pho-toelectric effect is accompanied by the emission of characteristic X-rays and Auger electrons by the atom. The cross section for photoelectric interaction is proportional to the 5th power for atomic number of the low atomic number materials. (b) R a y l e i g h S c a t t e r i n g . In Rayleigh (coherent) scattering, an incoming photon inter-acts with the whole atom. The change of its energy is negligible, only the direction Chapter 1. Introduction 17 is changed. (c) C o m p t o n S c a t t e r i n g . The Compton (incoherent) scattering is the most important interaction of photons in the human body. At an energy of 140 keV (corresponding to the decay of 9 9 m T c ) , 97% of interactions are due to Compton scatter interactions. In incoherent scattering, a photon interacts with an electron. During the interaction the photon changes its initial direction and part of its energy is transfered to the electron. 1.7.1 Linear Attenuation Coefficient The linear attenuation coefficient is defined as follows. When a monoenergetic narrow-beam (Fig. 1.6(a)) photon flux, $, passes through the attenuator, the fraction, d$, which is attenuated inside the a material of thickness dx is described by the formula: The proportionality factor \x is the linear attenuation coefficient. When photons pass through matter, any of the types of interactions described in Sec. 1.7, can take place independently, so the total ji is the sum of the linear coefficients from all three types of interaction. Integrating Eq. 1.4 for a non-homogeneous attenuator of thickness d, we obtain the photon flux $ which passes through the attenuator: d$ = -jj, • $ • dx. (1.4) where $o is the initial photon flux before attenuation (Fig. 1.6(b)). Chapter 1. Introduction 18 (a) dx d A f •mm i Source Figure 1.6: (a) Schematically described narrow beam geometry, (b) definition of the linear attenuation coefficient 1.8 Effect of Attenuation on Reconstructed SPECT Images Attenuation has a great impact on a reconstructed S P E C T image. Fil tered Backpro-jection which is a standard reconstruction technique used in clinical practice does not take into account the effect of attenuation of photons. Two kinds of attenuation artifacts can be found in reconstructed images. First is the apparent decrease of the reconstructed intensity in the areas lying deeper inside the body; the radiation coming from these areas is more attenuated than radiation coming from outer parts of the body. F i g . 1.7 shows that the average attenuation is much Chapter 1. Intra 19 60 50 ^ 40 o ffl 30 © CL 5 20 10 UJ 10 0 -15 -10 -5 0 5 10 15 Position on the diameter [cm] Figure 1.7: Fraction of radiation escaping from a water cylinder 30 cm in diameter con-taining a point source vs. position of the point source on the diameter. The fraction is an average calculated over 360° larger when a source of activity is positioned deeper inside the attenuating object than when it is closer to the surface of the object. For example, this difference causes the decrease of relative counts in the anterior wall of the heart for 2 0 1 TI perfusion studies [20]. It is apparent that such a decrease will vary depending on the anatomy of a patient (e.g. the existence of breasts in female patients and their absence in male patients) [21]. In Fig. 1.8 the effect of attenuation is shown in a reconstruction of a heart phantom. Although the heart and background were uniform in the phantom, they are not in the reconstructed image. Other artifacts caused by attenuation are the geometrical distortions of the shapes of reconstructed objects due to inconsistencies in the attenuated projections. In the absence of attenuating medium, a point source will be seen as having the same intensity in all the projections (Fig. 1.9(a)). FBP will reconstruct such an object without artifacts. When the point source is placed inside an attenuating medium, its intensity will be different in different projections due to attenuation (Fig. 1.9(b)). Reconstruction of such projections Chapter 1. Introduction 20 Figure 1.8: Reconstructed heart phantom with its profile. The bright area is the uniform heart. The reconstructed heart has darker and brighter areas due to attenuation, the part of the heart deeper in the phantom is darker. The reconstructed warm background also appears to be not uniform due to attenuation. wil l create shape distortions in the reconstructed image [22, 23]. F i g . 1.10 shows these artifacts in the F B P reconstruction of 1 cold and 2 hot spheres in a warm background. 1 . 9 C o r r e c t i o n f o r A t t e n u a t i o n As was demonstrated in Sec. 1.8, attenuation has a great impact on a reconstructed S P E C T image. There is a great amount of work being done to correct for attenuation errors and a large number of correction methods have been developed 1. The best results are obtained with methods where the reconstruction algorithms have information about density distribution (attenuation maps) and use it to correct for attenuation. F ig . 1.11 shows the improvement of the image when the attenuation correction is used; the cold spot is visible, the hot spheres are rounder, and the background is much more uniform. xFor further reference to these methods please see to chapter 7. Chapter 1. Introduction 21 Figure 1.9: The effect of attenuation on the number of counts in projections for (a) source without attenuation and (b) with attenuation. Figure 1.10: Attenuation artifacts in a reconstructed Jaszczak phantom. The Jaszczak phantom is a water cylinder, 22 cm in diameter filled with uniform background activity. There are 1 cold and 2 hot spheres inside the phantom. The cold spot is an air bubble; it is not seen at all due to severe attenuation artifacts. The shapes of the hot spheres are deformed, and the warm background of the phantom is not uniform due to streaks going from the hot spheres. 1 . 1 0 I m p o r t a n c e o f A t t e n u a t i o n C o r r e c t i o n As was shown in the example discussed in the last section, attenuation correction can greatly improve image quality by removing attenuation artifacts. Attenuation correction should be used also in clinical heart studies [24, 25, 26, 27]. It is hoped that quantitative S P E C T wi l l greatly improve diagnosis. A t the moment, S P E C T is only a qualitative technique. One of the main reasons for this is the lack of proper correction for attenuation. In order to do an accurate attenuation correction the density of the imaged object Chapter 1. Introduction 22 (a) (b) (c) Figure 1.11: Attenuation correction reconstruction of the Jaszczak phantom, (a) no correction (b) with correction (c) attenuation map used for the correction of this recon-struction. (attenuation maps) has to be known, and this work is about methods of acquiring these attenuation maps in a simple, inexpensive, and, at the same time, accurate way. Chapter 2 Different Methods of Measuring Attenuation Maps There are two approaches to the measurement of patient specific attenuation maps. The first is to use the information about tissue density determined by some technique different than S P E C T , like X-ray C T . The other is to use simultaneous or sequential acquisition of the transmission and the emission data with the S P E C T camera. The transmission scan is done using a source of radiation external to the scanned object which is placed on the opposite side of the patient from the detector. The detector measures the radiation transmitted through the patient. After the scan, attenuation maps can be reconstructed from the ratios of this transmitted radiation (transmission scan) to the radiation without the patient (blank scan). The following chapter contains a review of the methods which can be used to determine attenuation maps. Short descriptions of the geometry as well as the advantages and disadvantages are given for each of the different transmission source solutions. 23 Chapter 2. Different Methods of Measuring Attenuation Maps 24 2.1 Computed Tomography (CT) - X-rays C T employs X-rays for determining the density distribution maps of the human body [28, 29]. Advantages: • Computed Tomography provides very high resolution (about 1 mm) density distri-bution maps. Disadvantages: • C T density maps are in units of C T numbers; they have to be translated into attenuation coefficients. It is not an straightforward task because C T maps are obtained by polyenergetic X-rays and therefore measure average values for / i . As the X-rays pass through the patient low energies are attenuated more strongly (beam hardening); the energy spectrum of the beam changes with the depth which results in a different values for fi. • C T attenuation maps are acquired on different equipment than is used for SPECT, so they have to be coregistered with the SPECT image. Image coregistration is a complicated procedure and cannot always be done correctly (e.g. The position of the body of a patient can be different for the scans). A research system has been built for simultaneous C T and S P E C T [30]. Such a system, however, will not be widely available for some time [31]. • A C T scan represents a separate clinical test which is not always available. The following sections will discuss methods of determining attenuation maps with the use of transmission sources on nuclear medicine cameras. Chapter 2. Different Methods of Measuring Attenuation Maps 25 2.2 S h e e t S o u r c e This type of a transmission device is built with a uniform radioactive sheet source placed opposite to the detector equipped with a parallel hole collimator [32, 33, 34, 35, 36, 37, 38]. The sheet source can be used with and without a collimator in front of the source (Fig. 2.1). hole collimator Figure 2.1: Configuration of the sheet source (a) without collimator on the source and (b) with collimator. For the sheet source without collimator (a) a lot of scatter radiation can be detected changing the resulting attenuation coefficient (broad beam effect). This scatter is mostly eliminated by placing a collimator in front of the source (b) Advantages: • It is a simple solution; it does not require any special hardware. • The entire field of view ( F O V ) of the detector is used for detecting the transmission data. Even big objects can fit into F O V of this transmission source (no truncation Chapter 2. Different Methods of Measuring Attenuation Maps 26 artifacts). • The detector collecting the transmission data can also collect emission data in dif-ferent energy windows; it does not have to be dedicated to transmission data only. • The distance between the source and the object (patient) can be changed without any adjustment to the position of the transmission source. This is important for the simultaneous transmission/emission scans so that the detector can be placed as close to the patient as possible. D i s a d v a n t a g e s : • The transmission data measured by a sheet source without a collimator is contam-inated by a high level of scattered radiation (Fig. 2.1(a)). • The dose to the patient from a non-collimated source is much higher than from a collimated one. • The collimator added in front of the sheet source eliminates the previous disad-vantages and improves spatial resolution [39], but increases dramatically the size of the device and its activity. The dose to the technologist may be substantial, and handling of the source may be difficult. 2.3 Scanning Line Transmission Source The scanning line source [40] is currently manufactured by a few companies and it is considered to be the best commercially available transmission source geometry. As illustrated in Fig. 2.2 a single line source is moved along the detector in the direction Chapter 2. Different Methods of Measuring Attenuation Maps 27 moving collimated line source Figure 2.2: The scanning line source. During each S P E C T projection the whole area of the detector is scanned. A narrow area on the detector opposite the line source (shaded on the figure) is electronically windowed to detect only the transmission data; the rest of the detector acquires the emission data. perpendicular or parallel to the axis of rotation of the camera which has a moving elec-tronic window to detect only the transmission events. At each stop of the camera, the line source scans the whole area of the detector. Advantages: • The entire field of view (F0V) of the detector is used. • The detector collecting transmission data can also collect emission data in a different energy window. It does not have to be dedicated for transmission only. • Source collimation eliminates transmission scattered photons. • The windowing of the detector reduces contamination (cross-talk) of transmission Chapter 2. Different Methods of Measuring Attenuation Maps 28 data by emission radiation, and eliminates contamination of emission data by trans-mission radiation. Disadvantages: • The scanning line source requires additional electronic and mechanical hardware for moving the line and for windowing the detector. • Activity in the scanning line cannot be too high, otherwise there is a problem of dead time of the detector when there is no attenuating medium between the source and the detector. At the same time the source cannot be too weak because for big highly attenuating patients there would not be enough counts in transmission data to reconstruct good quality attenuation maps. It implies that the line source has to be replaced very often which increases the cost of the device. This problem can be diminished in some designs by using a graded filter system to reduce the flux when the source is hot and then to increase the flux (by reducing the filtration) as the source decays. 2 . 4 F a n a n d C o n e B e a m G e o m e t r y T r a n s m i s s i o n S o u r c e This kind of a transmission source can be implemented for any camera, but it has found the most success in triple head systems [41]. A line source is placed at the focal distance of the converging fan-beam detector collimator [26, 42, 43, 41, 44]. Also cone beam geometry can be used where a point source is positioned at the focal point of the converging cone beam collimator [45]. Advantages: • High efficiency of converging beam geometry due to low collimation of the source. Chapter 2. Different Methods of Measuring Attenuation Maps 29 Figure 2.3: The fan beam geometry of the transmission source. A single line source can be positioned in the centre (left), or to reduce truncation artifacts, off centre (right). The line source has to be located at the focal position of the converging collimator. • Attenuation maps obtained from a line source with a fan beam collimator have good spatial resolution. • A converging collimator reduces the scatter fraction in the transmission data [46]. D i s a d v a n t a g e s : • The truncation effect. As it can be seen in Fig. 2.3, the F O V of the convergent collimator is much smaller than the FOV of the detector with a parallel hole col-limator. It may result in severe truncation artifacts [47] in attenuation maps when large objects are being scanned. This effect is even more prevalent for the cone beam geometry [48]. • The line source has to be at a fixed distance from the detector at all times. 2 . 5 L i n e S o u r c e s p e r p e n d i c u l a r t o t h e a x i s o f r o t a t i o n A transmission source consisting of four line sources has been investigated in the Chapter 2. Different Methods of Measuring Attenuation Maps 30 past [49]. The sources were mounted perpendicular to the axis of rotation of the camera (Fig. 2.4). A x i s of Line sources U "*~ rotation Figure 2.4: The line sources perpendicular to the axis of rotation. A d v a n t a g e s a n d D i s a d v a n t a g e s are not fully known. This transmission source was never completely investigated, and the advantages and the disadvantages of that solution were not fully discussed. Because of the limited success of this design its further investigation, to the best of our knowledge, was dropped. 2 . 6 C o n c l u s i o n In this chapter the different types of geometry for a transmission source, together with their advantages and disadvantages were presented. The scanning line source is considered to be the most promising design, but it suffers from high cost and complicated design. In this dissertation a new transmission source which uses a series of collimated line sources will be discussed. Chapter 3 Multiple Line Transmission Source Geometry Our novel transmission source system [50, 51, 52, 53, 54, 55, 56, 57], utilizes a series of parallel collimated line sources. It is mounted opposite to a detector equipped with a standard parallel hole collimator and rotates together with it. The line sources are parallel to the axis of rotation of the camera. (Fig. 3.1) Advantages: • This design is simple and does not require any complicated hardware. • The F O V of the transmission source is the same as F O V for the emission scan. There are no truncation artifacts. • It allows the activity distribution in the system to be optimized to the density distribution in the human body; with the strongest sources placed at the centre 31 Chapter 3. Multiple Line Transmission Source Geometry 32 where there is the most attenuation and the weakest ones at the edges where there is the least. • Not all of the line sources have to be replaced as the activity of the source decays. If the ratio of the initial activity between the two neighbouring line sources is the same for all lines, then, when the sources decay by this ratio, restoring the initial activity requires only the two central strongest sources to to be replaced, and the other sources moved to the next, outside positions. • The distance between the source and the detector can be changed. • Source collimation minimizes scatter. Two types of multiple line sources were designed and tested. (a) (b) line sources axis of rotation • 0 0 - $ detector Figure 3.1: Geometry of multiple line sources. Collimated multiple line sources are placed parallel to the axis of rotation of a SPECT camera, (a) view perpendicular to the axis of rotation (b) view along the axis of rotation. Chapter 3. Multiple Line Transmission Source Geometry 33 3 . 1 C o l l i m a t e d L i n e S o u r c e s ( C L S ) The proposed CLS transmission source system consists of a series of highly collimated capillary tubes filled with activity. The whole system is mounted opposite to a detector equipped with a standard parallel hole collimator. The line sources are parallel to the axis of rotation of the camera. In order to minimize the number of line sources and optimize resolution, we use a unique geometry where the sources are highly collimated in the transaxial plane. Due to this, only half of the total area of the investigated object is illuminated at any time. Measurement of the complementary half is performed when the whole system is rotated by 180°. This geometry requires that the line source system be located so that the axis of rotation of the camera is on the edge of one of the central fan beams and that a full 360° is scanned (see Fig. 3.2). The activity in the central lines is the highest, and the ratio of activity between neighbouring lines is constant (Fig. 3.2). A special collimator on each line source ensures that the radiation forms a fan beam which falls well inside the acceptance angle of the detector collimator. The spacing between adjacent sources, the distance between the camera and the source and the source collimation are all closely related and are selected so that the fan beams illuminate the surface of the detector without overlapping. Since it is the source collimation which limits the width of the fan beams, any collimator with an acceptance angle larger than the source collimation angle may be used on the camera. The collimation angle of the line source is defined as the maximum angle between the direction of photons in the fan beam and the normal to the detector surface. The geometry used in our studies permits the use of all parallel hole collimators which are currently available in our nuclear medicine department. Figure 3.2: (a) The geometry of the collimated line sources (CLS) transmission source. The simulated (b) and experimental (c) profiles of the blank scan photon beam. Chapter 3. Multiple Line Transmission Source Geometry 35 3 . 2 M u l t i p l e L i n e A r r a y The multiple line array is a modified version of the CLS system. The principle of the MLA is similar to that of the CLS. The source is also built with a series of collimated line sources, but the distances between the lines vary, and the axis of rotation of the camera is not positioned half way between the source and the detector. In this design the collimator on the line sources has an collimation angle of less than 3° since then the photon beam is mostly perpendicular to the detector surface which allows for substantial elimination of scatter. The activity distribution in the M L A system follows the same principle as for CLS with the strongest sources in the centre and weaker ones in the outside positions. The activity in the lines increases exponentially going towards to the centre of the array, and the ratio of activity between the neighbouring lines is constant. The distances between the lines are chosen to give a smooth distribution of the beam density on the detector (Fig. 3.3). The scheme of the geometry of the SPECT camera with the M L A transmission source and the profile of the resulting photon beam are presented in Fig. 3.3. Unlike CLS where a full 360° scan was necessary, M L A requires only a 180° rotation. Figure 3.3: (a) The geometry of the multiple line array ( M L A ) transmission source (only 8 line sources are shown). The simulated (b) and experimental (c) profiles of the blank scan photon beam. Chapter 3. Multiple Line Transmission Source Geometry 37 3.3 S u m m a r y A novel geometry of the transmission source for S P E C T has been proposed. The collimated line source is built from a series of collimated line sources parallel to the axis of rotation. The sources are equally spaced, and the radiation coming from neighbouring lines do not overlap on the detector. This transmission source requires 360° rotation of the camera to acquire a full set of data needed for the reconstruction of the attenuation maps. The second design, the multiple line array also uses a series of line sources placed parallel to the axis of rotation of the detector. Here however sources are not equally spaced, and their fan beams do overlap on the detector producing a relatively smooth distribution on the detector. The multiple line source systems have many advantages, such as simplicity, low main-tenance, and the possibility of the activity in the source being tailored to the attenuation of the human body. Chapter 4 Methods of Reconstruction of the Attenuation Maps In this chapter, are described the approaches to the reconstruction of attenuation maps from data acquired by multiple line sources for both CLS and M L A source configurations. Since the CLS and M L A systems have slightly different geometries, they require different pre-reconstruction and reconstruction methods. 4 . 1 C o l l i m a t e d L i n e S o u r c e s - A t t e n u a t i o n M a p R e -c o n s t r u c t i o n The photon flux generated by CLS corresponds to the geometry of the multiple fan beams. Because the fan beams do not overlap on the detector, they produce a highly non-uniform distribution of counts on the detector (Fig. 3.2). There is a very low number of detected counts in the areas of the detector where the fan beams "touch" each other. Reconstruction of attenuation maps directly from the data with this low number of counts 38 Chapter 4. Methods of Reconstruction of the Attenuation Maps 39 would result in large errors due to low statistics. It may also happen that no counts are detected. This lack of counts creates a missing count problem. In order to avoid the problem of the missing counts for CLS, the following special rebinning procedure can be used. Although the goal of performing the rebinning procedure is to remove from the transmission data the bins with low count and bins where counts can be detected from unknown line sources (areas where the beams "touch"), the rebinning procedure controls also the resolution of attenuation maps. 4.1.1 Data Preprocessing - Rebinning The rebinning is performed only in the direction perpendicular to the line sources. Its principle is illustrated in Fig. 4.1. This special procedure is only applicable to the CLS system. During the scan, counts are collected into bins; the bin size is determined both by the matrix and the detector size used in acquisition. We rebin these initial data acquired in a standard matrix into a new set of bins. The position and size of these new bins are determined by the geometry of the transmission system. The size of a new bin corresponds to the detector area illuminated by one line source, or a fraction of this area equal to the size of this area divided by an integer. For example (Fig. 4.1) after collecting the data from 10 line sources into 128 bins, we can rebin the 128 bins into 10 new bins since there are 10 line sources, or we can rebin them into 20 new bins if we divide each area illuminated by one line source into 2 divisions, 30 if into 3 (as it is shown in Fig. 4.1) and so on. As a result of this transformation we have a new set of data with a running parameter {j} where j = {1,..., Ipr • Iunes • Idivision}, where Ipr is the number of projections, Iunes is the number of line sources and Idivision is the number of divisions. In our method, we rebin the counts recorded in detector bins for the blank scan B Chapter 4. Methods of Reconstruction of the Attenuation Maps 40 new bins Figure 4.1: Pre-reconstruction rebinning. In this example each fan beam was divided into 3 smaller fan beams (divisions). (scan taken in the absence of attenuator) and for the transmission scan I . Mathematically, the rebinning is described by the following: BTW = £ *kj • Bk 1J™> = £ Skj • h (4.1) k k The value of Skj is the fraction of the bin k, contained inside the new bin j . 4.1.2 Strip Geometry - Filtered Backprojection Filtered Backprojection is the easiest and fastest method which can be used for the reconstruction of attenuation maps obtained from the CLS system. In order to use FBP we have to make a few assumptions and further transform our data. We have to assume that the transmission data at the detector comes from a parallel strip. For each new bin, taking values of B^ew and Tjew we can calculate the log transform Chapter 4. Methods of Reconstruction of the Attenuation Maps 41 value Aj given by: Aj = log(Bfwmew) (4.2) The value of Aj corresponds to a line integral of the attenuation taken along the line perpendicular to the detector and originating at the bin j. In order to reduce the error which we are making doing this approximation, we in-crease the density of new bins by combining each projection with its opposite projection (see Fig. 4.2). The size of the new bin is "squeezed" to half of its size, and grouped at the centre of its fan beam. Gaps created are filled with values of A from the "squeezed" new bins from the opposite projection. Because of this transformation, a 360° initial data set becomes a set of 180° data, and the number of projections is halved. Next, the attenuation maps are reconstructed by FBP from the transformed A data. In Fig. 4.2 the real shape of the radiation coming from each line source (fan beams) is compared with our approximation (parallel strips). We can see that the parallel geometry approximation is worst near the source and detector and best at the centre, where scanned object are actually positioned. This approximation is valid for any number of divisions. Chapter 4. Methods of Reconstruction of the Attenuation Maps 42 (a) (b) Figure 4.2: Strip geometry in F B P . (a) The first two images show the two opposite projections which wi l l be combined into one data set. The shaded strips represent our approximation of the fan beams (solid lines). The third image is the combination of these two opposite projections. The dot in the middle represent the axis of rotation. This procedure can be performed on data where the ini t ial beam is divided into any number of subdivisions, (b) Rebinning procedure where the beam is divided into 2 divisions. Chapter 4. Methods of Reconstruction of the Attenuation Maps 43 4.1.3 Exact Geometry - Iterative Approach The attenuation maps from the CLS system may also be reconstructed using iterative methods which can take into account the exact geometry of the fan beams. In all iterative methods we need to define the projection operation. This operation calculates from the current estimate of the reconstructed image the resulting projection elements which can then be used in iterative reconstruction. In other words, the projection operation creates the relation between the reconstructed image and the measured projection elements. We will first define the projection operation for emission S P E C T , before discussing transmission imaging, since the emission case is simpler and more intuitive. I t e ra t ive R e c o n s t r u c t i o n o f E m i s s i o n S P E C T Images The goal of the reconstruction in emission S P E C T is to find the density of the radio-active material inside the patient, emission density, from acquired projections. Let A be a digitized estimate of the true emission density A. If Wij is the probability that a photon emitted from pixel i is detected in bin j , the resulting value of projection element pj can be expressed by: Pj = ^2wi,i- ^i- (4-3) i The matrix of these probabilities, w, is the system matrix which describes the conditions of the data acquisition (system geometry, detector collimator, Compton scatter, attenu-ation). In a simplistic model of 2D SPECT data acquisition (without taking into account collimator blurring, scatter and attenuation) for a parallel geometry, w can be calculated as a value proportional to the area of overlap of the pixel i and a strip corresponding to the width of bin j (see Fig. 1.5 on page 15). Chapter 4. Methods of Reconstruction of the Attenuation Maps 44 Figure 4.3: Method of calculating system matrix w for multiple fan beams. I t e ra t ive R e c o n s t r u c t i o n o f T r a n s m i s s i o n Images for C L S S y s t e m Let us derive the relation between the counts I measured in the transmission scan, and the attenuation map, / i . In order to find this relation an additional scan has to be performed, a blank scan B . Let $(x) be the photon density per unit length of the detector at position x (Fig. 4.3) with attenuation, and $Q(X) be the photon density with no attenuation. Using the basic Eq. 1.6 on page 17 for photon attenuation and digitizing this formula we have $(x) = $o(x) • exp(-/,-(x) • m) = $0(x) • exp(^ -U{x) • m) (4.4) i i where U(x) is the length of the intersection of the line between point x on the detector and the line source, S, with the pixel i. If we define A(x) = \og($>0(x)/<&(x)), Eq. (4.4) becomes: Chapter 4. Methods of Reconstruction of the Attenuation Maps 45 A(x) = Ylli(x)-(ii. (4.5) i Calculating an average value of A(x) inside bin j where bin j extends from x\ to x2, we obtain: 1 f%2 Aj = / 22 li(x) ' • dx. (4.6) X2 — Xi Jxi • By exchanging the order of integration and summation, and replacing x with x' (simple scaling by ^ (see Fig. 4.3)) we obtain: D p ' a Ai = £ / * t -j- I -U{x') • dx' j X2 — Xi dij Jx'i (4.7) where D is the distance from the line source to the detector, and dij is the distance from the middle of pixel i to the line source S. Because of the small angle of the fan beam we will use the approximation X 2 ~ a. The integral f £ i * -li(x') • dx' is overlap a,-fj of the pixel i with the fan beam j (fan beam j is presented in Fig. 4.3 as a triangle (xi, x2, S)). Taking this into account we have a final expression for Aj, where: Aj can be approximated by: AJ = £ wi,j • Vi (4.8) B • A , - ~ l o g - ^ (4.10) where Ij and Bj are the number of counts in bin j for scan with and without attenuation respectively. Verification of 4.10 is included in Appendix D. Using Eq. 4.8 and Eq. 4.10 we now have Chapter 4. Methods of Reconstruction of the Attenuation Maps 46 Ij = Bj • exp(- (4.11) i where the system matrix w is as defined in Eq. 4.9. Eq. 4.11 gives the relation between our measurement I and the value of u. which is being reconstructed. Eq. 4.11 also gives the equation for calculation of projection elements I from the current estimate of the attenuation map, /}, which will be used later in this chapter. Calculation of the matrix w is a challenging task because it is difficult to find a general expression for the overlap between a square pixel and a triangular fan beam. The size of w is large (50Mb) but it can be compressed because of its sparseness. The method and the code for calculation and compression of w are included in Appendix A. R e c o n s t r u c t i o n M e t h o d s - M a x i m u m L i k e l i h o o d The maximum likelihood [18, 58, 59] reconstruction method uses Poisson statistics to model the transmission data. Let Ij be the number of counts detected in the jth bin of the transmission scan, and Ij be the mean of Ij. Since Ij is a Poisson variable and Ij is unknown, the value of Ij has to be estimated in order to write down the likelihood function for the transmission measurements. Using Eq. 4.11, this estimate can be written; Ij = Bj • exp(- wi,j • Vi) (4-12) i where Bj is the mean of the counts acquired during the blank scan. Bj is also unknown, but it can be well approximated by the number of counts Bj measured during the blank scan. The approximation is good because Bj is usually a large number, so the Poisson fractional deviation of Bj from the mean is small. Assuming that: Bj = Bj (4.13) Chapter 4. Methods of Reconstruction of the Attenuation Maps 47 we get Ij = Bj • exp(- Y • hi)- (4.14) i Eq. 4.14 gives an equation for the calculation of the mean of the Poisson variable Ij. It is known [60] that if a given process can be described by a Poisson probability distribution, the probability that Ij events will take place with a mean Ij is equal to: I-Ij Pj(Ij I h) = exp(-/,-) • -jj. (4.15) Because Ij is dependent on Bj, a constant, and on \x (which is to be found) we can express Eq. 4.15 as: Wi I ft = Wi I A*)- (4-16) The likelihood function for the system will be: I(I|/0 = nW;l^)- ( 4- 1 7) 3 This likelihood function expresses the probability of recording I counts in the detector, where I is a vector of Poisson variables. This probability is dependent on p. The maximum value of likelihood corresponds to a configuration of \i which gives the values of I with the greatest probability. Eq. 4.17 is complicated; its full form based on Eqs. 4.14, 4.15, and 4.17 is: 1(1 | H) = IW-S; • e x P ( - X>* • A*.")) • ( ^ " e X P ( ~ ^ ^ " ^ ) ) J j - (4-18) Eq. 4.18 can be simplified by taking its logarithm which leads to the log-likelihood function 1(1 | H) = In L(N | fi) = }Z(-B3 • e x P ( - wh3 • ^) ~ h Yl wi,j • ^ + h MBi) ~ l n ( 7 r ))• 3 i i (4-19) Chapter 4. Methods of Reconstruction of the Attenuation Maps 48 The last two terms in Eq. 4.19 are constant. Finally, neglecting the two constant terms, the maximum likelihood objective function for a transmission scan is: FMM = 52(-Bj • exp(- Wij • Pi) Wi,j • Pi). (4.20) 3 i i The maximum of the log-likelihood function can be estimated giving a maximum likelihood solution of p. R e c o n s t r u c t i o n M e t h o d s - L e a s t S q u a r e s Another method used for the reconstruction of attenuation maps is least squares (LS)[8, 61]. In this method the objective function, the function which is being minim-ized, is the sum of the squares of differences between the measured values and the values projected (through the projection operation) from the current estimate of the image. The objective function is thus: ^ ) = E(1°S(^)-E^-^) 2- (4-21) i *i i Eq. 4.21 does not take into account the statistical character of the measured data. The value Ij is measured with Poisson error. In order to take this error into account, the contribution to the objective function from a given detector bin can be weighted with the reciprocal of variance of log(y1). Assuming that Ij is a Poisson variable and neglecting error of Bj we have: Var(\og(B3/I3)) = 1 (4.22) The weighted least squares objective function (WLS) becomes: Chapter 4. Methods of Reconstruction of the Attenuation Maps 49 FWLS{H) = T, h -(M^-E^-^)2- (4-23) Weight P e n a l t y T e r m Due to the presence of noise in the measured data, as well as some geometrical incon-sistencies (for example 2D reconstruction from 3D data acquisition), the inverse problem of Eq. 4.11 is ill conditioned. The unregularized methods of reconstruction produce noise levels which increase with increasing numbers of iterations [62, 63]. To remedy this prob-lem, several regularization methods can be used to impose smoothness constraints on the image estimate. One such method is to incorporate a smoothness penalty or a prior [64, 65, 66]. This approach is straightforward with M L and WLS methods. The penalized objective function has a form: FP(p) = F(p) + B • R(p) (4.24) where F(p) is non-penalized objective function, R(p) is a penalty term, and 3 is a penalty parameter. Generally the effect of the penalty term R(u) is to discourage big differences between neighbouring pixels while F(p) tends to make an agreement with measured data. Usually these are conflicting goals. The penalty parameter 3 controls the trade off between these two tendencies. Many different priors have been investigated [67, 68, 69, 70]. In this dissertation a quadratic penalty was used: RM = \]£E*,k\(l«-l*k)2. (4.25) The value of the matrix a is described in Fig. 4.4. Chapter 4. Methods of Reconstruction of the Attenuation Maps 50 a u = a i , 3 = a i > 6 = a i r 0 J l ai,2=aiA=airauTl'° Figure 4.4: A schematic presentation of the matrix a^k- Horizontal and vertical neigh-bours contribute to the penalty with weight 1, diagonal with 0.71, and 0 otherwise. 4 . 2 M u l t i p l e L i n e A r r a y - M a p R e c o n s t r u c t i o n The reconstruction of the attenuation maps from the multiple line array is simpler than for CLS because of the simpler geometry. The reconstruction is based on parallel geometry which means that all lines of response, for each bin, are assumed to be parallel to each other and perpendicular to the surface of the detector. We cannot use the exact approach as for the CLS case because of the beams' overlap on the detector. For M L A , we have to approximate the geometry of multiple overlapping fan beams to the parallel geometry. 4.2.1 Data Preprocessing - Missing Counts Correction Algorithm The transmission projections for M L A may suffer from missing counts. The problem may be encountered when large patients, who are highly attenuating, are being scanned and no counts are recorded for some bins. The problem of missing counts causes serious artifacts in the attenuation maps. Our method of correction for missing counts uses an interpolation of the data in cubic volumes of 27 pixels. These 3D calculations use the projection data and are performed using X, Y and 0 coordinates. The algorithm "walks" through the data, starting from the 1 mm 3 a i .5 6 H 8 Chapter 4. Methods of Reconstruction of the Attenuation Maps 51 boundaries of the object where there is little attenuation and many counts are recorded, and moving towards the centre, where the attenuation is greater and the data are more scarce. It replaces the content of negative1 and/or low counts pixels with the the average of the 27 pixel cube. Although described procedure does not preserve the number of total counts, it was found that this method performed better than other standard methods [56] that can be used to remove bins with zero/negative counts. Counts are not preserved because for some experimental studies with large attenuators when the scatter correction was applied there were some large areas in the projections with average count less then zero, and preservation of counts in these areas had to broken in order to calculate the log-transform. 4.2.2 Filtered Backprojection As for the CLS case, filtered backprojection can be used for reconstruction of the attenuation maps. Before performing filtered backprojection, Ij and Bj counts acquired with and without attenuation are transformed by the log transform operation. The value of Aj = log Bj/Ij is the value of the line integral of the attenuation coefficient calculated along a line perpendicular to the detector and through the projection element j. The maps of the attenuation coefficients can be reconstructed from these line integral values Aj by a standard FBP. 4.2.3 Iterative Methods The same iterative methods were investigated for the M L A system as for the CLS system, namely the maximum likelihood and least squares methods. The only difference was that the system matrix w, which describes the geometry of the transmission source, Negative counts are possible when scatter correction is applied. See Section 7.1.2. Chapter 4. Methods of Reconstruction of the Attenuation Maps 52 was different for the M L A and CLS systems. In this section only the calculation of the system matrix for M L A will be described, for a description of iterative methods please refer to Sec. 4.1.3. System matrix The derivation of the system matrix w for M L A is very similar to the derivation of w for CLS. Let $(x) again be the photon density with attenuation per unit length of the detector at position x (Fig. 4.5), and $o(x) be the photon density detected with no attenuator. Using the digitized basic Eq. 1.6 for the photon attenuation we have: $(x) = $„(* ) II e x p H i ( x ) • (4.26) Since we have assumed parallel geometry of the system, the lines of response are per-pendicular to the detector surface. As for CLS, we define A(x) = l o g ( $ 0 ( x ) a n d Eq. 4.26 becomes: A{x) = Y^li(x)-fXi. (4.27) i Calculating an average value of A(x) inside a bin j , Aj (bin j extends from x\ to X 2 ) we arrive at: Aj = / Y li(x) • Vi • d x - (4-28) X2 — X\ Jxi • Interchanging the order of summation and integration and using the fact that -/,-(x) • dx is the area a;j of the pixel i inside the parallel strip coming from the bin j we have the final expression: = Y wiJ • V* (4-29) i where ai,j BinSize (4.30) Chapter 4. Methods of Reconstruction of the Attenuation Maps 53 Figure 4.5: The geometry of the M L A source. Schematically shown are one pixel of the reconstruction area, i, and one projection bin, j. and where BinSize is the length of the projection bin (x2 — xi). Now using again Eq. 4.10, we have the projection operation defined for M L A as: Ij =. Bj • exp(- ]T witj • m) (4.31) i with w defined by Eq. 4.30. Eq. 4.10 for the case of M L A is a much better approximation than it was for CLS because the bin size is much smaller for the M L A . There is one further point worth mentioning. The system matrix w is calculated exactly. The strip path integrals are used with no approximations being made during its calculation. Some simplified algorithms, such as the line path integral where W{j is approximated by the length of the intersection of the strip's centre line with the pixel, introduce errors during the projection operation [71]. Chapter 4. Methods of Reconstruction of the Attenuation Maps 54 4.3 S u m m a r y a n d C o n c l u s i o n s In this chapter the methods used for the reconstruction of the attenuation maps for both CLS and M L A were presented; filtered backprojection and iterative. These methods can be divided into two categories. Filtered backprojection is simple and very fast. Its drawback is that it does not account for the statistical character of transmission data. Iterative methods are more sophisticated than FBP. Effects like collimator blurring, attenuation, scatter, and even septal penetration can be incorporated into them. The effect of Poisson noise can also be taken into account during reconstruction, either in the log-likelihood objective function in the maximum likelihood method or by using weights in the WLS algorithm. Generally the image quality reconstructed by an iterative method should be better than that from FBP reconstruction providing enough iterations are performed [61, 7, 72]. The drawbacks of iterative methods are that they carry a much larger computational burden than FBP, and during the reconstruction, noise can blow up due to the ill conditioned nature of the tomographic reconstruction. The second problem can be diminished by using priors, but this adds even more computation time to the reconstruction. All of the iterative attenuation map reconstructions presented in this work use an accurate system matrix which can be precalculated for any given system and stored in a compressed form on a disk. It was shown [73] that mismatch in the resolution of attenuation maps and emission image can cause artifacts when doing attenuation correction especially in the regions with rapidly changing attenuation coefficients (e.g. myocardium and lung). Because of that reason, the resolution of attenuation maps reconstructed by the methods described in this chapter should be equal to the resolution of emission images. As it was mentioned before the rebinning procedure for the CLS system as well as choice of the reconstruction filter for Chapter 4. Methods of Reconstruction of the Attenuation Maps 55 FBP reconstruction or normalization factor for iterative methods trade off the resolution for variance reduction. Thus, the resolution of the attenuation maps can be adjusted by usage of the appropriate reconstruction method to resolution required to match this of the emission images. Chapter 5 Simulation Experiments Simulation experiments were performed for both source configurations. The goal of the tests was to check the feasibility of the proposed attenuation map generation techniques using both transmission sources, as well as to check the performance of the reconstruction methods. In simulations we have control of all simulated parameters of the system. We also know the exact geometry of the modelled object, which allows us to determine the accuracy of the reconstruction. 5.1 G o a l s o f t h e S i m u l a t i o n s We focused our simulations on the investigation of the following parameters of the transmission systems: • geometry, • reconstruction methods, 56 Chapter 5. Simulation Experiments 57 • total activity in the source, • number of divisions in the preprocessing (CLS only), • thickness of the line sources (CLS only), • distribution of activity in the source, and • spatial resolution. 5 . 2 M e t h o d s 5.2.1 Analytic Phantom The object used in the analytic simulations was a cylindrical torso phantom with an elliptical cross section of 30x22 cm. It contained an ellipsoidal lung insert, 5x4.5x6 cm, a cylindrical spinal cord with diameter 2.5 cm, and three small spherical objects in the centre with diameters of 1 cm, 1.5 cm, and 2 cm. Simulated attenuation coefficients were 0.15, 0.02, 0.30, 0.00 c m - 1 for the body, lungs, spinal cord and spherical inserts, respectively. The central slice of the phantom is presented in Fig. 5.1. Figure 5.1: Central slice of the analytic phantom. The two square areas are the regions of interest (ROIs) used later in the analysis. Chapter 5. Simulation Experiments 58 5.2.2 Data Generation All of the simulation experiments were done using the analytic transmission projection generator (ATPG) (Appendix C), a fast computer program which generates transmission projections from the transmission sources. It calculates the number of counts expected in a given detector bin and then adds to it Poisson noise of the desired level. A T P G simulates the detector collimator blurring, detector intrinsic resolution, physical thickness of the line sources (CLS only), and the transmission source collimation. The Compton scatter, septal penetration, and source collimator penetration were not simulated. For more details about A T P G please refer to Appendix C. In all the simulations the Compton scatter effect was not included. In experiment, some of the transmission photons can be scattered in the attenuator and still be detected due to poor energy resolution of Anger camera detector. The investigation of this effect was done in our group, and it showed that this effect was relatively small (about 3-5%). Because of this small amount, the scatter effects were not included in the simulations. The results of that scatter fraction estimation will be included in separate paper. S P E C T D a t a G e n e r a t i o n for C L S A system of 10 line sources placed 65 cm from the 40x30 cm detector was modelled for a Sophy DST camera. The spacing between the sources was 4 cm. The detector was equipped with a low energy high resolution (LEHR) collimator with a 3° acceptance angle. A total of 64 projections was simulated with a single head rotation over 360°. Data was acquired into a 128x128 matrix with a 0.34 cm bin size. Three different levels of total activity in the transmission source were simulated producing different count statistics in the projections. The Low statistics image corresponded to 1.4 GBq (37 mCi), the Medium to 4.8 GBq (130 mCi), and the High to 11.8 GBq (320 mCi). The simulated count levels Chapter 5. Simulation Experiments 59 in the projections corresponded to a 20 s/projection acquisition protocol. S P E C T D a t a G e n e r a t i o n for M L A Twenty line sources were placed 100 cm from the detector. Three different levels of activity were simulated corresponding to a total activities in the source of 1.4 GBq (37 mCi), 2.8 GBq (75 mCi), and 9.3 GBq (250 mCi). A set of 64 projections was taken over 180°. Each projection was acquired for 20 seconds. Counts were collected into a 128x128 matrix with bin size 0.48 cm, the bin size for a Diacam Siemens Camera. The collimator used on the detector corresponded to an L E H R collimator with a 3° acceptance angle. 5.2.3 Attenuation Map Reconstruction The data from the simulations was preprocessed as described in Sec. 4.1.1 for CLS and in Sec. 4.2.1 for M L A . The filtered backprojection (FBP), maximum likelihood (ML,PML), and least squares (LS, WLS, PWLS) (Chapter 4) methods were tested for attenuation map reconstruction. 5.2.4 Error Analysis The accuracy and precision of the reconstructions were measured using two regions-of-interest (ROI) (6x6 pixels) placed in the central slice of the reconstructed object (Fig. 5.1). ROI-1 was placed in the heart area, and ROI-2 was placed near the edge of the phantom to monitor changes in the resolution of the reconstructed image. The local accuracy of the reconstruction was measured using the following formula: Et=i \vi - hi\ Ef = 1 hi (5.1) Chapter 5. Simulation Experiments 60 where //; and (ii are the values of attenuation coefficient assumed in the phantom and obtained from reconstruction, respectively, and N is the number of pixels inside the ROI. A value of a equal to 0 corresponds to the reconstructed ROI values being exactly equal to the phantom values [74]. The precision was measured by the % root mean square (% rms) uncertainty using the equation: u = 100 • ^ (5.2) where p, is the average value over the ROI calculated from: 1 N iV i=l and v is the variance over the ROI: The global accuracy of the reconstruction A, calculated for the entire slice, was meas-ured by the Cartesian distance between the phantom and the reconstructed image, cal-culated using the formula: A = J T , (^-Pi)2 (5.5) V «',Mi>0 5.3 Results - Simulations of C L S 5.3.1 Geometry of the Simulated CLS One of the transmission source design goals was to have as few line sources as possible while maintaining a high density of counts on the detector and illuminating the entire area of the detector. In order to find an optimal number of sources, transmission images without an atten-uator were simulated. The distance between the detector and the source was chosen to be Chapter 5. Simulation Experiments 61 Position on the detector [cm] Figure 5.2: Simulated profiles of a one line source, (a) no collimation on the source (b) low collimation (c) high collimation 65 cm since this distance was later used in experiments with a Sophy DST camera. The profiles of the transmission fan beams (Fig.5.2) for a non- colli mated and collimated line source as measured by the detector equipped with a L E H R parallel hole collimator were simulated by A T P G . Fig. 5.2(a) shows that one non-collimated line source can illuminate about 6 cm of the detector surface, but that the number of counts in areas further than 2 cm from the centre of the beam is less than 25% of the peak. In the final design the distances between the line sources were chosen to be about 4 cm since larger distances would result in large areas with low counts, leading to problems with map reconstruction. The 4 cm distance between the line sources means that 10 equally spaced line sources are needed to illuminate a 40 cm detector. The source collimation has to be designed to create a series of non-overlapping fan beams. Such collimation would create the profile of the fan beam presented in Fig. 5.2(b). We can see (Fig. 5.2(b)) that due to the finite thickness of the line source, a penumbra is created on the detector (value of counts is not 0 at the position 3 cm and 7 cm). This penumbra will always be present in the experimental data, not only because of the definite thickness of the line sources, but also due to septal penetration. The effect might be decreased by more restrictive source collimation, as Chapter 5. Simulation Experiments 62 shown in Fig. 5.2(c), while still keeping 4 cm between the lines, but this would decrease the number of detected counts, and create some areas with no counts at all. Therefore, in our simulation, the source collimation was kept at 2° (Fig. 5.2(b)) producing a 0.5 cm overlap of the fan-beams. 5.3.2 Reconstruction Methods for CLS For each method the precision and, the local and global accuracies of the reconstructed attenuation maps were measured. Three divisions (see. 4.1.1) were used for the rebinning of the data; later (Sec. 5.3.3) it will be shown that 3 is the optimal number of divisions. F i l t e r e d B a c k p r o j e c t i o n The FBP reconstruction method with a Butterworth filter (order 5) [75, 76] and a cutoff frequency varying from 0.4 to 1.0 times the Nyquist frequency was used. Results of the reconstructions are shown in Fig. 5.3. As expected, increasing the activity in the source improved the accuracy and the precision of the reconstructed image. When a lower cutoff frequency was used, the image was smoother. There were opposing tendencies in the accuracy of the two ROIs. With a smaller cutoff frequency, the accuracy and precision of ROT1 increased because the image was smoother. On the other hand, a smoother image decreased the accuracy of ROT2 which was placed near the edge of the phantom. Some ringing artifacts due to overlap of the penumbras can be seen on the reconstruc-tions with cutoffs 0.85 and 1.00. These artifacts were masked by using stronger smoothing filter (Fig. 5.4). The presence of these artifacts was also seen in our earlier, less advanced, 2D simulations [52] which are not included in this thesis. Chapter 5. Simulation Experiments 63 l-J , , , , , , , , , , , | i • 1 ' 1 1——| 1 | 1 1 r \ 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.4 0.5 0.6 0.7 0.8 0.9 1.0 C u t o f f f r e q u e n c y C u t o f f f r e q u e n c y Figure 5.3: Accuracy and precision of the ROIs plotted vs. cutoff frequency. The top row corresponds to low activity in the CLS source, the middle, and the bottom rows to middle and high activities respectively. Chapter 5. Simulation Experiments 64 Figure 5.4: Central slices of the attenuation maps from the medium CLS source activity reconstructed using FBP with the Butterworth filter (order 5), with the cutoff frequencies of 0.45, 0.65, 0.85, and 1.00 respectively for images (a), (b), (c) and (d). Chapter 5. Simulation Experiments 65 ROI 1 0.0,-| . , . , . , . , . , . 1 4-| , , , , , , , , , , . 1 0 5 10 15 20 25 30 0 5 10 IS 20 25 30 ITERATION ITERATION Figure 5.5: The local accuracy and % uncertainty, plotted vs. number of iterations for the maximum likelihood algorithm using medium activity in the source. The 4 lines correspond to M L and P M L with different penalty parameter reconstructions. The value of penalty parameter is given (in the brackets) in the legend. M a x i m u m L i k e l i h o o d The transmission maximum log-likelihood (ML) objective function and penalized max-imum log-likelihood (PML) (Sec. 4.1.3) were maximized by the conjugate gradient method [77]. The changes of the local accuracy and precision for both ROIs vs. number of itera-tions are shown in Fig. 5.5. Resulting reconstructed images (20th iteration) are presented in Fig. 5.7. Fig. 5.5 shows that after about 20 iterations of the conjugate gradient method the preci-sion and accuracy changes for the two ROIs were negligible. Also the changes of likelihood were negligible with the number of iterations higher than 25 (Fig. 5.6). It strongly suggests that 20-25 iterations should be used for the reconstruction. Fig. 5.6 shows that addition Chapter 5. Simulation Experiments 66 -2835400 - , Figure 5.6: The value of log-likelihood plotted vs. number of iteration for maximum like-lihood reconstruction. The figure corresponds to a total activity in the CLS source equal to 4.8 GBq. The 4 lines correspond to M L and P M L with different penalty parameter reconstructions. NB: The plotted function is not an M L objective function but only the value of the log-likelihood. The value of penalty parameter is given (in the brackets) in the legend. Chapter 5. Simulation Experiments 67 (a) (b) ^ ^ ^ ^ ^ -60.72 • -51.94 • Figure 5.7: Central slices of the attenuation maps reconstructed by the maximum likeli-hood method with (a) 20 iterations of non-penalized maximum likelihood, (b) 20 iterations of PML with penalty parameter a = 0.4, (c) 20 iterations of P M L with a = 1, and (d) 20 iterations of P M L with a = 2. The total activity in the CLS source was 4.8 GBq. of the penalty term to the objective function decreased the reconstructed value of log-likelihood, and that, for a higher penalty parameter, conjugate gradient converged after only 40 iterations, unlike the case for the unconstrained (no penalty) algorithm. Some ring artifacts due to penumbra overlap are seen in the attenuation maps reconstructed with no or a small penalty parameter (Fig. 5.7(a)(b)). They are decreased by using a stronger penalty parameter (Fig. 5.7(d)). Chapter 5. Simulation Experiments 68 Least Squares The objective functions for the least squares (LS), weighted least squares (WLS) and penalized weighted least squares (PWLS), as defined in Sec. 4.2.3, were minimized by the conjugate gradient algorithm (CG) [77]. The maximum number of C G iterations was 100 unless the algorithm converged earlier. The algorithm was assumed to be at its minimum if the difference in objective function between iterations was less than 10 - 1 5 . Figs. 5.8 and 5.9 show the changes of the accuracy and the precision during the iterations. The figures showed that about 20 iterations of conjugate gradient method were enough to minimize the least squares objective function. Further iterations produce only minor changes in the accuracy and precision in the reconstructed images. For non-regularized objectives (LS and WLS), the accuracy and precision of the reconstruction decreased with increasing number of iterations. Again, reconstructed images with no penalty functions revealed some ring artifacts due to penumbra overlap on the detector. Chapter 5. Simulation Experiments 69 Figure 5.8: Local accuracy and % uncertainty of ROI-1 plotted vs. number of iterations. Rows a), b), and c) correspond to the total activities in the CLS source equal to 1.4, 4.8, and 11.8 GBq respectively. The value of penalty parameter is given (in the brackets) in the legend. Chapter 5. Simulation Experiments 70 ROI 2 0.01 H , , , , , , , , , , , 1 OH , { , 1 . , . , . , , 1 0 5 10 15 20 25 30 0 5 10 15 20 25 30 ITERATION ITERATION Figure 5.9: Accuracy and % uncertainty of ROI-2 plotted vs. iteration number. Rows a), b), and c) correspond to different activities in the CLS source equal to 1.4, 4.8, and 11.8 GBq respectively. The value of penalty parameter is given (in the brackets) in the legend. Chapter 5. Simulation Experiments 7\ (a) (b) Figure 5.10: Central slices of the attenuation maps reconstructed by the least squares method with a total activity in the C L S source of 4.8 G B q . (a) 10 iterations of LS , (b) 10 iterations of W L S , (c) 20 iterations of P W L S with the penalty parameter a = 5, and (d) 20 iterations of P W L S with penalty parameter a = 10. Chapter 5. Simulation Experiments 72 Comparison of the Reconstruction Methods Table 5.1: Comparison of the different methods of the reconstruction for CLS Reconstruction Method ROI 1 ROI 2 Global Accuracy Accuracy Uncertainty Accuracy Uncertainty FBP cutoff 0.4 0.008 1.0 0.092 12.3 1.47 FBP cutoff 0.7 0.044 5.6 0.058 7.9 1.35 LS 10th iter. 0.065 7.6 0.056 7.3 1.46 WLS 10th iter. 0.064 7.6 0.053 6.9 1.42 PWLS2.0 20th iter. 0.040 5.1 0.058 6.6 1.32 PWLS5.0 20th iter. 0.031 4.1 0.062 7.5 1.35 PWLS10.0 20th iter. 0.030 3.5 0.067 8.8 1.37 ML 20th iter. 0.061 7.2 0.055 7.2 1.60 PML0.4 20th iter. 0.053 6.4 0.054 6.8 1.40 PML1.0 20th iter. 0.045 5.5 0.055 6.6 1.37 PML2.0 20th iter. 0.036 4.5 0.060 6.7 1.36 The results of the reconstruction tests are summarized in Table 5.1. Although the PWLS gave the best results, the differences between the methods, as presented in Table 5.1 and in the figures showing the reconstructed attenuation maps (Figs. 5.4, 5.7, and 5.10 on pages 64, 67, and 71) , were small. All the techniques could reduce the ring artifacts created by penumbra overlaps by using more smoothing. The PWLS method was better suited for the reconstruction of the CLS attenuation maps than M L because it converged faster than M L and the weighting of the data made the penumbra overlap less important. The overlapped regions contributed little to the objective function. The quality of the attenuation maps obtained by FBP was comparable to that of other methods, but the time required for calculations was substantially shorter. Therefore, unless specified otherwise, the FBP was selected as a method of choice for attenuation map reconstruction in the remaining part of the chapter. Unless otherwise specified, the Butterworth filter, order 5, cutoff 0.5, was used. Chapter 5. Simulation Experiments 73 5.3.3 Number of Divisions One of the most important parameters in map reconstruction for the CLS is the number of divisions used in rebinning of the data. There is a trade off between resolution and noise when the transmission beams are divided. If there are more divisions then more information is available for reconstruction, and approximation done during rebinning (Eq. 4.10) is better, but there are fewer counts in each new bin, and therefore the statistical error increases. To find the optimal number of divisions, we tested attenuation maps obtained with fan beams divided into 1 to 6 divisions. The accuracy and precision were compared for the ROIs, and the images were inspected visually. The reconstructions were again done using both FBP and the iterative methods. Since the results of the tests for all reconstruction methods were similar, only the FBP reconstructions are discussed here. Table 5.2: Accuracy in the ROIs for FBP reconstructed maps for low and medium activity in CLS source for different number of divisions. Divisions ROI 1 ROI 2 Global Accuracy low medium low medium low medium 1 0.122 0.107 0.320 0.332 2.35 2.35 2 0.055 0.017 0.089 0.092 1.70 1.72 3 0.061 0.017 0.104 0.089 1.42 1.48 4 0.072 0.047 0.110 0.054 1.29 1.50 5 0.085 0.066 0.089 0.040 1.31 1.68 6 0.138 0.096 0.106 0.056 1.40 1.95 The global accuracies shown in Table 5.2 and visual inspection of the maps show that there is a big improvement when using 3 divisions rather than 1 or 2 divisions. Changes for a higher number of divisions (4-5) are much smaller. We can conclude that 3 is the smallest number of divisions needed for reconstruction. It is also optimal because it is less vulnerable to changes in the penumbra overlap than higher number of divisions, and Chapter 5. Simulation Experiments 74 the results are less sensitive to the noise in the data. (1) (2) (3) Figure 5.11: The attenuation maps obtained from data rebinned into different number of divisions and reconstructed by FBP for CLS medium activity. The number at each figure indicates the number of divisions used in data processing. 5 . 3 . 4 Beam Profile, Thickness of the Lines One of the main features of CLS is its non-uniform activity distribution. The central lines have the highest activity (Fig. 5.12(b)(c)). This section will discuss the advantage of using a non-uniform activity distribution rather than a uniform (Fig. 5.12(a)) activity distribution in the transmission source system. The diameter of the line source will also be discussed here. The larger the diameter of each particular line, the larger the penumbra created on the detector, and the larger the area being illuminated by the source. To decrease the penumbra effect stricter collimation of the source has to be used, but fewer counts will be detected (Fig. 5.12(c)). Chapter 5. Simulation Experiments 75 (a) (b) (c) Figure 5.12: Profiles through the source activity on the detector, (a) 1 m m line sources with uniform activity distribution across the sources, (b) 1 m m and (c) 3 m m diameter line sources with non-uniform activity distribution. As expected, and as is shown in Table 5.3, there was an advantage in using an ex-ponential source activity distribution over the uniform distribution. W i t h the same total activity in the source system the exponential shape gave more accurate attenuation maps. The reconstructed attenuation maps are shown in F i g 5.13, but the differences between them cannot be clearly seen. Chapter 5. Simulation Experiments 76 Table 5.3: The values of accuracy and precision (expressed in % uncertainty) in the ROIs and the global accuracy for different beam profiles, reconstructed by F B P with total medium activity in the C L S system. Profile Type R 0 I 1 R O I 2 Global Accuracy Accuracy % Uncertainty Accuracy % Uncertainty Exponential 1mm 0.017 2.14 0.089 11.27 1.42 Uniform 1mm 0.072 4.14 0.100 11.64 1.47 Exponential 3mm 0.038 4.66 0.092 11.38 1.46 (a) (b) (c) • • • p^pp CiQ^Jp1 f j^^ p Figure 5.13: Attenuation maps acquired using C L S with different beam profiles: (a) exponential with 1 m m , (b) uniform with 1 m m , and (c) exponential with 3 m m in diameter line sources. Chapter 5. Simulation Experiments 77 5 .3 .5 Spatial Resolution - CLS In order to investigate the resolution of the attenuation maps obtained with the C L S system, an object consisting of 27 attenuating rods (/i=0.30 c m - 1 ) each 1 cm in diameter, spaced by 10 cm was simulated. The transmission scan was done using A T P G with the same protocol as before (Sec. 5.2.2). The attenuation maps were reconstructed using F B P . The results of this simulation are shown in F i g . 5.14. Figure 5.14: The reconstructed attenuation map of the set of 27 bars. The image is 256x256 pixels with pixel size of 0.24 cm (the area shown is 61.4x61.4 cm). The outmost bars displayed truncation artifacts. There was a small difference in the the resolution between the bars at different positions due to the collimator blurring effect; the reconstructed maximum values of fx for each bar were very close to each other. The reconstructed value of u- for the most central bar was smaller than for the others. This effect was caused by the particular position of this bar which in , al l projections, was exactly at the edge of the central fan beam. Thus, for all the projections, the photon flux through that bar was weak which resulted in a smaller reconstructed fi, as shown in Chapter 5. Simulation Experiments 78 Fig. 5.14. This effect is closely related to the approximation which was made in Chapter 4 in Eq.4.10. Such an artifact is related to the small size of the object and it will not be present when extended objects are imaged. Chapter 5. Simulation Experiments 79 5.4 R e s u l t s - S i m u l a t i o n s o f M L A 5.4.1 Geometry of the Simulated MLA Source A set of 20 non-uniformly spaced line sources was used for the simulation of the MLA system. The distances between the lines were chosen to simulate an approximately smooth beam profile on the detector. The goal was to produce a beam profile which corresponds to the shape of attenuation of a water cylinder 40 cm in diameter (Fig. 5.15). Such an activity distribution should create a signal to noise ratio that is approximately constant across the detector for the majority of patient studies. This profile can be achieved when using pairs of line sources placed symmetrically about the centre. Each pair of lines was weaker than its more interior neighbouring pair by a constant factor. The factor is related to the decay of the source over a fixed time interval, say six months. The variations in the spacing of the lines were chosen so that the idealized activity profile could be approximated. The scheme permits a regular replenishment of the source (e.g. every six months) by disposing of the weakest (outmost) pair, introducing a new pair of full-strength lines in the centre, and moving the other lines outward. The distances between the lines are shown in Fig. 5.15. Chapter 5. Simulation Experiments 80 (a) (b) fc) -40 -30 -20 -10 0 10 20 30 40 Position on the detector [cm] Figure 5.15: Activity distribution for the M L A . (a) the distances between the lines are shown in cm. Comparison of the beam profiles for a large 60 cm water cylinder obtained by a simple calculation (b) and using simulation of a blank projection by A T P G is presented in (c). The calculated line is produced by; Profile = Norm • exp(0.15 • d) where Norm is a factor normalizing Profile to the A T P G simulation. Chapter 5. Simulation Experiments 81 5.4.2 Reconstruction Methods for M L A F i l t e r e d B a c k p r o j e c t i o n The attenuation maps generated with the M L A system were reconstructed using FBP with the Butterworth filter (order 5) with the cutoff varying from 0.3 to 1.0 of the Nyquist frequency. Results of these reconstructions for three different activities in the source are shown in Fig. 5.16. As for the CLS system, increase in the total activity in the source improved the accuracy and the precision of the reconstructed image. There was a trade-off between accuracies in the two ROIs; smoother image had im-proved accuracy in the centre, but worse at the edge. M a x i m u m L i k e l i h o o d The transmission maximum log-likelihood (ML) and penalized log-likelihood (PML) objective functions (Sec. 4.2.3) were maximized by the conjugate gradient method [77]. The changes of the accuracy and precision (% uncertainty) vs. number of iterations are shown in Fig. 5.18. The reconstructed images (20th iteration) are presented in Fig. 5.20. The results were similar to the results for CLS reconstruction. The accuracy and precision of non-regularized ML quickly decreased with increasing number of iterations. This effect was removed by use of a penalized objective function. Chapter 5. Simulation Experiments 82 Figure 5.16: Accuracy and precision of ROIs in the reconstructed image plotted vs. cutoff frequency. The top row corresponds to low activity in the source, the middle one to medium, and the bottom one to high activity in the M L A system. Chapter 5. Simulation Experiments 83 Figure 5.17: The central slices of the attenuation maps of the phantom reconstructed by FBP with the Butterworth filter (order 5) and different cutoff frequencies equal to 0.45, 0.65, 0.85, and 1.00 for medium activity of M L A system. Chapter 5. Simulation Experiments 84 ROI 1 I T E R A T I O N I T E R A T I O N Figure 5.18: Accuracy and precision plotted vs. iteration number for the maximum likelihood algorithm, for a total activity in the M L A source equal to 2.8 GBq. The value of penalty parameter is given in the legend. Chapter 5. Simulation Experiments 85 N o p e n a l t y o c = 0 . 2 Figure 5.19: The value of log-likelihood vs. the number of iterations for maximum like-lihood reconstruction. The graph is for a total activity in the M L A source equal to 2.8 GBq. The lines correspond to different penalty parameters given in the legend. NB: The plotted function is not an M L objective function but the value of the log-likelihood. Chapter 5. Simulation Experiments 86 (c) (d) Figure 5.20: Central slices of the attenuation maps of the phantom reconstructed by the maximum likelihood method with (a) 20 iterations of non-penalized (ML) (b) 20 iterations of P M L with penalty parameter a = 0.4, (c) 20 iterations of P M L with a = 1, and (d) 20 iterations of P M L with a = 2 for medium activity of M L A source. Chapter 5. Simulation Experiments 87 Least Squares The least squares (LS), weighted least squares (WLS), and penalized weighted least squares (PWLS) methods with up to 100 iterations and with different penalty parameters were investigated for the M L A system. Figures 5.21 and 5.22 show the values of accuracy and precision for the LS recon-struction. Again the accuracy and the precision of the non-regularized methods (LS and WLS) decreased with increasing number of iterations. All the penalized reconstructions produced much better results and found the minimum of the objective function after about 20 iterations. Further iterations produced only minor changes in accuracy and precision. Chapter 5. Simulation Experiments 88 ROI 1 c) .-•••rsiiM" . • W ~ A A A A A A A A A A A A • • • • • • • • • t 1 5 • 0> —•—LS — • — WLS —A— PWLS(2.) — • — PWLS(5.) — • — PWLS(10.) t » t » £ » j i » $ W W W WW www 10 1S 20 25 30 15 20 25 30 ITERATION ITERATION Figure 5.21: Accuracy and precision for ROI-1 vs. the number of iterations. Rows (a), (b), and (c) correspond to activities in the M L A source of 1.4 GBq, 2.8 GBq, and 9.3 GBq, respectively. The value of penalty parameter is given (in the brackets) in the legend. Chapter 5. Simulation Experiments 89 a) b) ROI 2 20 25 - . • • • • • • • { • • • • • • l l l l i i 20 25 15 20 25 30 C) -T 1 r-- • — LS - • — WLS -A— PWLS(2.) —T— PWLS(5.) —•—PWLS(10.) I T E R A T I O N S I T E R A T I O N S F i g u r e 5 .22: A c c u r a c y a n d p r e c i s i o n f o r R O I - 2 vs. t h e n u m b e r o f i t e r a t i o n s . R o w s a ) , b ) , a n d c) c o r r e s p o n d t o a c t i v i t i e s i n t h e M L A s o u r c e o f 1.4 G B q , 2 .8 G B q , a n d 9 .3 G B q , r e s p e c t i v e l y . T h e v a l u e o f p e n a l t y p a r a m e t e r i s g i v e n ( i n t h e b r a c k e t s ) i n t h e l e g e n d . Chapter 5. Simulation Experiments 90 (a) (b) * ' i \ft^ f IUI, V»M vrr\T -65.45 J i f mm « i 33.53 ^ • (c) (d) Figure 5.23: Central slices of the attenuation maps of the phantom reconstructed by the least squares method with a) 10 iterations of least squares, b) 10 iterations of WLS, c) 20 iterations of PWLS with penalty parameter a = 5, and d) 20 iterations of PWLS with a = 10 for medium activity of the M L A source. Chapter 5. Simulation Experiments 91 C o m p a r i s o n o f the R e c o n s t r u c t i o n M e t h o d s Table 5.4: Comparison of the different methods of reconstruction for M L A with 2.8 GBq total activity in the source. Reconstruction Method ROI 1 ROI 2 Global Accuracy Accuracy Uncertainty Accuracy Uncertainty FBP cutoff 0.4 0.031 8.1 0.063 8.1 1.57 FBP cutoff 0.7 0.067 11.1 0.067 8.2 1.75 LS 10th iter. 0.118 14.5 0.107 12.7 1.94 WLS 10th iter. 0.055 5.7 0.070 8.6 1.75 PWLS2.0 20th iter. 0.035 4.4 0.063 8.0 1.60 PWLS5.0 20th iter. 0.021 2.3 0.061 7.8 1.60 PWLS10.0 20th iter. 0.009 1.2 0.065 8.1 1.68 M L 10th iter. 0.032 3.8 0.072 7.1 1.68 PML0.4 20th iter. 0.037 4.5 0.066 8.1 1.63 PML1.0 20th iter. 0.020 2.3 0.061 7.8 1.62 PML2.0 20th iter. 0.018 1.3 0.062 7.8 1.66 Table 5.4 shows that basically there is no difference between the two iterative methods. Both methods give good attenuation maps after about 20 iterations. FBP gave slightly worse numerical analysis results, and also the attenuation maps were not as good quality as those obtained by iterative methods. This can be seen by considering the results in the Table 5.4 and by visual inspection of the reconstructed attenuation maps (Figs. 5.20, 5.17, and 5.23 on pages 86, 83, and 90). 5.4.3 M L A B e a m P r o f i l e One of the advantages of the multiple line transmission source is that the beam profile can be tailored to the attenuation of a human body. As was for the case of CLS system, a simulation was performed to check the advantage of using a non-uniform activity distribution in the M L A source. Two levels of activity Chapter 5. Simulation Experiments 92 were used: the low corresponding to 1.4 GBq, and the high corresponding to 9.3 GBq of the total activity in the source. For these two activities, simulations were performed for uniform and non-uniform activity distribution. The results of these simulations, presented in the Table 5.5, clearly show the advantage of using the non-uniform activity distribution. Table 5.5: The value of accuracy and precision in the ROIs for different beam profiles, reconstructed by FBP with the Butterworth filter (order 5, cutoff 0.4) for the M L A system. Source Profile ROI 1 ROI 2 Global Activity Type Accuracy %Uncertainty Accuracy %Uncertainty Accuracy Low Exponential 0.032 6.35 0.058 7.62 1.59 Uniform 0.037 5.69 0.094 11.01 1.88 High Exponential 0.025 4.50 0.058 7.44 1.54 Uniform 0.031 4.86 0.081 11.63 1.72 5.4.4 S p a t i a l R e s o l u t i o n - M L A In order to investigate the resolution of the attenuation maps obtained with the M L A system, a set of 27 bars, 1 cm in diameter, was modelled. The acquisition protocol was the same as for previous simulations, so the data was acquired by a 180° camera rotation. Reconstruction was done by FBP. The outermost bars displayed truncation artifacts and there was a difference in the resolution due to the collimator blurring effect. The resolution decreased when the average distance between the object and the detector increased. Due to the 180° acquisition, this effect was much stronger for M L A than it was for CLS, where 360° acquisition was used. Again, as it was for the CLS, the middle bar was more blurred the other bars. This was due to the special positioning of the bar which for all the projections is exactly between Chapter 5. Simulation Experiments 90° 270° Figure 5.24: Reconstructed image of 27 bars, each 1 cm in diameter. The image contains 128x128 pixels of 0.48 cm. Acquisition was over 180°, from 0° to 180°. two fan beams, which lowered the photon flux going through the bar. It was a very similar to that for CLS which was discussed in (Sec. 5.3.5). It could be seen only when small objects were being scanned. 5.5 S u m m a r y a n d C o n c l u s i o n s This chapter covered the simulation experiments performed for both multiple line source systems. All the simulations were generated using the analytical transmission projection generator. Throughout the experiments, the same phantom was used as a attenuating object; it was a torso phantom with a lung insert, and three small non-attenuating spheres placed in the centre of this 3D phantom. Additionally, a set of 27 bars was used to investigate resolution issues. For both systems, CLS and M L A , 64 Chapter 5. Simulation Experiments 94 projections were acquired, the difference being that for the CLS they were collected over 360° of camera rotation, whereas for the M L A only over 180°. Three levels of activity were simulated for both sources; for CLS they corresponded to 1.4, 4.8, and 11.8 GBq and for M L A to 1.4, 2.8, and 9.3 GBq. The transmission data were acquired into a 128x128 matrix and then reconstructed using 2D reconstruction. The central slice was quantitatively analyzed using accuracy and precision parameters calculated over the 2 ROIs. In addition, the global accuracy over the whole slice was evaluated, and reconstructed images of the attenuation maps were visually inspected. It was found that both transmission sources produced attenuation maps without ser-ious artifacts. The small ring artifact due to penumbra overlap for CLS could easily be removed by regularization techniques. Iterative techniques gave better quality of the re-constructed maps than FBP. The penalized weighted least square (PWLS) method gave the best results. Although the iterative methods were superior to FBP, the reconstruction times needed for iterative methods were much higher (about 50-100 times more) than the time required for the FBP reconstruction It was demonstrated that the methods of preprocessing of the CLS transmission data before reconstruction work well and that the optimal number of divisions is 3 (Sec. 5.3.3). The thickness of the line source is an important issue for the CLS; the thicker the line source, the larger the penumbra and tighter the source collimation must be used in order to prevent large overlap of the fan beams. Stricter collimation results in loss of count density in projections of the CLS system (Fig. 5.12(c)), which leads to poorer attenuation maps (Table 5.3). For both sources it was shown that a non-uniform activity distribution in the source is better than uniform distribution (Sec. 5.3.5, 5.4.4). Spatial resolution of the attenuation maps and accuracy of reconstructed yu's were investigated for both sources. For CLS, the Chapter 5. Simulation Experiments 95 reconstructed bars had almost uniform maximum /i-value over the reconstruction area. The M L A reconstruction gave different values of u. for different bar positions. This was because the M L A acquisition was not symmetrical (180° rotation), so the bars which were closer to the detector during the scan had better resolution than the others. This effect is caused by collimator blurring similar to emission S P E C T studies where the resolution decreases with the increase of the distance from the detector [78], and can only be seen for small objects. To conclude this chapter it must be said that both attenuation sources were proven to give good quality attenuation maps. The reconstruction techniques for CLS and M L A dis-cussed in Chapter 4 result in accurate reconstruction of attenuation maps of the simulated phantom. Chapter 6 Experiments - Generation of Attenuation Maps 6 . 1 G o a l o f t h e E x p e r i m e n t s The final step of the investigation of new transmission sources was to perform experi-ments with prototypes of the transmission devices. All the experiments were done in the Nuclear Medicine Department of Vancouver Hospital and Health Science Centre. The goal of the experiments was to confirm the feasibility of the multiple lines geometry and to demonstrate that transmission sources using systems of multiple line can indeed be used for experimental determination of patient specific attenuation maps. 6 . 2 M e t h o d s 6.2.1 Prototype of the CLS The prototype of the CLS was designed to be used with a Sophy 2-head DST cam-era. The transmission device was built at T R I U M F (Tri-University Meson Facility) in 96 Chapter 6. Experiments - Generation of Attenuation Maps 97 Vancouver. G e o m e t r y o f T h e C L S P r o t o t y p e The CLS prototype was designed to be mounted on one of the 2 heads of the DST Sophy detector which was employed solely as the support for the transmission source while the other was used to acquire the transmission data. The source consisted of 10 line sources separated by 4 cm distances. The distance between the source and the detector was the same for all the experimental trials and equal to 65 cm. The CLS was collimated only in one direction by collimators placed parallel to the lines. The radioactive line sources were 1.5 mm in diameter, 16 cm in length, and filled with 1 5 3 G d which decays emitting gamma rays with half live 242 days. The activities in the sources were 200 MBq (5.4 mCi) for the hottest, central lines and 52 MBq (1.58 mCi) for the coldest, peripheral ones. Total activity in the source was 1.16 GBq (31.4 mCi). The ratio of the activity between the neighbouring lines was constant and equal to 0.7. Although the distance between the line sources could have been changed (each line source was mounted separately), for all these experiments the distance was kept the same (4 cm). Fig 6.1 shows an experimental setup of the CLS system built only from 6 line sources (the results presented are from a source system containing 10 lines). D e s i g n o f the L i n e S o u r c e We wanted to have the possibility of changing the parameters of the CLS prototype, such as collimation, position of the sources etc, and thus we decided to built the trans-mission source as a set of separate units. The position and the collimation of each line source could be changed. Figure 6.2 shows the detailed design of each line source along with a photograph. Chapter 6. Experiments - Generation of Attenuation Maps 98 Figure 6.1: Picture of the prototype of CLS. The Sophy detector with 2 heads is shown. The transmission source containing 6 line sources is mounted on one of the heads (the upper one); the other head is used for data acquisition. Between the heads, the Jaszczak Phantom is being scanned. Each line source contained two collimators. The protection collimator prevented the radiation from going to the sides of the line source and the beam collimator formed the fan beam. The collimation of the line source could be altered by changing the distance between the protection and beam collimators (Fig. 6.2). Figure 6.2: The design of line source for the CLS. A detailed scheme of the design is shown in (a) and a photograph of the line source prototype in (b). Chapter 6. Experiments - Generation of Attenuation Maps 100 6.2.2 Design of the M L A The prototype of the M L A was built by Siemens Medical Systems in Chicago. It was designed for use with a one headed Siemens Diacam camera. G e o m e t r y o f the M L A p r o t o t y p e The prototype M L A transmission source was mounted opposite to the detector head in the Siemens Diacam camera. It was mounted on a steel bar, at approximately 100 cm distance from the detector. It contained 20 1 5 3 G d line sources each 16 cm in length and 3 mm in diameter. The activity of each of the hottest lines, the central ones, was 407 MBq (11 mCi), and the coldest, the outermost lines, were of activity 17 MBq (0.46 mCi) each. The total activity in the source was 2.78 GBq (75 mCi). The line sources were placed inside an aluminum box, shielded with lead, with one side open for the source collimator (Fig. 6.4). The size of the box was 59x21x4.5 cm. A photograph of the experimental setup of the M L A is presented in Fig. 6.3. C o l l i m a t i o n o f the M L A In the M L A prototype, the source was collimated by placing a parallel hole collimator on top of the line sources. It collimated the beam in both directions, perpendicular and parallel to the direction of the line sources. The collimation angle was approximately equal to the acceptance angle of the parallel hole collimator which was used (Fig. 6.4). Chapter 6. Experiments - Generation of Attenuation Maps 101 Figure 6.3: Experimental setup for the M L A . The single head Siemens Diacam camera is used. The picture shows data acquisition for a patient heart perfusion scan. The M L A transmission source is attached to the gantry by a steel bar. Opposite the M L A , behind the patient, the detector head can be seen. Figure 6.4: Collimation of the M L A source. A block of a standard high sensitivity parallel hole collimator (a) is placed on top of the aluminum, lead shielded box containing line sources (b) (only 5 are drawn). Chapter 6. Experiments - Generation of Attenuation Maps 102 Figure 6.5: Prototype of the M L A . The line sources are shown on the bottom of the aluminum box. The source collimator used is shown on top of the box. Chapter 6. Experiments - Generation of Attenuation Maps 103 6.2.3 Phantoms Used in Experiments Four different objects were scanned in our experiments. For the CLS only the data for the lung phantom are presented. For M L A , the Jaszczak phantom, the lung phantom, the lung phantom with bags and patients were scanned. The larger variety of objects used for the M L A is due to the better geometrical flexibility of the Siemens Diacam camera. Simply, larger objects like the lung phantom with bags and patients did not fit into the camera on which the CLS prototype was mounted. The Jaszczak phantom was also used for experiments with the early CLS prototype. The early prototype of CLS was built from only 6 line sources and the experiments gave unsatisfactory attenuation maps which are not comparable to those presented in this thesis. The results of these early tests were presented in [52]. The Jaszczak phantom is a 22 cm diameter cylinder filled with water. The lung phantom is a elliptical water cylinder with diameters of 32.0 cm and 23.5 cm. There are 2 low attenuating inserts inside this phantom and a high attenuating, 37 mm diameter bar, imitating lungs and spinal cord respectively. The lung phantom with bags is the lung phantom with attached water bags to increase the size of the object to simulate big, highly attenuating, patients. 6.2.4 Data Acquisition and Reconstruction For the CLS system, 64 projections over 360° camera rotation were taken. The dur-ation of each projection was 40 s. For the M L A , 64 projections of 20 s each were taken over 180° camera rotation. The duration of each projection for M L A system was half that of CLS because the activity in the CLS system was approximately half that used in the M L A system. For both, the data was acquired in a 128x128 matrix, and number of slices was reduced from 128 to 63 by summing counts from adjacent slices (the first and the last Chapter 6. Experiments - Generation of Attenuation Maps 104 slices were discarded). The pixel size was 0.34 cm for the CLS and 0.48 cm for the M L A . The transmission data for both sources were preprocessed as was described in chapter 4 and reconstructed with the penalized weighted least squares (PWLS) algorithm with different penalty parameters. The patient attenuation maps (Sec. 6.3.2) were acquired during a simultaneous trans-mission/emission scan and required correction for emission down-scattered photons. This correction is described in Chapter 7. 6.3 Results 6.3.1 Experimental Results - CLS Transmission Scans The first row of Fig. 6.6 shows reconstructed bolts in the lid of the phantom, each bolt being 2 cm in diameter. These are displayed to demonstrate the resolution of the attenuation map. The other two rows show reconstructed attenuation maps of the lung phantom; the middle row shows the slice from the lung region and the bottom row shows the uniform region of the phantom (containing only spinal cord insert). All the maps presented in Fig. 6.6 have good quality. In the least smoothed images (first column) we can see ring artifacts due to penumbra overlap (this effect was also seen in simulations (Chapter 5), but as expected, they are decreased by using stricter smoothing. Since the phantom was filled with water we expected that quantitatively the resulting attenuation coefficient for the uniform regions of the phantom would be equal to the linear attenuation coefficient of water for 100 keV; the radiation energy of 1 5 3 G d . That coefficient according to the tables [79] is equal to 0.171 c m - 1 . Table 6.2 shows the summary of our assessment of the value of fi in the 10x10 pixel ROIs defined as in Fig. 6.6 and measured Chapter 6. Experiments - Generation of Attenuation Maps 105 in the lung and the uniform parts of the phantom. Quantitatively, the value of \x measured in the ROI agreed with the expected values to within 2-3%. Table 6.1: Values of fi [cm - 1] and their standard. deviations found in ROI for the lung phantom scanned using transmission CLS system. Data was reconstructed with 20 iter-ations of PWLS with different values of the penalty parameter a. Region = 2 a = = 5 a = : 10 h a fi cr h a Lung .1709 .0109 .1763 .0070 .1813 .0049 Uniform .1786 .0085 .1719 .0047 .1747 .0034 Chapter 6. Experiments - Generation of Attenuation Maps 106 Figure 6.6: Slices from the reconstructed attenuation maps of the lung phantom obtained by the CLS prototype. The rows correspond to the bolts on the lid of the phantom, a slice from the lung region, and a slice from the uniform region, respectively. In each row the columns correspond to the 20 iterations of PWLS reconstructions with different values of the penalty parameter (a = 2, 5, and 10). The ROIs from which the values of \i were calculated are marked in the first column. Chapter 6. Experiments - Generation of Attenuation Maps 107 6.3.2 Experimental Results - M L A Transmission Scans Using the prototype of the M L A transmission source described previously, several scanned. M e t a l B a r s The problem of spatial resolution was addressed already in Sec. 5.4.4 using simulated data. It was shown that although the reconstructed thin bars did not show geometrical artifacts, the resolution changes depending on the position of the bar. We decided to investigate this issue further using experimental data. Three long aluminum bars each 2 cm in diameter were used. They were positioned parallel to the axis of rotation, and as much off-centre as possible. The distances of the bars from the axis of rotation were 20.0, 20.5, and 26.5 cm (Fig. 6.7). tests were performed. The tests varied in terms of the complexity of the object being 9 Bars 20cm 26.5cm Axis of rotation Figure 6.7: Position of metal bars in experiment with M L A system. 128 projections over 360° were acquired. From the 360° data set, four sets of 64 projections each taken over 180° data were extracted, starting at 0°, 90°, 180°, and 270°, Chapter 6. Experiments - Generation of Attenuation Maps 108 and ending at 180°, 270°, 360°, and 450°, respectively. The resulting 4 sinograms are shown in Fig. 6.8. In this figure the discontinuity in the sinogram of the largest radius is seen. This discontinuity is caused by close distance between the rod and transmission source for some projections. For the projection just before the discontinuity the rod is imaged by one line source, in next projection it moves to "dead area" between 2 fan beams (no image of the rod is visible on the detector), and in next projection the rod is imaged by different line source, and the image of the rod on the detector is positioned discontineusly. This effect is only present when object is positioned very close to the transmission source and never will occur in normal data acquisition. Other two objects which are just 6 cm. closer to the centre do not shown this discontinuity. Figure 6.8: The four sinograms of log-ratio for 3 aluminum bars. Each sinogram shows 64 projections of 180° data starting at 0°,90°,180° and 270° (see text). The bars were reconstructed by PWLS (a = 2.0) for each subset, and then the 4 reconstructed images were added together (Fig. 6.9). As the result an image of 12 bars was produced. Such an image corresponded to the reconstruction of 12 bars from the 180° acquisition taken from 0° to 180° (see the marks in the Fig. 6.9). From now on it Chapter 6. Experiments - Generation of Attenuation Maps 109 will be referred to as a reconstructed image of 12 bars. This procedure was done in order to increase the number of experimental points. 90° 270° Figure 6.9: Reconstructed metal bars with PWLS. The positions 0° and 180° correspond to the first and the last projection. The 4 outermost bars displayed truncation artifacts (Fig. 6.8, 6.9) because the bar positioned the farthest away from the centre of rotation was out of the field of view for some of the projections. These 4 bars were not taken into account. Some ringing artifacts can be seen in the Fig. 6.9; they are due to differences in the stationary blank scan (scan without attenuator) and the transmission scan. Due to the gravitational distortion of the gantry there is a slight shift in the relative positions of the source and detector for views taken at different angles, which will create these rings. These ring artifacts are multiplied in the shown image because 4 images are added together, so the ring artifact is 4 times more pronounced than it is for regular data acquisition. The ring artifacts are not seen when extended objects are being scanned. Chapter 6. Experiments - Generation of Attenuation Maps 110 The other 8 bars differed in the maximum value of due to detector collimator blur-ring. Taking the maximum value of the reconstructed attenuation coefficient for a given bar as a measure of the spatial resolution, we could see that resolution decreased linearly with the average distance of the bar from the detector surface during the data acquis-ition. Fig. 6.10 shows the relation of such measured spatial resolution vs. the average distance from the detector for experiment and for the computer simulations described in Sec. 5.4.4. In this figure the linear relation between the resolution and average distance to the detector can be seen. The slope of the lines is similar because the same parallel hole collimator was used for the experiment and simulation. The lines do not overlap because the value of /J, in simulation was equal to 0.3 c m - 1 , while the value of fi for the metal bars was much higher. Chapter 6. Experiments - Generation of Attenuation Maps 111 Figure 6.10: The relation of the maximum \i reconstructed for a given bar vs. the average distance of that bar from the detector for the M L A . The actual value of y, was different for experiment and simulation. Chapter 6. Experiments - Generation of Attenuation Maps 112 Phantoms Fig. 6.11 shows slices of the reconstructed attenuation maps of the phantoms. These results showed that the method gave good attenuation maps. The maps were free of any major artifacts and, quantitatively, the values of the attenuation coefficient \x agreed with the previous experiment using the CLS system (Table 6.1), and with the expected value of \x for 100 keV photons. Table 6.2: Values of /i in [cm - 1] and their standard deviations found in ROIs for phantoms scanned using the transmission M L A system. The reconstruction was done by 20 itera-tions of PWLS with different values of the penalty parameter a. Phantom OL = 0.5 a = 2.0 a = 5.0 fi cr fi a fi a Jaszczak .1794 .0053 .1737 .0083 .1726 .0116 Lung Phantom .1795 .0032 .1728 .0058 .1737 .0074 Lung Phantom + Bags .1794 .0016 .1703 .0048 .1748 .0054 Chapter 6. Experiments - Generation of Attenuation Maps 113 Figure 6.11: Slices from the reconstructed attenuation maps of the phantoms. The rows correspond to the experiments with the Jaszczak, lung phantom, and lung phantom with bags. In each row the columns are the PWLS reconstructions with different values for the penalty parameter, o = 0.5, 2.0, and 5.0. The white square in column 1 designates the region in which ji is determined. Chapter 6. Experiments - Generation of Attenuation Maps 114 Patients Three patients were selected for the emission heart perfusion study with the simul-taneous transmission scan. The patients were chosen for their size, and will be referred to as small, medium, and large. The weights of these patients were 75, 100, and 140 kg respectively. In Fig. 6.12 we can see that attenuation maps are almost of the same quality regardless the size of the patient. Also Table 6.3 provides similar attenuation coefficients /i and their standard deviations for all three patients. Table 6.3: Values of fi [cm'1] and their standard deviations found in ROIs for patients scanned using transmission M L A system. Data was reconstructed with PWLS with different values for the smoothing parameter a. Patient a = 0.5 a 2.0 a = 5.0 fi a /' a fi a Small .1675 .0104 .1679 .0056 .1677 .0038 Medium .1611 .0147 .1629 .0083 .1635 .0052 Big .1546 .0121 .1711 .0073 .1710 .0055 Figure 6.12: Slices of the attenuation maps of patients taken from the lung region. Rows correspond to different patients: small, medium, and large. Columns correspond to different values of the penalty parameter in PWLS reconstruction, a=0.5, 2.0, and 5.0. Chapter 7 Example of the Application of the Attenuation Maps Attenuation maps themselves are not clinically useful unless they are used, for ex-ample, in attenuation correction1 of S P E C T images. In this chapter we present an ex-ample of the acquisition protocol for simultaneous transmission/emission scanning with the M L A and an example of the application of the acquired attenuation maps in attenu-ation correction. It was shown in Sec. 1.8 that an effect of photon attenuation, if uncorrected, is to produce artifacts in emission S P E C T images. This chapter will demonstrate that atten-uation correction, using the attenuation maps determined by the collimated line sources, removes these artifacts and gives improved emission images. 1 Attenuation maps can also be used in an accurate patient specific scatter correction, or for a cross-talk correction for dual-isotope imaging [80]. 116 Chapter 7. Example of the Application of the Attenuation Maps 117 7 . 1 M e t h o d s 7.1.1 Simultaneous Transmission/Emission Data Acquisition Simultaneous transmission/emission scans were performed. Different isotopes were used for the transmission and emission scans. In the example presented in this chapter 1 5 3 G d was used for transmission and 9 9 m T c for emission. The energies of photons emitted by these isotopes were 100 keV and 140 keV, respectively. 7.1.2 Cross Talk Correction With the simultaneous transmission/emission scan, when the transmission energy is lower than the emission energy, many scattered in the object emission photons are detected in the transmission energy window. Without correction, these photons will contribute to the transmission data, changing the resulting values of the /i-map. In our experiments, as stated, we used energies of 100 keV ( 1 5 3 Gd) for transmission, and 143 keV ( 9 9 m Tc) for emission. A schematic representation of the energy spectra of photons on the detector is shown in Fig. 7.1. In this work we used the window subtraction method for down-scatter correction. This method was employed in estimation of the scatter fraction for planar imaging [81] and later for S P E C T [82, 83]. This method uses counts acquired in additional images corres-ponding to the scatter energy window, scales them by a factor k, and subtracts the scaled images from the primary window (Eq. 7.1,7.2). The scatter data before subtraction can be convolved with a smoothing filter in order to reduce noise. In our case the factor k describes the ratio of the number of scattered emission photons detected in the transmis-sion and scatter energy windows. Although the method is not completely correct since the spatial distribution of the scattered Chapter 7. Example of the Application of the Attenuation Maps 118 Emission photopeak Energy Emission scattered photons Figure 7.1: Schematic of the energy spectrum of the simultaneous transmission/emission acquisition. Spectrum of photons scattered in Nal is not shown. The number of these photons is small compared to number photons scattered in the object. photons is different for the different energies [84, 85], it is fast and easy to implement. The errors introduced are often negligible, especially for noisy data [86]. For the simultaneous acquisition of 9 9 m T c and 1 5 3 G d isotopes, the following energy windows were used in our experiments: window-1 for emission data at 140.0 ± 7.5% keV, window-2 for transmission data at 100 ± 10% keV, and window-3 for scatter data made from the sum of counts in two energy bands 86 ± 5% keV and 116 ± 5% keV. The window positions are shown schematically in Fig. 7.2. The window subtraction scatter correction method can be expressed by the following I, B , I s and B s are the counts in transmission, blank, scatter, and blank scatter energy windows. I c o r r and B c o r r are scatter corrected values of transmission and blank scans. c o r r I - k • ( I s * G) (7.1) g c o r r = B - A: • ( B s * G) (7.2) Chapter 7. Example of the Application of the Attenuation Maps 119 Transmission Emission Energy [keV] Scatter Figure 7.2: Energy window positions and width for the simultaneous transmis-sion/emission ( G d / T c ) M L A transmission source. G is a smoothing kernel. Although during the blank scan, B , there are no emission photons, B has to be corrected for scatter as in E q . 7.2. This is because the scatter energy window picks up non-scattered transmission photons (10% of the total number of photons detected in trans-mission window). This is again due to poor energy resolution and the close positioning of the scatter windows relative to the transmission window (Fig. 7.2). Effectively then, in E q . 7.1 we not only subtract scattered emission photons, but also this small 10% fraction of transmission photons detected in the scatter window. To account for that effect we also had to correct for this small fraction in the blank scan (Eq. 7.2). 7.1.3 Emission Reconstruction There are many techniques for the reconstruction of emission images with attenuation correction using attenuation maps [87, 88, 89, 90, 91, 92]. In this thesis, the maximum likelihood expectation maximization ( M L E M ) [17] and the ordered subset expectation maximization ( O S E M ) [19] methods were used as they are becoming a gold standard in image reconstruction in nuclear medicine. Chapter 7. Example of the Application of the Attenuation Maps 120 In order to use an iterative method, such as M L E M , the attenuated projector /back-projector operation has to be defined. The reason for that has already been discussed in Sec. 4.1.3. The projection operation for parallel beam geometry and the fast implementation of the attenuated projector/backprojector with the exact system matrix (the same as for trans-mission) were used. The algorithm is described in detail in Appendix E . In order to assess the effect of the attenuation correction, the emission images were reconstructed with and without attenuation correction. 7 . 2 R e s u l t s 7.2.1 Estimation of the Ratio Parameter k The ratio, k, of counts between window-2 and window-3 as well as the dependence of that ratio on projections was found. The results for 5 different objects are summarized in Table 7.1. They were very similar for all objects. The standard deviation of k for different projections was very small, indicating that it was reasonable to use the same value of k for all the projections. Table 7.1 also shows that k did not depend on the object being scanned. The values of k were directly calculated from the number of counts detected in window-2 and window-3 for the emission only scan. 7.2.2 Effect of Attenuation Correction The S P E C T reconstruction of a lung phantom was done with 10 iterations of OSEM with 16 subsets. The reconstructed emission image was first filtered with a 3D (27 points) median filter and then reoriented along the short axis of the heart. Examples of a transverse slice of the emission image before reorientation are shown in Fig. 7.3, and reoriented along the short axis of the heart in Fig. 7.4. Chapter 7. Example of the Application of the Attenuation Maps 121 Table 7.1: Values of the ratio k of scatter counts between the transmission and scatter energy windows. The standard deviation, cr, of k between the projections is also presented. Object k cr Lung Phantom 1.121 0.026 lung region 1.121 0.026 Lung Phantom uniform region 1.116 0.022 Big Patient 1.099 0.021 Medium Patient 1.124 0.021 Small Patient 1.113 0.019 Fig. 7.5 shows emission images of the same lung phantom but this time scanned with the cold water bags (the attenuation maps of a similar phantom were already presented in Sec. 6.3.2). The emission reconstruction was carried out using 50 iterations of M L E M . The figure 7.5 shows reconstructed transverse slices of that object with and without attenuation correction. 7.3 D i s c u s s i o n a n d C o n c l u s i o n In this chapter the reconstructions of the SPECT emission images were presented. Reconstruction with and without attenuation correction was performed for two phantoms. The attenuation maps acquired by simultaneous transmission scan were used for the attenuation correction. In the first example (Fig. 7.4(a)) the background in image reconstructed without attenuation correction was not uniform throughout the phantom, and the inferior wall of the heart appeared to contain less activity than the anterior wall (30% less activity). This is a well known effect of attenuation which was described in the introduction (Sec. 1.8). Fig. 7.4(b) shows that the attenuation artifacts were removed by applying the attenuation Chapter 7. Example of the Application of the Attenuation Maps 122 (a) Spinal cord L i d o f t n e h e a r t phantom f j (b) Figure 7.3: Example slices of the lung phantom which is schematically presented in (a), and reconstructed with no attenuation correction (b) and with attenuation correction (c). Reconstruction was performed by OSEM. correction based on attenuation maps acquired using our source. The other example (Fig. 7.5(a)) shows that non-compensated attenuation creates geo-metrical distortions in the object. These artifacts were also described in Sec. 1.8. In Fig. 7.5(a), these distortions appeared as streaks going from the phantom. Again, by applying the attenuation correction, they were removed from the image Fig. 7.5(b). The removal of similar artifacts was also shown for a different phantom in Fig. 1.11 on page 22 in the introduction. Chapter 7. Example of the Application of the Attenuation Maps 123 (a) (b) Figure 7.4: Example slices of emission reconstruction with O S E M of the lung phantom, reoriented along the short axis of the heart, and a profile, (a) with no attenuation correc-tion and (b) with attenuation correction Chapter 7. Example of the Application of the Attenuation Maps 124 Figure 7.5: Emission reconstructions with M L E M of the lung phantom with bags, (a) no attenuation correction and (b) attenuation correction. The contrast of both images was adjusted to show more clearly the streak artifacts. Chapter 8 Conclusions This thesis presents a new geometry for the transmission source to be used in simul-taneous or sequential S P E C T imaging in order to obtain attenuation maps of the scanned object. In this final chapter of the thesis, the work is summarized, advantages of our transmission source are emphasized, and some possible clinical applications are outlined. 8 . 1 S u m m a r y o f t h e W o r k The first chapter of the thesis introduced the fundamentals of nuclear medicine imaging. The photon attenuation effect and the effect that it has on reconstructed S P E C T images were discussed. It was shown that artifacts created by attenuation can be removed by applying an attenuation correction based on object specific attenuation maps. Chapter 2 presented methods of acquiring the attenuation maps using techniques which have already been developed. Short descriptions of these methods were presented and their advantages and disadvantages were discussed. Although a few techniques exist for obtaining attenuation maps, all of them have disadvantages. In Chapter 3 of this thesis, a 125 Chapter 8. Conclusions 126 novel geometry of a transmission source used for obtaining attenuation maps for SPECT imaging was proposed. The new source uses a series of collimated line sources positioned parallel to the axis of rotation of the camera, opposite to the detector. Two versions of this system were discussed in the thesis: (i) collimated line source and (ii) multiple line array. The CLS system uses 10 equally spaced collimated line sources. The fan beams originating from each line source do not overlap on the detector. This geometry requires a 360° rotation of the camera. The M L A system is built from 20 unequally spaced collimated line sources. In this version, the transmission device is placed further away from the detector than it was for the CLS case, and the fan beams coming from the line sources do overlap on the detector. Only a 180° rotation is needed to acquire complete data in this case. The two versions of the system have different geometries and therefore require different approaches for the reconstruction of attenuation maps. These approaches were described in Chapter 4, starting from data preparation for reconstruction (rebinning for CLS and missing counts correction for MLA) and ending with the derivation of the objective func-tions for iterative reconstructions. Filtered backprojection and iterative techniques were tested for the reconstruction of the attenuation maps. For CLS, the usage of FBP required an approximation of the fan beams by parallel strips. The iterative methods used the ex-act geometry of the system, so the fan shapes of the beams were taken into account. The exact geometry for the iterative techniques was encoded in the system matrix w which was precalculated and stored. The FBP and iterative techniques were tested also for the reconstruction of the attenuation maps for the M L A , but in this case no approximation was required. The performance of the different reconstruction methods was evaluated in Chapter 5 using simulation experiments. These tests confirmed that attenuation maps could be obtained with our transmission sources using the reconstruction methods described in Chapter 8. Conclusions 127 Chapter 4. Based on these simulations, the optimal values of some essential parameters of the sources were established, such as the number of line sources, their collimation, and the distances between them. Chapter 6, the experimental part of thesis, presented technical details of the designs of the prototypes. Both sources were tested with phantoms and the conclusions from the simulation part of the thesis were confirmed. Accurate and artifact free attenuation maps were obtained. Patient attenuation maps were obtained for 3 patients using the MLA system. The quality of these maps was also good, regardless of the size of the patients. Finally, an example of the application of multiple line sources in attenuation correction of the emission data was presented in Chapter 7 where the transmission data were acquired simultaneously with the emission scans. The emission images were reconstructed with attenuation correction based on the attenuation maps obtained by the multiple line source transmission system. It was shown in this chapter that by using an attenuation correction the artifacts due to photon attenuation can be removed from the reconstructed image. 8 . 2 A d v a n t a g e s o f M u l t i p l e L i n e s T r a n s m i s s i o n S o u r c e s In this section we would like to attempt to answer a question the reader may be asking. Why are multiple line sources better than any other transmission source? We think that answer to that question is in: T h e " S m a r t " G e o m e t r y With all the other techniques of obtaining attenuation maps as described in Chapter 2, the signal to noise ratio (SNR) is dramatically low for the central part of the detector where attenuation is high compared with the side parts of the detector where attenuation is low, and SNR high. A low SNR may create problems during the reconstruction of Chapter 8. Conclusions 128 attenuation maps due to missing counts. The value of the SNR can be improved by increasing the activity of the source, but this will increase the dose to the patient and technologist and will also increase the maintenance cost of the equipment. By using the multiple line source, this problem can be partially avoided by placing higher activity in central lines where high attenuation is expected and low activity in peripheral lines where frequently there is no or, at most, low attenuation. Such a setup helps to equalize the SNR over the area of the detector and makes much more efficient use of the radioactivity in the transmission source (Sec. 5.3.4, 5.4.3). As the radioactivity decreases the whole set of radioactive sources does not have to be replaced. Only the two central lines are replaced, and all the other sources are moved to the next more outside positions to replenish the transmission source to its full strength. T h e S i m p l i c i t y o f the D e s i g n Unlike some other types of transmission sources, our geometry is very simple. It does not require any special hardware or special collimators. T h e F e a s i b i l i t y o f the T r a n s m i s s i o n S o u r c e s The work included in this thesis proves that our transmission sources do give good attenuation maps which can successfully be used in attenuation correction. 8.3 Clinical Application We think that, because of its advantages, multiple line transmission sources can be used clinically for obtaining attenuation maps. In this thesis it was shown, using the prototype of the M L A system, that patient specific attenuation maps can be found. Fig. 8.1 showed the difference between the patient heart perfusion image when they are corrected and not corrected for attenuation. The M L A system will be used commercially as the transmission source for the Siemens Chapter 8. Conclusions 129 Medical Systems ecam camera. The a-version of the M L A system is currently being tested in a few hospitals in Canada, the USA and Europe. A series of 21 patients were scanned using the a-version of M L A also in our hospital. The acquired transmission and emission data were reconstructed using the methods described in this thesis and will be analyzed in a separate paper. Figure 8.2 shows a prototype for the new Siemens ecam camera which will have two M L A systems, one for each head. Chapter 8. Conclusions 130 (a) (b) Figure 8.1: Slice of emission image of the heart. Reconstruction was performed by OSEM (a) without and (b) with attenuation correction. The presented image is reoriented along short axis of the heart. Figure 8.2: Prototype for the new Siemens 2-head ecam camera with 2 "wings" containing two M L A transmission sources attached to the gantry by a single mounting bar. Chapter 8. Conclusions 132 8.4 F i n a l W o r d The goal of this work was to present a new geometry for a transmission source for SPECT imaging, and to check its performance for obtaining attenuation maps. We have shown that both implementations of the geometry, the collimated line sources and mul-tiple line array, produce good quality attenuation maps. We have developed methods to reconstruct the attenuation maps from the multiple line sources, and we have optimized the geometry of the transmission sources for the best performance. We believe that both implementations can be used clinically to improve S P E C T ima-ging by obtaining attenuation maps which can be used for accurate attenuation correction and as a crucial element of an accurate scatter correction technique. Bibliography H .O . Anger. Scintilation camera with multichannel collimators. J. Nuc. Med., 5:515-531, 1964. H.O. Anger and J . McRae. Transmission scintiphotography. J. Nuc. Med., 9:267-269, 1968. G.N. Ramachandran and A . V . Lakshminarayanan. Three-dimensional reconstruc-tion from radiographs and electron micrographs: application of convolution instead of Fourier transforms. Proc. Nat. Acad. Sci., 68(9):2236-2240, 1971. J .A. Parker. Image Reconstruction in Radiology. C R C Press, Inc., Boston, 1990. R . A . Brooks and G. Di Chiro. Principles of computer assisted tomography (cat) in radiographic and radioisotopic imaging. Phys. Med. Biol, 21(5):689-732, 1976. S.R. Deans. The Radon transform and some of its applications. John Wiley & Sons, Toronto, 1983. D.R. Gilland, R . J . Jaszczak, K . L . Greer, and R .E . Coleman. Quantitative SPECT reconstruction of iodine-123 data. J. Nuc. Med., 32:527-533, 1991. T.F. Budinger, G.T. Gullberg, and R . H . Huesman. Emission computed tomography Image Reconstruction From Projections. Springier-Verlag, New York, 1979. S. Kawata and 0 . Nalicioglu. Constrained iterative reconstruction by the conjugate gradient method. IEEE Trans. Med. Imag., MI-4:65-71, 1985. L. Kaufman. Maximum likelihood, least squares, and penalized least squares for pet. IEEE Trans. Med. Imag., 12(2):200-211, 1993. Jorge Llacer. Tomographic image reconstruction by eigenvector decomposition: Its limitations and areas of applicability. IEEE Trans. Med. Imag., MI-l(l):34-42, 1982. M . F . Smith, C E . Floyd, R.J . Jaszczak, and R .E . Coleman. Reconstruction of S P E C T images using generalized matrix inversion. IEEE Trans. Med. Imag., 11:165— 175, 1992. R. Gordon, R. Bender, and G.T. Herman. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J. Theor. Biol., 29:471-481, 1970. 133 Bibliography 134 [14] R. Gordon. A tutorial on ART. IEEE Trans. Nuc. Sci., NS-21:78-93, 1974. [15] K. Kouris, H. Tuy, A. Lent, and G.T. Herman R . M . Lewitt. Reconstruction from sparsely sampled data by A R T with interpolated rays. IEEE Trans. Med. Imag., MI-1(3):161-165, 1982. [16] P.L. Combettes. The convex feasibility problem in image recovery. In P. Hawkes, ed-itor, Advances in imaging and electron physics, volume 95, pages 155-270. Academic Press, New York, 1996. [17] L .A. Shepp and Y. Vardi. Maximum liklihood reconstruction for emission tomo-graphy. IEEE Trans. Med. Imag., MI-L113-122, 1982. [18] K. Lange and R. Carson. E M reconstruction algorithms for emission and transmis-sion tomography. J. Computer Assisted Tomography, 8:306-316, 1984. [19] H.M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag., 13:601-609, 1994. [20] R.L. Eisner, M.J . Tamas, and K. Cloninger. Normal S P E C T thallium-201 bull's-eye display: gender differences. J. Nuc. Med., 29:1901-1909, 1988. [21] S.H. Manglos, F.D. Thomas, G . M . Gagne, and B.J. Hellwig. Phantom study of breast tissue attenuation in myocardial imaging. J. Nuc. Med., 34:992-996, 1993. [22] S.H. Manglos, R.J. Jaszczak, C E . Floyd, L.J . Hahn, K . L . Greer, and R.E. Cole-man. Nonisotropic attenuation in SPECT: phantom tests of quantitative effects and compensation techniques. J. Nuc. Med., 28(10):1584-1591, 1987. [23] S.H. Manglos, R.J. Jaszczak, C E . Floyd, L.J . Hahn, K . L . Greer, and R.E . Coleman. A quantitative comparision of attenuation-weighted backprojection with multiplicat-ive and iterative postprocessing attenuation compensation in SPECT. IEEE Trans. Med. Imag., 7:128-134, 1988. [24] E . G . Depuey and E.V. Garcia. Optimal specificity of thallium-201 S P E C T through recognition of imaging artifacts. J. Nuc. Med., 30:441-449, 1989. [25] F .L. Datz, G.T. Gullberg, P.E. Christian, G.L. Zeng, Y . L . Hsieh, C H . Tung, S. Val-divia, R. Ahluwalia, C M . Anderson, and K . A . Morton. Roc comparision of simultan-eous transmission emission tomography to routine S P E C T for diagnosis of coronary artery disease on thallium imaging. J. Nuc. Med., 34:71P, 1993. [26] F .L . Datz, G.T. Gullberg, G.L. Zeng, Y . L . Hsieh, C H . Tung, P.E. Christian, A. Welch, and R. Clack. Application of convergent-beam collimation and simultan-eous transmission emission tomography to cardiac single photon emission computed tomography. Sem. Nucl. Med., XXIV:17-34, 1994. Bibliography 135 [27] E.P. Ficaro, J.A. Fessler, P.D. Shreve, J.N. Kritzman, P.A. Rose, and J.R. Corbett. Simultaneous transmission/emission myocardial perfusion tomography, diagnostic accuracy of attenuation-corrected 9 9 mTc-sestamibi single-photon emission computed tomography. Circulation, 93:463-473, 1996. R.J. Fleming. A technique for using C T images in attenuation correction and quan-tification in SPECT. Nucl. Med. Comm., 10:83-97, 1989. B.H Hasegawa, T . F . Lang, and J.K. Brown. Object-specific attenuation correction of S P E C T with correlated dual-energy x-ray C T . IEEE Trans. Nuc. Sci., 40:1242-1252, 1993. T . F . Lang, B.H. Hasegawa, S.C. Liew, J .K. Brown, S.C. Blankespoor, S.M. Reilly, E .L. Gingold, and C E . Cann. Description of a prototype emission-transmission computed tomography imaging system. J. Nuc. Med., 33(10):1881-1887, 1992. M.A. King, B.M.W. Tsui, and T.P. Pan. Attenuation compensation for cardiac single-photon emission computed tomographic imaging: Part 1. impact of attenu-ation and methods of estimating attenuation maps. J. Nuc. Card., 2(6):513-524, 1995. H. Maeda, H. Itoh, and Y. Ishii. Determination of the pleural edge by gamma-ray transmission computed tomography. J. Nuc. Med., 22:815-817, 1981. J.A. Malko, R.L. Van Heertum, G.T. Gullberg, and W.P. Kowalsky. S P E C T liver imaging using an iterative attenuation correction algorithm and an external flood source. J. Nuc. Med., 26:701-705, 1985. D. L. Bailey, B.F. Hutton, and P.J. Walker. Improved S P E C T using simultaneous emission and transmission tomography. J. Nuc. Med., 28:844-851, 1987. K . L . Greer, C C Harris, and R.J. Jaszczak. Transmission computed tomography data acquisition with S P E C T system. J. Nuc. Med. Tech., 15:53-56, 1987. B.M.W. Tsui, G.T. Gullberg, and E.R. Edgerton. Correction of nonuniform attenu-ation in cardiac SPECT imaging. J. Nuc. Med., 30:497-507, 1989. J.R. Gait, S.J. Cullom, and E .V. Garcia. SPECT quantitation: a simplified method of attenuation and scatter correction for cardiac imaging. J. Nuc. Med., 33:2232-2237, 1992. E. C. Frey, B .M.W. Tsui, and J.R. Perry. Simultaneous acquisition of emission and transmission data for improved thallium-201 cardiac S P E C T imaging using a technetium-99m transmission source. J. Nuc. Med., 33:2238-2245, 1992. Bibliography 136 [39] Z.J. Cao and B.M.W. Tsui. Performance characteristics of transmission imaging using uniform sheet source with parallel-hole collimator. Med. Phys., 19:1205-1212, 1992. [40] P. Tan, D. Bailey, S. Meikle, S. Eberl, R. Fulton, and B. Hutton. A scanning line source for simultaneous emission and transmission measurements in SPECT. J. Nuc. Med., 34:1752-1760, 1993. [41] H.T. Morgan, B .G. Thornton, D.C. Shand, J.S. Ray, and P.J. Maniawski. A sim-ultaneous transmission-emission imaging system: description and performance. J. Nuc. Med., 35:P193, 1994. [42] C H . Tung, G.T. Gullberg, G.L. Zeng, P.E. Christian, F .L . Datz, and H.T. Morgan. Non-uniform attenuation correction using simultaneous transmission and emission converging tomography. IEEE Trans. Nuc. Sci., 39:1134-1143, 1992. [43] R.J. Jaszczak, D.R. Gilland, M.W. Hanson, S. Jang, K . L . Greer, and R.E. Cole-man. Fast transmission C T for determining attenuation maps using a collimated line source, rotatable air-copper-lead attenuators and fan-beam collimation. J. Nuc. Med., 34:1577-1586, 1993. [44] B.J. Kemp, F.S. Prato, R.L. Nicholson, and L. Reese. Transmission computed tomo-graphy imaging of the head with S P E C T system and a collimated line source. J. Nuc. Med., 36:328-335, 1995. [45] S.H. Manglos, D.A. Bassano, F.D. Thomas, and Z.D. Grossman. Imaging of the human torso using cone-beam transmission C T implemented on rotating gamma camera. J. Nuc. Med., 33:150-156, 1992. [46] M . Ljungberg, S-V. Strand, N. Rajeevan, and M.A. King. Monte carlo simulation of transmission data for improved thallium-201 cardiac S P E C T imaging using a technetium-99m transmission source. J. Nuc. Med., 41:1577-1584, 1994. [47] D.J.Kadrmas, R.J.Jaszczak, and R.E.Coleman. Fb-tct truncation suppression by combining two shifted scans with extrapolation techniques. In 1994 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, pages 1247-1250, Oc-tober 1994. [48] S.H. Manglos, G . M . Gagne, and D.A. Bassano. Quantitative analysis of image trun-cation in focal-beam ct. Phys. Med. Biol, 38:1443-1457, 1993. [49] S.A. Larsson, S. Kimiaei, and T. Ribbe. Simultaneous S P E C T and C T with shutter controlled radionuclide source and parallel collimator geometry. IEEE Trans. Nuc. Sci., 40:1117-1122, 1993. Bibliography 137 A. Celler and A. Sitek. Transmission SPECT scans using multiple collimated line sources. In 1995 IEEE Nuclear Science Symposium Conference Record, pages 1121— 1125, October 1995. A. Sitek, A. Celler, and R. Harrop. Multiple line sources for S P E C T transmission imaging. J. Nuc. Med., 37:120P, 1996. A. Celler, A. Sitek, and R. Harrop. Reconstruction of multiple line source attenuation maps. In 1996 IEEE Nuclear Science Symposium Conference Record, pages 1420-1424, November 1996. A. Sitek, A. Celler, and R. Harrop. Methods of attenuation map reconstruction for collimated line sources. J. Nucl. Med, 38:216P, 1997. A. Celler, A. Sitek, E . Stoub, D. Lyster, C. Dykstra, D. Worsley, and A. Fung. Development of a multiple line source attenuation array for S P E C T transmission scans. J. Nucl. Med., 38(5):215P, 1997. A. Celler, A. Sitek, and R. Harrop. Reconstruction of multiple line source attenuation maps. IEEE Trans. Nuc. Sci., 44(4):1503-1508, 1997. A. Celler, A. Sitek, E . Stoub, D. Lyster, C. Dykstra, D. Worsley, and A. Fung. Development of a multiple line source attenuation array for S P E C T transmission scans. J. Nucl. Med., 1997. submitted. A. Sitek and A. Celler. Different approaches to reconstruction of multiple lines source attenuation maps. IEEE Trans. Med. Imag., 1997. To be submited. K. Lange, M . Bahn, and R. Little. A theoretical study of some maximum likelihood algorithms for emission and transmission tomography. IEEE Trans. Med. Imag., MI-6:106-114, 1987. J.M.Ollinger. Maximum-likelihood reconstruction of transmission images in emis-sion computed tomography via the E M algorithm. IEEE Trans. Med. Imag., 13:89-101, 1994. J.A. Sorenson and M . E . Phelps. Physics in Nuclear Medicine. W.B. Sauders Com-pany, 1987. J.A. Fessler. Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Trans. Med. Imag., 13(2):290-300, 1994. D.L. Snyder and M.I. Miller. The use of sieves to stabilize images produced with the E M algorithm for emission tomography. IEEE Trans. Med. Imag., 35:3864-3871, 1985. Bibliography 138 [63] D.L. Snyder, M.I. Miller, L .J . Thomas, and D.G. Polite. Noise and edge artifacts in maximum-likelihood reconstruction for emission tomography. IEEE Trans. Med. Imag., 6:228-238, 1987. [64] S. Geman and D.E. McClure. Bayesian image analysis: An application to single photon emission tomography. Proc. Stat. Corny. Sec. of Amer. Stat. Assoc., pages 12-18, 1985. [65] Z. Liang, R.J. Jaszczak, and K. Greer. On Bayesian image reconstruction from projections: Uniform and nonuniform a priori source information. IEEE Trans. Med. Imag., 8:227-335, 1989. [66] J.A. Fessler, N.H. Clinthorne, and W.L. Rogers. Regularized emission image recon-struction using imperfect side information. IEEE Trans. Nucl. Sci., 39:1464-1471, 1992. [67] T. Ffebert and R. Leahy. A generalized E M algorithm for 3-D Bayesian reconstruc-tion from Poisson data using Gibbs priors. IEEE Trans. Med. Imag., 8:194-202, 1989. [68] P.J. Green. Bayesian reconstruction from emission tomography data using a modi-fied E M algorithm. IEEE Trans. Med. Imag., 9:84-93, 1990. [69] I.G. Zubal, M . Lee, A. Rangarajan, C.R. Harrel, and G. Gindi. Bayesian reconstruc-tion of S P E C T images using registered anatomical images as priors. J. Nuc. Med., 33:963-972, 1992. [70] D.S. Lalush and B.M.W. Tsui. A generalized Gibbs prior for a maximum a posteriori-EM S P E C T reconstruction with fast and stable convergence properties. J. Nuc. Med., 33:832-840, 1992. [71] Shih-Chung B. Lo. Strip and line path integrals with a square pixel matrix: A unified theory for computational C T projections. IEEE Trans. Med. Imag., 7:355—363, 1988. [72] E .U. Mumcuoglu, R. Leahy, S.R. Cherry, and Z. Zhou. Fast gradient-based methods for Bayesian reconstruction of transmission and emission P E T images. IEEE Trans. Med. Imag., 13(4):686-701, 1994. [73] A. Chatziioannou and M.Dahlbom. Detailed investigation of transmission and emis-sion data smoothing protocols and their effects on emission images. In 1994 IEEE Nuclear Science Symposium Conference Record, pages 1568-1572, 1994. [74] G.T. Herman. Image Reconstruction from Projections. Academic, New York, 1980. Bibliography 139 [75] R.K. Otnes and L. Enochson. Digital Time Series Analysis. John Wiley, New York, 1972. R.W. Hamming. Digital Filters. Prentice-Hall, Englewood Cliffs,N.J., 1977. W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T . Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, New York, 1989. P.J Early and D.B Sodee. Principles and Practice of Nuclear Medicine. Mosby-Year Book, Inc., St. Louis, 1995. H. Johns and J. Cunningham. The Physics of Radiology. Charles C. Thomas, Springfield, Illinois, fourth edition, 1983. R.G. Wells, A. Celler, and R. Harrop. Experimental validation of an analytical method of calculating S P E C T projection data. IEEE Trans. Nuc. Sci., 44(3):1283-1290, 1997. P. Bloch and T. Sanders. Reduction of the effects of scattered radiation on a sodium iodine imaging system. J. Nuc. Med., 14:610-614, 1973. R.J. Jaszczak, K . L . Greer, C E . Floyd, C C Harris, and R.E . Coleman. Improved SPECT quantification using compensation for scattered photons. J. Nuc. Med., 25(8):893-900, 1984. M . Smith and R. Jaszczak. Generalized dual-energy-window scatter compensation in spatially varying media for SPECT. Phys. Med. Biol., 39:531-546, 1994. C E . Floyd, R.J. Jaszczak, K . L . Greer, and R.E. Coleman. Deconvolution of Compton scatter in SPECT. J. Nuc. Med., 26(4):403-408, 1985. G. Hademenos, M . Ljungberg, M . King, and S. Glick. A Monte Carlo investigation of the dual photopeak window scatter correction method. IEEE Trans. Nuc. Sci., 40:179-185, 1993. C. Floyd, R. Jaszczak, C. Harris, and R. Coleman. Energy and spatial distribution of multiple order Compton scatter in SPECT: a Monte Carlo investigation. Phys. Med. Biol., 29:1217-1230, 1984. S. Bellini, M . Piacentini, C. Cafforio, and F. Rocca. Compensation of tissue absorp-tion in emission tomography. IEEE Trans. ASSP, 27(3):213-218, 1979. G.T. Gullberg and T . F . Budinger. The use of filtering methods to compensate for constant attenuation in single-photon emission computed tomography. IEEE Trans. Biomed. Eng., 28:142-147, 1981. Bibliography 140 [89] E. Tanaka, H. Toyama, and H. Murayama. Convolutional image reconstruction for quantitative single photon emission computed tomography. Phys. Med. Biol., 29:1489-1500, 1984. [90] L .T . Chang. A method for attenuation correction in radionuclide computed tomo-graphy. IEEE Trans. Nuc. Sci., NS-25:638-643, 1978. [91] T. Morozumi, M . Nakajima, K. Ogawa, and S. Yuta. Attenuation correction methods using the information of attenuation distribution for single photon emission C T . Med. Imag. Technoi, 2:20-28, 1988. [92] G.L. Zeng, G.T. Gullberg, B.M.W. Tsui, and J. Tarry. Three-dimensional iterative reconstruction algorithms with attenuation and geometric point response correction. IEEE Trans. Nuc. Sci., NS-38:693-702, 1991. Appendix A Generation of the System Matr ix for the Collimated Line Source This chapter describes the most important numerical part of the iterative reconstruc-tion: the numerical creation of the system matrix w for the collimated line source geometry. The value of the matrix element calculated in the Sec. 4.1.3 is: The challenging task was to calculate a,^ -, which is the area of overlap of pixel i and the fan-beam j. In general, fan-beam j is a subsection of the fan-beam originating from one line source (see Chapter: 4- Methods of Reconstruction). The other terms in Eq. A . l were easy to find. The parameter a was the collimation angle of the fan beam divided by the number of divisions, and dij was the distance from the centre of the pixel to the line source from which the fan beam j originates. 141 Appendix A. Generation of the System Matrix for the Collimated Line Source 142 A . l Calculation of ajj The area a;j is the intersection area of a pixel i and the fan-beam j. To calculate a,-j, the program calculated the parameters of the equations of two lines defining the boundaries of the fan beam y = mi • X + bx and y = m2 • X + 62 (Fig. A.3) and then found the distances from the pixel i to these lines. Taking into account these distances, there were 15 possible relative positions of these lines across the pixel (Fig. A.4). Considering these 15 possible ways and assuming that lines are parallel to their average direction the algorithm found the the area a 8 j . The assumption that lines are parallel to the average direction introduced a negligible error due to the small angle between the two lines. For more details of the calculation please refer to the code. fan beam j y=m 2 X+b2 Figure A.3: Mathematical formulation of the calculation of a t )j A .2 Compression of Sparse Matrix w Since the calculated matrix w was sparse, only the non zero elements were stored in order to save space. The storage element was a structure which stored in its first 2 bytes Appendix A. Generation of the System Matrix for the Collimated Line Source 143 Figure A.4: All the possibilities of the relative position of two parallel lines to the pixel. Each has to be treated differently because the area of the pixel between them has to be calculated from different formulas. the index of the pixel and in the next 4 bytes the value of the matrix element. All the non-zero matrix elements corresponding to the given projection element j (row of matrix w) were stored as a series, with a binary number proceeding the series indicating the number of elements in the series. 1 Number of storage elements for following series corresponding to storage element storage element storage , element i storage element given projection element Next Series Figure A.5: Scheme of matrix Wij compression. Appendix A. Generation of the System Matrix for the Collimated Line Source 144 A . 3 C o d e f o r G e n e r a t i o n o f t h e S y s t e m M a t r i x f o r C L S /* Code f o r c a l c u l a t i o n of the System Matrix f o r CLS */ /* by Arkadiusz Sitek */ /* l a s t modified June 1997 */ #include <stdio.h> ' #include <string.h> #include <math.n> #define OK 1 #define BE 0 #define PI 3.14159265 /* Two kinds of cooridinates */ /* For case 1 */ #define cooxl(A) ( ( f l o a t ) A-0.5-(float)(BINS/2))*BinLenght #define cooyl(A) ((float)-A+0.5+(float)(BINS/2))*BinLenght /* For case 2 */ #define coox2(A) ((float)-A+0.5+(float)(BINS/2))*BinLenght #define cooy2(A) ((float)-A+0.5+(float)(BINS/2))*BinLenght #define min(A,B) ((A)<(B)?A:B) #define max(A,B) ((A)>=(B)?A:B) #define endx(A) *( A) #define endy(A) *(A+1) #define lenght(A) *(A+2) #define direction(A) *(A+3) /* transmission system parameters */ #define PROJECTIONS 64 #define SourceDistance 65.0 #define LINES 10 /* Nr of l i n e sources */ #define Distance 4.0 /* Distance between the l i n e s */ #define DIVISION 5 /* number of d i v i s i o n s */ #define BINS 128 /*MUST BE EVEN?, not necessarely, but better be*/ #define BinLenght 0.34 typedef struct /* The d e f i n i t i o n of the storage element * { /* Note: i t i t s l i g h l t y d i f f r e n t than described i n the */ short i,j; /* th e s i s . Difference i s not important, simply the p o s i t i o n of the p i x e l i s f l o a t w; /* coded in 4 bytes insted of 2 */ } matrix; struct b /* The d e f i n i t i o n of the series */ { /* Number of series i s the same as number of */ matrix m; /* projection elements = PR0JECTI0NS*LINES*DIVISI0N */ struct b *next: } *mat[PROJECTIONS][LINES][DIVISION]; int num[PROJECTIONS][LINES][DIVISION]={0}; /* Number of elements i n the series */ struct b *MakeP(int, i n t , int ); f l o a t GetW(int, i n t , i n t . i n t , i n t ) ; int SaveMatrix( char* ); int SaveInfo( char* ); f l o a t *Vector(double,double,double,double); f l o a t memory[BINS][BINS]={0.0}; int main( i n t argc, char *argv[] ) { char name[30]; i n t i , j , k ; Appendix A. Generation of the System Matrix for the Collimated Line Source 145 if(argc<=l) { e x i t ( l ) ; } ; strcpy(ftname[0],argv[l]); p r i n t f ( " o k \ n " ) ; for(i=0;KPROJECTIONS;i++){ p r i n t f ("new: '/.d \ n " , i ) ; for(j=0;j<LINES;j++){ for(k=0;k<DIVISION;k++){ mat[i] [j] [k] = MakeP( i , j , k ); } } } SaveHatrix(name); pr i n t f ( " o k mat\n"); Savelnfo(name); return( OK ); > struct b* MakeP( i n t pr, i n t l i , i n t d i ) { int i b . j b ; f l o a t r et; struct b *point, *hold, *holdp; int FLAG=1; for(ib=l;ib<=BINS;ib++){ for(jb=i;jb<=BINS;jb++){ i f ( ( r e t = GetW(ib,jb.pr, l i , d i ) ) > 0.0) { if(FLAG){ hold = (struct b*)malloc(sizeof(struct b));FLAG=0;point=hold;holdp = point; > else hold = (struct b*)malloc(sizeof(struct b)); i f ( h o l d == NULL) { p r i n t f ( " a l l o c problems, exiting\n");exit(OK);> holdp->next = hold; hold->m.i = ib; hold->m.j = jb; hold->m.w = r e t ; holdp = hold; num[pr] [ l i ] [di]++; } } } if(IFLAG) hold->next = NULL; i f (num[pr] [ l i ] [di]==0) return(NULL); return(point); > /*=========== ========= = =================^^ f l o a t GetW(int i/*row*/, i n t j/*column*/, i n t pr, i n t l i , i n t d i ) { f l o a t X,Y; /* coordinates of p i x e l */ f l o a t XS.YS, XSC, YSC; /* source coordinates */ double atem,btem,atem2,btem2; /* temporary values of a and b */ double al,bl,a2,b2; /* l i n e s d e s c r i p t i o n */ double t e t a , tetaprim, tetaprimb; f l o a t c o , s i ; f l o a t d l , d i r l , d 2 , d i r 2 , *point; f l o a t T1.T2; /* treshholds */ f l o a t mind.maxd, d, dp, si p , cop, retvalue; in t licz=0; f l o a t diff=0.0,bin_angle, coef; f l o a t b i g , s m a l l , d i f f l , d i f f 2 ; double temp; double dteta; /* ****MARK 1 ******** */ te t a = 2*PI + PI/2. - (double)pr*2.0*PI/(double)PR0JECTI0NS; Appendix A. Generation of the System Matrix for the Collimated Line Source i f ( t e t a >= 2. * PI ) t e t a -= 2.* PI; i f ( t e t a < PI/4. ||(teta > 3*PI/4 && t e t a < 5+PI/4) I I t e t a > 7+PI/4) X = cooxl(j ) ; Y = c o o y l ( i ) ; co = ( f l o a t ) c o s ( t e t a ) ; s i = ( f l o a t ) s i n ( t e t a ) ; /* ********MARK 2********/ /•antisymetric version*/ YSC = ((float)(LINES/2 - l i ) ) * Distance - 0.25*Distance; XSC = -SourceDistance/2.0; temp = 0.5 - (double)di/(double)DIVISION; atem = (temp*Distance)/SourceDistance; btem = YSC - atem * XSC; temp = 0.5 - (double)(di+l)/(double)DIVISION; atem2= (temp*Distance)/SourceDistance; btem2= YSC - atem2 *XSC; if(atem > 100. II btem > 100. || atem2 > 100. II btem2 > 100. ) printf("atem = '/,f btem =*/,f atem2='/,f btem2='/,f \h", atem ,btem, atem2,btem2); /*********MARK 3 ***********/ YS = YSC*co + XSC*si; XS = XSC*co - YSC*si; a l = (atem *co + si)/(co-atem*si); b l = btem/(co - atem*si); a2 = (atem2 *co + si)/(co-atem2*si); b2 = btem2/(co - atem2*si); /* point = Vector(al,bl,(double)X,(double)Y); dl=lenght(point); d i r l = d i r e c t i o n ( p o i n t ) ; point = Vector((float)a2,(float)b2,X,Y); d2 = lenght(point); dir2 = d i r e c t i o n ( p o i n t ) ; mind = min(dl,d2); maxd = max(dl,d2); /* dteta = (0.5*Dist ance -((double)di+0.5)*Distance/(double)DIVISION)/SourceDistance */ dteta =0.0; tetaprim = fmod((teta-dteta),(PI/2.0)); i f ( t e t a p r i m > PI/4.) tetaprim -= PI/2.; tetaprim = fabs(tetaprim); si p = ( f l o a t ) s i n ( t e t a p r i m ) ; cop = (float)cos(tetaprim); TI = BinLenght * 0.707107 * cos(PI/4.0 - tetaprim); T2 = BinLenght * 0.707107 * sin(PI/4.0 - tetaprim); i f ( t e t a p r i m == 0.0) } Appendix A. Generation of the System Matrix for the Collimated Line Source i TI = BinLenght / 2.0; T2 = BinLenght / 2.0; } big = max(dirl,dir2); small = m i n ( d i r l , d i r 2 ) ; if(big<=PI) d i f f = big - small; else if(small<=PI) { d i f f l = big - small; d i f f 2 = small+2.0*PI-big; d i f f = m i n ( d i f f 1 , d i f f 2 ) ; else if(small>PI){ d i f f = big - small; > else { printf("ERROR"); } else X = coox2(i); Y = cooy2(j); t e t a -= PI/2; co = ( f l o a t ) c o s ( t e t a ) ; s i = ( f l o a t ) s i n ( t e t a ) ; /* ****MARK 2 ******** */ /•antisymetric version*/ YSC = ((float)(LINES/2 - l i ) ) * Distance- 0.25*Distance; XSC = -SourceDistance/2.0; temp = 0.5 - (double)di/(double)DIVISION; atem = (temp*Distance)/SourceDistance; btem = YSC - atem * XSC; temp = 0.5 - (double)(di+i)/(double)DIVISION; atem2= (temp*Distance)/SourceDistance; btem2= YSC - atem2 *XSC; if(atem > 100. II btem > 100. II atem2 > 100. II btem2 > 100. ){ p r i n t f ("atem = '/,f btem ='/,f atem2='/,f btem2='/,f \n", atem ,btem, atem2,btem2); > YS = YSC*co + XSC*si; XS = XSC*co - YSC*si; /* ****MARK 3 ******** */ al = (atem *co + si)/(co-atem*si); b l = btem/(co - atem*si); a2 = (atem2 *co + si)/(co-atem2*si); b2 = btem2/(co - atem2*si); /* ****MARK 4 ******** */ point = Vector(al,bl,(double)X,(double)Y); d l = lenght(point); d i r l = d i r e c t i o n ( p o i n t ) ; point = Vector((float)a2,(float)b2,X,Y); d2 = lenght(point); Appendix A. Generation of the System Matrix for the Collimated Line Source 148 dir2 = d i r e c t i o n ( p o i n t ) ; mind = minCdl,d2); maxd = max(dl,d2); /* dteta = (0.5*Distance-((double)di+0.5)*Distance/(double)DIVISI0N)/SourceDistance */ dteta =0.0; /* This 3 l i n e s work as a reminder */ tetaprim = fmod((teta-dteta),(PI/2.0)); i f ( t e t a p r i m > PI/4.) tetaprim -= PI/2.; tetaprim = fabs(tetaprim); sip = ( f l o a t ) s i n ( t e t a p r i m ) ; cop = (float)cos(tetaprim); TI = BinLenght * 0.707107 * cos(PI/4.0 - tetaprim); T2 = BinLenght * 0.707107 * sin(PI/4.0 - tetaprim); i f ( t e t a p r i m == 0.0) •C TI = BinLenght / 2.0; T2 = BinLenght / 2.0; } i f j d i r l < 0. II dir2 < 0.) printf('7.f '/.f \n", d i r l , dir2); > big = max(dirl,dir2); small = min(dirl,dir2); if(big<=PI) d i f f = big - small; else if(small<=PI) •c d i f f l = big - small; d i f f 2 = small+2.0*PI-big; d i f f = m i n ( d i f f l , d i f f 2 ) ; else if(small>PI) d i f f = big - small; > else •C printf("ERROR"); > / * * / /* Considering a l l the possible cases */ /******** MARK 5*******/ i f ( ( ( ( d i f f ) < P I / 2 . 0 ) ) ) { if(mind > TI) {retvalue = 0.0;> else if(mind > T2 ){ d = TI - mind; if(maxd > TI) { retvalue = d*d/(BinLenght*BinLenght*2.0*sip*cop); > else { dp = T i - maxd; retvalue = (d*d-dp*dp)/(BinLenght*BinLenght*2.0*sip*cop); > Appendix A. Generation of the System Matrix for the Collimated Line Source 149 } else { d = T2 - mind: if(maxd > T i ) { retvalue = 0.5 * sip/cop + d/(BinLenght*cop); > else if(maxd > T2) { dp = TI - maxd; retvalue =0.5*sip/cop+d/(BinLenght*cop)-dp*dp/(BinLenght*BinLenght*2.0*sip*cop); else { dp = T2 -maxd; retvalue = (d - dp)/(BinLenght*cop); } } i f ( r e t v a l u e > 0.5) printf("More then 0.5\n"); } else { if(mind > TI) { retvalue = 1.0; else i f ( mind > T2){ d = TI - mind; if(maxd > T l ) { retvalue = 1.0 - d*d/(BinLenght*BinLenght*(float)sin(2.0*tetaprim)); } else { dp = TI - maxd; retvalue = 1.0 - (d*d+dp*dp)/(BinLenght*BinLenght*(float)sin(2.0*tetaprim)); else { d = T2 - mind; i f ( maxd > TI ){ retvalue = 1.0 - 0.5* sip/cop - d /(BinLenght * cop); } else if(maxd >T2) { dp = TI - maxd; retvalue = 1.0 - 0.5* sip/cop - d /(BinLenght * cop) - dp*dp/(BinLenght*BinLenght*(float)sin(2.0*tetaprim)); } else { dp = T2 - maxd; retvalue = 1.0 - sip/cop - (d + dp)/(BinLenght * cop); } > > i f ( r e t v a l u e — 0 . 0 ) return(retvalue); else •C bin_angle = Distance/((float)DIVISION*SourceDistance); coef = (BinLenght*BinLenght/(bin_angle*sqrt((XS-X)*(XS-X)+(YS-Y)*(YS-Y)))); retvalue *= coef; i f ( r e t v a l u e > 1.5) { p r i n t f ("ALERT: i= '/.d ,j= '/.d, '/.. 2f\n" , i , j .retvalue); } return(retvalue); /*************************************************************************************/ int SaveMatrix(char *name) { Appendix A. Generation of the System Matrix for the Collimated Line Source FILE *wsk,*wsk_test; int i , j , k,t; short number; int l i c z n i k ; struct b *m; f l o a t area; wsk = fopen(name, "w"); wsk_test = fopen("test","w"); for(i=0;i<PROJECTIONS:i++){ for(j=0; j<LINES; j++H for(k=0;k<DIVISION:k++){ number = (short)num Ci] [j] DC-Hi = mat [i] [j] [k] ; t fwrite(ftnumber,sizeof(short),l.wsk); licznik=0; area = 0.0; if(number > 0){ do{ fwrite(&(m->m), sizeof(matrix) , 1, wsk); if(j==l&&k==0) memory[m->m.i-l][m->m.j-l] += m->m.w; area+=m->m.w; m = m->next; licznik++; }while(m != NULL); } p r i n t f ( " l i c z n i k '/,d '/,. 2f\n", l i c z n i k , area); i f ( l i c z n i k != number) {. p r i n t f ( " E r r o r occured in saving\nExiting\n"); exit(OK); > } } fwrite(memory, s i z e o f ( f l o a t ) , BINS+BINS, wsk_test); fclose(wsk);fclose(wsk_test); return(OK); / * = = = = = = = = = = = r r = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ========================GE0METRIC FUNCTI0NS=============================== = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = * / f l o a t •Vector(double a, double b, double x, double y) { s t a t i c f l o a t r e t [ 4 ] ; double vx,vy,kat; ret[0] = ( f l o a t ) ( ( y*a - a*b + x) / (a*a+l)); r e t [ i ] = ( f l o a t ) ( ( y*a*a + x*a + b ) / (a*a+l)); / • p r i n t f ("'/.f */.f '/.f '/.f \n",a, b, ret [0] ,ret [1] ); •/ ret [2] = sqrt ( ( x - r e t [ 0 ] ) * ( x - r e t [ 0 ] ) + ( y - r e t [ l ] ) + ( y - r e t [ l ] ) ) ; vx = ret [0] - x; vy = retCi] - y; if(vx!=0.0) kat = (float)atan((double)(vy/vx)); else if(vx==0.0fc&vy==0.0) •C ret[3] =0.0; pr i n t f ( " v x =0 vy = 0\n"); re t u r n ( r e t ) ; Appendix A. Generation of the System Matrix for the Collimated Line Source } else < (vy>0.0)?(kat=PI/2.0):(kat=3.*PI/2.0); ret[3] = kat; re t u r n ( r e t ) ; i f ( v x > 0. && vy >=0.) ret [3] = kat; else i f ( v x > 0. && vy < 0.) ret[3] = kat + 2.*PI; else i f ( v x < 0. && vy >=0.) ret[3] = kat + PI; else i f ( v x < 0. && vy < 0.) ret[3] = kat + PI; else printf("Imposible\n"); r e t u r n ( r e t ) ; /* „/ Appendix B Generation of the System Matr ix for the Multiple Line Array Although the multiple line array geometry is simple, and the system matrix can be generated more easily than for the CLS, we decided to include the description of the method of the generation of the system matrix, because it will also be important for constructing the attenuated projector/backprojector in Appendix E . B . l C a l c u l a t i o n o f t h e M a t r i x E l e m e n t s The calculation of the matrix element was done in a different way than for the CLS since the simple parallel geometry made the problem easier. During the calculation, the pixels of reconstruction area was rotated and the 2 parallel lines which define the beam remained stationary for all projections. After rotation, taking into account the relative position of these lines and the rotated pixel, the area was found. The algorithm took advantage of the symmetry of the pixel. 152 Appendix B. Generation of the System Matrix for the Multiple Line Array 153 Since the M L A system matrix was sparse, it was compressed with the same algorithm as we used for the compression the CLS system matrix. B.2 Ordering of the Elements in the Series The ordering is very important in constructing the attenuated projector for emission reconstruction; for the M L A transmission system it does not play any role. An important part of the calculation of the matrix was that the elements in the series were ordered. They were ordered accroding to the distance of a given pixel to the detector surface. The smaller the distance of the pixel to the detector, the earlier in the series this element is placed. This ordering is shown schematically in Fig B.6. '•••1 2\ Ordering 5>-'l-P 9\ 3 '• , 4 8 '• , 7 '•.6 7 8 5'\ 4 •'. 3 i'Q ''•2 1-, Detector Figure B.6: The ordering of the elements in the series in the system matrix for M L A . The closer the element is to the detector, the earlier in the series it is positioned after the ordering. Appendix B. Generation of the System Matrix for the Multiple Line Array B . 3 C o d e f o r t h e G e n e r a t i o n o f t h e S y s t e m M a t r i x M L A #include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <stdlib.h> #define TRUE 1 #define FALSE 0 #define OK 1 #define BE 0 #define coox(A) ((float)A-0.5-(float)(BINS/2))*BinLenght #define cooy(A) ((float)(BINS/2) - (float)A + 0.5)*BinLenght #define BinPos(i) ((float)(CAMBIN/2 - i - 0.5))*CamBinLenght #define min(A.B) ((A)<(B)?A:B) #define max(A.B) ((A)>=(B)?A:B) /* define direction(A) *(A+3) */ #define PROJECTIONS 64 #define WHOLERot FALSE #define CAMBIN 128 #define CamBinLenght 0.48 #define BINS 128 /*MUST BE EVEN?, not neseserly, but better be*/ #define BinLenght 0.34 #define start_angle 90.0 #define d i r e c t i o n 1. /* 1 - CW -1-CCW */ /* // Start angle i s +90 degree - start_angle Siemens */ #define DISCRETE 0 f l o a t tresh; /* f l o a t PI; */ /* C a r e f u l l , change was made from pr to -pr correct t h i s l a t e r */ typedef struct { unsigned short i ; f l o a t w; } matrix; struct b { matrix m; struct b *next; > *mat[PROJECTIONS][CAMBIN]; struct dan { unsigned short i , j ; f l o a t d i s t ; } sorted[BINS*BINS]; int num[PROJECTIONS][CAMBIN]={0}; struct b *MakeP(int, i n t ); f l o a t GetW(int, i n t , i n t , i n t ) ; Appendix B. Generation of the System Matrix for the Multiple Line Array int SaveMatrix( i n t , char* ); void s o r t ( i n t ) ; f l o a t memory[BINS][BINS]={0.0}; int comp(const void *aa, const void *bb); f l o a t Normalizer[2][CAMBIN]; int main( i n t argc, char *argv[] ) •c char name[30]; int i , j , k ; FILE *wsk_test; if(DISCRETE) tresh =0.5; else tresh = 0.0; /* PI = 2.0 * asin( l . O ) ; */ wsk_test = fopen("Abel:DIAGNOSTIC:testl","w"); /* // if(argc<=l) { e x i t ( l ) ; } ; // p r i n t f ( " G i v e me the name of the output f i l e \n"); //strcpy(ftname[0], a r g v [ l ] ) ; */ p r i n t f ( " G i v e the filename \n"); scanf ("'/,s", name); for(i=0;i<PR0JECTI0NS;i++) s o r t ( i ) ; for(j=0;j<CAMBIN;j++) { Normalizer[0][j] =0.; Normalizer[1] [j] = 0.; } for(j=0;j<CAMBIN;j++) { mat[i][j] = MakeP(i, j ) ; } SaveMatrix(i,name); > ^return( OK ); /*********************************************************************/ i n t comp(const void *aa, const void *bb) { int ret; struct dan *a,*b; a = (struct dan *)aa; b = (struct dan *)bb; if(a->dist < b->dist) ret = 1; else i f ( a - > d i s t > b->dist) ret = -1; else ret=0: re t u r n ( r e t ) ; /*********************************************************************/ void s o r t ( i n t pr) { f l o a t t e t a , s i , co, X, Y; unsigned short i , j ; Appendix B. Generation of the System Matrix for the Multiple Line Array f l o a t dteta; dteta = ( f l o a t ) d i r e c t i o n * (-PI) / (float)PROJECTIONS; if(WHOLERot) t e t a = start.angle * PI / 180. + (f l o a t ) ( p r ) * ( 2 . 0 ) * d t e t a ; else t e t a = start_angle * PI / 180. + ( f l o a t ) ( p r ) * d t e t a ; co = ( f l o a t ) c o s ( - t e t a ) ; /* minus becouse we have to rotate the bins */ s i = ( f l o a t ) s i n ( - t e t a ) ; for(i=l;i<=BINS;i++){ f o r ( j = l ; j<=BINS; j++X X = coox(j) ; Y = cooy(i); sorted[(i-l)*BINS+ j-l].i=i; sorted[(i-l)*BINS+ j-l].j= j ; sorted[(i-l)*BINS+ j-l].dist= X*co - Y * s i ; > } qsort(&sorted[0],BINS*BINS, s i z e o f ( s t r u c t dan), comp); } struct b* MakeP( i n t pr, int b i ) •c int ib; f l o a t ret; struct b *point, *hold, *holdp; int FLAG=1; for(ib=l;ib<=BINS*BINS;ib++){ ret = GetW(sorted[ib-l].j/*col*/,sorted[ib-l].i/*row*/,pr,bi); Normalizer[0] [bi] += re t ; i f ( r e t > tresh*BinLenght*BinLenght) { Normalizer[1][bi] += re t ; if(FLAG){ hold = (struct b*)malloc(sizeof(struct b)); FLAG=0; point=hold; holdp = point; } else hold = (struct b*)malloc(sizeof(struct b)); i f ( h o l d == NULL) { exit(OK); > holdp->next = hold; hold->m.i = ( ( s o r t e d [ i b - l ] . i ) * (unsigned short)256); hold->m.i += (sorted [ i b - 1 ] . j ) ; hold->m.w = re t ; holdp = hold; num[pr][bi]++; } > if(!FLAG) hold->next = NULL; if(num[pr][bi]==0) return(NULL); return(point); Appendix B. Generation of the System Matrix for the Multiple Line Array 157 } f l o a t GetW(int i / * c o l * / , int j /*row*/, int pr, i n t bi) { f l o a t X.Y; f l o a t t e t a , tetaprim, co, s i , LineL, LineR, Be, Bs, dL, dR, same, valL, valR, temp; f l o a t dteta; X = coox(i); Y = cooy(j); dteta = ( f l o a t ) d i r e c t i o n * (-PI) / (float)PROJECTIONS; if(WHOLERot) t e t a = start.angle * PI / 180. + ( f l o a t ) ( p r ) * ( 2 . 0 ) * d t e t a ; else t e t a = start.angle * PI / 180. + ( f l o a t ) ( p r ) * d t e t a ; co = ( f l o a t ) c o s ( - t e t a ) ; s i = ( f l o a t ) s i n ( - t e t a ) ; tetaprim = t e t a - floor(teta/(PI/2.0)) * PI/2.0; /* / / i f ( t e t a p r i m < 0.0 ) tetaprim = PI/4.0 - tetaprim; */ i f (tetaprim < 0.0) p r i n t f ("Teta prim '/,f\n", tetaprim); i f ( t e t a p r i m > PI/4.0) tetaprim = PI/2.0 - tetaprim; Y = Y*co + X * s i ; LineL = BinPos(bi) - CamBinLenght/2.0; LineR = BinPos(bi) + CamBinLenght/2.0; Be = 0.7071068 * BinLenght * cos(tetaprim - PI/4.0); Bs = 0.7071068 * BinLenght * sin(PI/4.0 - tetaprim); co = (float)cos(tetaprim); s i = ( f l o a t ) s i n ( tetaprim); dL = LineL - Y; dR = LineR - Y; if(dR <= -1.0*Bc II dL >= Be ) return(O.O); if(dR >= Be && dL <= -1.0 * Be) return(BinLenght*BinLenght); if(dL*dR >= 0.0) { ^ same = TRUE; else { same = FALSE; dL = fabs(dL); dR = fabs(dR); if(dR < dL) { temp =dR; dR = dL; dL = temp; > i f (dR <= B c H i f ( d L < Bs){ valL = BinLenght * dL / co; else i f ( d L < Bc){ valL = (BinLenght / (2.0*co))*(2.0*Bc*dL - dL*dL - Bs*Bs)/(Bc-Bs); Appendix B. Generation of the System Matrix for the Multiple Line Array } else { valL = BinLenght*BinLenght/2.0; } if(dR < B s H valR = BinLenght * dR / co; > else if(dR < Be) { valR = (BinLenght / (2.0*co))*(2.0*Bc*dR - dR*dR - Bs*Bs)/(Bc-Bs); } else {. valR = BinLenght*BinLenght/2.0; } if(same){ return(valR-valL); else{ return(valR+valL); > else if(dR >Bc){ i f ( d L <= Bs){ valL = BinLenght * dL / co; } else i f ( d L < Bc){ valL = (BinLenght / (2.0*co))*(2.0*Bc*dL - dL*dL - Bs*Bs)/(Bc-Bs); > else printf("Wh00ps\n"); if(same) •C return(BinLenght*BinLenght/2.0-valL); } else •C return(BinLenght*BinLenght/2.0+valL); } } else { printf("NO POSSIBLE\n");} } int SaveMatrix(int i.ehar *name) s t a t i c FILE *wsk; int j , k,t; int number; unsigned short a,b; int l i c z n i k ; struct b *m; struct b *oldm; f l o a t area; s t a t i c FIRST=1; for(j=0; j<BINS*BINS; j++) memory [j/BINS] [j'/.B INS] = 0.; i f (FIRSTH if((wsk = fopen(name, "wb"))==NULL) p r i n t f ( " f i l e z l y " ) ; FIRST =0; for(j=0;j<CAMBIN;j++){ number = num[i][j]; m = mat[i] [j] ; fwrite(ftnumber,sizeof(int), 1,wsk); licznik=0; area = 0.0; if(number > 0){ do{ Appendix B. Generation of the System Matrix for the Multiple Line Array if(Normalizer[1][j] != 0.) m->m.w *= Normalizer[0][j]/Normalizer[1] fwrite(&(m->m), sizeof(matrix), 1, wsk); a = (m->m.i)/(unsigned short)256-(unsigned short)1; b = (m->m. i)'/,(unsigned short)256-(unsigned short) 1; area+=m->m.w; memory[a][b] = licznik + j ; oldm = m; m = m->next; free(oldm); licznik++; }while(m != NULL); printf ("licznik '/.d V,.2fW .licznik, area); i f ( l i c z n i k != number) { printf("Error occured in saving\nExiting\n"); exit(OK); > return(OK); } Appendix C Analytical Transmission Projection Generator (ATPG) For each bin in a projection, the probability of detecting a photon was calculated using the following formula: PJ = A • E / A c t • C A F • S C F • exp(-AttLength) • SolidAngle • dl (C.2) ; • J along line tines " where: Pj - probability of detecting a photon in the detector bin j , A - normalization constant (determined experimentally), Act - term proportional to the activity in the line source, J2iines faiongiine ' summation performed over all line sources and integration along each particular line source required because a photon can originate from any line and from any position in a line, CAF - collimator acceptance function (depends on the angle of the incoming photon), 160 Appendix C. Analytical Transmission Projection Generator (ATPG) 161 SCF - source collimation function (describes the collimation of the source), AttLength - attenuation term equal to the integral of the attenuation coefficient along the path from the position on the line source / to the bin, and SolidAngle - term proportional to the solid angle created by bin "seen" from the line source. C. l Collimator Acceptance Function, CAF In all the simulations, the collimator on the detector was modelled analytically. The hexagonal shape of each hole in the collimator was approximated by a circle with the same area as the hexagonal shape. The C A F was calculated as the fraction of the area of the hole at the bottom of the collimator "seen" by the incoming photon. This area depends on the incident angle 6 of the photon. Fig. (C.7) Upper hole photon Lower hole y Figure C.7: Definition of the collimator acceptance function (CAF). Shaded is the area "seen" by photons incoming with angle 6 The C A F can be expressed by a phenomenological equation: Appendix C. Analytical Transmission Projection Generator (ATPG) 162 FCAF{9) — P(upper) • P(lower\upper) (C.3) where P(upper) is the probability that photon goes through the upper hole, and P(lower\upper) is the probability that the incoming photon goes through the lower hole under the assump-tion that it has gone through the upper one. P(upper) is the fraction (fh) of the area of the holes to the total area of the collimator. Usually that ratio is about 80%. The value of P(lower\upper) was calculated analytically with assumption of small angle 0. The total C A F is: FCAF(0) = { fh • (TT - sin(2arcsin(f)) - 2 arcsin(f ) ) / T T , if|0| < 0 0, if|0| > 0 (C.4) Here 0 was the acceptance angle of the collimator. Using Eq. C.4 the C A F can be plotted for collimators commonly used in nuclear medicine. For low energy high resolution (LEHR) collimator and low energy all purposes (LEAP) collimator manufactured by Sophy, the parameters of the collimators are shown in Table C . l . Table C . l : Parameters of the collimators. Data for the Sophy DST Camera. Collimator Hole Septum Septum Acceptance h type Diameter [cm] Thickness [cm] Length [cm] Angle [°] L E A P 0.19 0.02 3.2 3.6 0.83 L E H R 0.19 0.02 3.8 3.0 0.83 The collimator acceptance function is plotted in Fig. C.8. Appendix C. Analytical Transmission Projection Generator (ATPG) 163 Figure C.8: The collimator acceptance function (probability that photon will be detected) vs. the incidence angle 0. C.2 Source Collimation Function, SCF The SCF function is different for the M L A and for the CLS. C L S The CLS was collimated by two lead strips parallel to the line source. Such a design of collimator created a penumbra. The SCF for the CLS was the ratio of the cross section area of the line source visible from a given angle (Fig. C.9), to the total cross section area of the line. M L A The source collimation of M L A was done by placing a sheet L E H R collimator on top of the line sources, so the SCF was effectively equal to the C A F . C .3 Attenuation Length, AttLength The attenuating medium was defined as a number of ellipsoidal objects. Each object Appendix C. Analytical Transmission Projection Generator (ATPG) 164 I I Source Collimator SCF I Corresponds to / shaded area Figure C.9: The SCF for the CLS. was described by 7 parameters: 3 coordinates for the centre of any given ellipsoid, 3 radii, and the attenuation coefficient. A cylinder could be defined as an ellipsoid with radius along the axis of rotation equal to oo. Objects might contain each other but they cannot overlap. A T P G finds AttLength by calculating the equation of the line which, in 3D, connects a bin with some point in the line source, and then analytically calculating the intersection length between this line and the ellipsoid (or ellipsoids if more than one object is used). These distances are weighted by the attenuation coefficient then added together (taking into account their relative position) to give the value of AttLength. The factors inside the summation and integration in Eq. C.2 gave a count distribution on the detector, but these counts were not normalized properly to the activity in the source because the efficiency of the detector was not taken into account and Act was just C . 4 N o r m a l i z a t i o n c o n s t a n t , A . P o i s s o n N o i s e Appendix C. Analytical Transmission Projection Generator (ATPG) 165 a term proportional to the activity in a source. To find a proper normalization we used the easiest and the most accurate way. We calculated the shape of the beam for one line source using Eq. C.2 and than we fitted this shape to that experimentally obtained image from single line source with known activity; this yielded to a proper value of normalization constant A. Having the value of the counts in a given bin, Poisson noise was added to the data using the algorithm described in [77]. Appendix D Verification of Eq . 4 .10 A simple simulation was performed in order to verify Eq. 4.10. This equation was used for the derivation of the system matrix for CLS and M L A in Chapter 4. We here verify this equation for CLS. The equation has the form: Ij ~ Bj • exp(-A j ) (D.5) where Ij and Bj were the number of counts detected in the jth detector bin for the transmission and the blank scan, respectively. The Aj defined in Sec. 4.1.3 was an average of A(x) (= \og(B(x)/I(x))) over the projection bin which extended from position xi on the detector to position x2. It was expressed by Aj = ^ — ^ (D.6 x2 - Xi ' B(x) and I(x) are the photon flux densities at the detector position x for the blank and the transmission scans, respectively Transforming Eq. D.5 we have A ^ l n ( ^ ) (D.7) 166 Appendix D. Verification of Eq. 4.10 167 The Aj in Eq. D . 7 is the approximation of Aj expressed by Eq. D . 6 . In this appendix the noiseless projection data of a phantom were generated analytically in 2D. The created projections were then binned creating values of Aj with 3 levels of approximation of the exact value Aj. They are: E x a c t In a real scan, the exact contiguous transmission projection data will never be available because we acquire data in bins, but in this simulation we created con-tinuous exact projection data (B(x),I(x)), and then calculated projection elements by using the definition of Aj from Eq D . 6 : Aj = J * ' l n ^ ( D . 8 ) x 2 - x i T h e B e s t A v a i l a b l e In this level of approximation we found Aj for all the bins using Eq. D . 7 and then we rebinned the value of Aj into new bins corresponding to CLS geometry (see Sec. 4.1.1). This approximation is called the best available because given data acquired in bins, it is the best approximation of the exact continuous distribution of counts. The calculation of Aj was done using the formula A . = ZkSkj-HBk/h) ( D 9 ) J2k sk,j where Ik and Bk are the average values of counts over the bin k in the transmission and blank scans (Fig. D .10), and where Sk,j is a fraction of length of the bin k inside the new bin j . U s e d The method used in this thesis which uses the approximation of Eq. D.5. Although worse than the other, the method substantially reduces noise by adding counts from several bins. Here Aj was again calculated from Eq. D . 7 but the values of counts in this equation, Ij and Bj, were the values of the counts in the new bins. A h ( D . 1 0 ) Appendix D. Verification of Eq. 4.10 168 Schematically the approximations are shown in Fig. D.10. ( c) Detector bins I 1 1 1 1 1 1— New bins Figure D.10: The approximations made during the binning of ID projection data, (a) Exact projection uses continuous values acquired in transmission and blank scans (b) Data is acquired in bins; the value of the bin is calculated as an average over the bin from the continuous distribution, (c) further rebinning of (b) into new bins. Method (c) was used in this thesis for CLS system. The M L A data was calculated as described in (b). Using the data obtained by these three levels of approximation the image of analytical phantom was reconstructed and the accuracies of reconstructions from different approx-imation levels were compared. Fig. D . l l shows the result of the least squares reconstruction of the 2D phantom. The reconstruction method was previously described (Sec. 4.1.3). E R R O R is defined by Appendix D. Verification of Eq. 4.10 169 Eq. 5.5 on page 60 and is a Cartesian type distance of the reconstructed image from the true phantom. In Fig. D . l l the Exact calculation of the projection elements leads to the reconstructed image converging1 to the true values image when the number of iterations increases. The ERROR for the other methods behaves similarly; increasing with number of iterations as expected, but the rate of the increase is approximately the same, which means that Used does not introduce many more inconsistencies than are already in Best Available. Although the reconstructed image from Best Available is closer to the true image because it uses finer binning and it is a better approximation of Eq. D.5, Used behaves almost equally well. -0.50 -0.55 ^ -0.60 O C O C C -0.65 C C LU D) -0.70 O -0.75 -0.80 i • r Exact The Best Available Used _i L 20 40 60 Iterations 80 100 Figure D . l l : Error in least squares reconstruction of the simulated projection data binned with different methods. 1The value of ERROR will never reach zero because the simulated phantom was mathematical whereas the calculation of ERROR is made on the digital (pixelized) data. Appendix E Implementation of Attenuated Pr o j ector / B ackpr o j ector This appendix describes the implementation of the fast attenuated projection operation which was used in chapter 7. It is based on the parallel system matrix (Appendix B) also used for reconstruction of the attenuation maps for M L A . In S P E C T , the element of the system matrix, corresponds to the probability that a photon emitted in pixel i is detected in bin j. A system matrix with elements proportional to this probability was already calculated (Appendix B). With the attenuating medium, this probability has to be modified to allow for attenuation between pixel i and bin j. Using the attenuation maps and matrix u>;j already created, we could very easily modify the matrix taking attenuation into account because of the special ordering of the elements in this matrix (Appendix B). The modification goes according the following equation: i-i 1 = Whj - eXP(X] Wk,j • Vk + 7y ' wi,j ' Vi) ( E . l l ) k=l Z where u>,-j was the attenuation modified system matrix and the values of /i's came from 170 Appendix E. Implementation of Attenuated Projector/Backprojector 171 the attenuation map. Now by using tujj instead of Wij in the projector/backprojector equation (Eq. 4.3) we get an attenuated projector/backprojector. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items