University of Southern California Dissertations and Theses

Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

(USC Thesis Other)

Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content FAST ITERATIVE IMAGE RECONSTRUCTION FOR 3D PET AND ITS EXTENSION TO TIME-OF-FLIGHT PET by Sanghee Cho A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) December 2008 Copyright 2008 Sanghee Cho Dedication To my mother who must be watching over us in heaven. ii Acknowledgments My first thanks go to my Ph.D advisor, Richard M. Leahy. His inspiration, thoughtful guidance and support made this thesis possible. His positive attitude and kind supervi- sion always kept my spirits up, and there was always his insightful advice whenever I had problems on both academic and personal levels. Tremendous help and disscussion were provided by Quanzheng Li and Sangtae Ahn. I cannot thank them enough. They were mentors who had guided me since I joined the lab, and also great colleagues who gave me a lot of help with insightful ideas and valuable suggestions. I also would like to thank the dissertation committee members, Prof. C.-C Jay Kuo, Prof. Manbir Singh, Prof. Arion Chatziioannou, for their valuable suggestions. I would like to thank Anne M. Smith, Bing Bai and Mu Chen at Siemens Preclin- ical Solutions. It was a great honor for me to have my work used in their commercial software. While working with them, I learned many important things for PET systems in real world. I also acknowledge their financial support. Next, I acknowlege the con- tributions made by Erkan Mumcuoglu, Jinyi Qi, Bing Bai, Quanzheng Li and Evren Asmain in developing our PET image reconstruction code. I was a real beneficiary of the accrued contributions. I also would like to thank my country, Korea, for providing 2 years scholarship. I feel very fortunate to finish up my Ph.D work in such a friendly, and undoubtedly excellent research lab. I specially thank Anand Joshi and Sangeetha Somayajula for iii sharing an office from the EEB basement to the current office in a new building, and also thank other current officemates, Joyita Dutta and Zheng Li. I appreciate their kind and generous support. I also would like to thank all other former and current lab members, Dimitrios Pantazis, Belma Dogdas, Abhijit Chaudhari, Evren Asma, Esen Kucukaltun- Yildirim, Felix Darvas, Hua Hui, Juan L. P. Soto, Yu-Teng Chang and David Wheland, Syed Ashurafulla, Shailesh Chiplunkar, Ran Ren and Yanguang Lin. Personally, it was very tough to pursue the Ph.D program while supporting a family. However, they were my real energy source, and their emotional support was the most powerful thing for me to stay in the right track during those years. I truly thank my wife, Yong-Chun, and three adorable kids, Priscilla (Young-seo), Caleb (Kyoung-chan) and Stephanie (Kyoung-eun). I especially cannot thank my wife enough for her hard work in taking care of all the kids and family matters. I also deeply thank my parents, uncles and aunts, brothers and sisters, in-laws, and all other relatives for their support and encouragement. iv Table of Contents Dedication ii Acknowledgments iii Abstract xiv Chapter 1 Introduction 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Contribution and Dissertation Organization . . . . . . . . . . . . . . . 3 Chapter 2 Background 6 2.1 Physical Principle of PET . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 3D PET Data Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Iterative PET Image Reconstruction . . . . . . . . . . . . . . . . . . . 11 2.4 Fourier Rebinning for Fast 3D PET Image Reconstruction . . . . . . . . 15 2.5 Time-of-Flight (TOF) PET . . . . . . . . . . . . . . . . . . . . . . . . 18 Chapter 3 Fast Iterative 3D PET Image Reconstruction 25 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.1 Factored System Model . . . . . . . . . . . . . . . . . . . . . . 27 3.2.2 Two Complicating Factors Encountering When Using Fourier Rebinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Fast 3D Projector Using Exact Inverse Rebinning . . . . . . . . . . . . 29 3.3.1 2D Projection: 3D Image to 2D Direct Sinograms . . . . . . . . 31 3.3.2 Inverse Rebinning: 2D Direct Sinograms to Pseudo 3D Sinogram 33 3.3.3 Post Corrections: Pseudo 3D Sinogram to 3D Sinogram . . . . . 35 3.3.4 IRB Projector Pair: From 3D Image to 3D Sinogram and Back . 37 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4.1 Computation Cost and Accuracy . . . . . . . . . . . . . . . . . 39 3.4.2 Application to MAP Reconstruction: Resolution Properties . . . 42 3.4.3 Application to MAP Reconstruction: Phantom and in vivo Studies 49 v 3.4.4 Application to MAP Reconstruction: Reconstruction Time Com- parison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Chapter 4 Unified Framework for Fourier Rebinning 55 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Analytical Data Model and Generalized Projection Slice Theorem . . . 56 4.2.1 Line Integral Model with Spatially Invariant TOF Kernel Function 57 4.2.2 Taking Fourier Transform in TOF Bin Direction . . . . . . . . . 58 4.2.3 Analytical Properties of 2D TOF PET Data . . . . . . . . . . . 60 4.2.4 Generalized Projection Slice Theorem . . . . . . . . . . . . . . 68 4.3 All Exact Mappings Between Different PET Data Sets . . . . . . . . . . 70 4.3.1 Mapping D: 3D non TOF and 2D non TOF PET Data . . . . . . 70 4.3.2 Mapping B: 3D TOF and 2D TOF PET Data . . . . . . . . . . . 73 4.3.3 Mapping C: 3D TOF and 2D non TOF PET Data . . . . . . . . 76 4.3.4 Mapping A: 3D TOF and 3D non TOF PET Data . . . . . . . . 77 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Chapter 5 Fourier Rebinning of 3D TOF PET Data for Fast Reconstruction 81 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.2 Rebinning of 3D Time of Flight (TOF) Data to 2D TOF Data . . . . . . 83 5.2.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.2.2 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.3 Rebinning of 3D TOF Data to 2D and 3D non TOF Data . . . . . . . . 90 5.3.1 Non TOF Fourier Rebinning and Approximation . . . . . . . . 91 5.3.2 TOF Rebinning Mappings and Approximations . . . . . . . . . 93 5.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Chapter 6 Summary and Future Work 109 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 References 112 vi List of Tables 2.1 Sinogram dimension and 3D PET data size for four Siemens PET scan- ners. Note that * indicates a sinogram bin size of 4 byte integers, and ** is a sinogram bin size of 2 byte integers. . . . . . . . . . . . . . . . . . 12 2.2 Property comparison of various PET scintillators. . . . . . . . . . . . . 19 3.1 MicroPET Focus 220 Scanner Specifications . . . . . . . . . . . . . . . 40 3.2 Comparison of computation times for forward and back projections with different projectors. For the IRB projectors, the run times for inverse rebinning and 2D projection are also shown separately. See text for details. 43 3.3 Comparison of the MAP reconstruction time (Unit: Minutes). . . . . . 53 4.1 Summary of all mapping equations derived from the generalized projec- tion slice theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.1 Simulated 3D TOF cylindrical PET system parameters . . . . . . . . . 85 5.2 Exact and approximate points in(ω x ,ω y ) plane for the inverse rebinning for mappings C and A. . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.3 Simulated 3D TOF cylindrical PET system parameters . . . . . . . . . 100 vii List of Figures 2.1 Physical principle of PET imaging. An unstable neutron deficit isotope undergoes a proton decay producing a neutron, a positron and a neu- trino. The two products of the decay, positron and neutrino, are emitted from the isotope, and the emitted positron is annihilated with a nearby electron. Due to the positron-electron annihilation, two 511 keV pho- tons are produced, and emitted in opposite directions of approximately 180 ◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 The emitted photon pair from the positron-electron annihilation in Fig- ure 2.1 is detected by PET detector pair A and B. The photon pair is called a coincident event when their detected time difference is smaller than the preset coincidence timing window,T . . . . . . . . . . . . . . 7 2.3 A line ( or tube) of response (LOR) in 3D PET acquisition mode between two detectors in two different rings. Each histogrammed projection data bin contains the number of detected coincident events within the line of response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 The detectors are paired in the way to give parallel LORs for a fixed angle to form a parallel projection data. . . . . . . . . . . . . . . . . . 9 2.5 (a) The 3D PET sinogram is the combination of 2D direct sinograms and all oblique sinograms. Notice that the number of oblique sinograms for a fixed oblique angle is less than the number of 2D direct sinograms due to the finite axial aperture of the scanner. (b) Axial mashing is performed to reduce the 3D PET data size. The parameter ‘span’ determines the number LORs which will be combined. . . . . . . . . . . . . . . . . . 10 2.6 Example of the Michelogram which is useful for illustrating 3D PET data reduction by span and maximum ring difference (MRD). The Mich- elogram is shown for a 16-ring scanner, where the span is 5 (2 + 3), and MRD is 12 (13-1) so that the sinograms marked as circled asterisks are discarded. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 viii 2.7 Fast 3D PET data reconstruction using Fourier rebinning. 3D PET data is rebinned to 2D stacked direct sinograms so only a slice-by-slice 2D reconstruction is required after Fourier rebinning. Fourier rebinning is an algorithm to reduce the dimensionality of 3D PET data. . . . . . . . 14 2.8 (a) Top and (b) side view of a cylindrical 3D PET scanner. . . . . . . . 16 2.9 Time-of-flight (TOF) information can be utilized for better image recon- struction. In (a), when backprojected, the uncertainty in the measured data is distributed along the whole LOR, but in (b), uncertainty is reduced by weighted backprojection with the TOF kernel function. . . . . . . . 21 2.10 TOF PET data can be modelled using a line integral weighted by a TOF kernel function which is a probability density function of the measure- ment error associated with the time difference of two detected photons at the detector pair. It is often assumed spatially invariant, and approxi- mated as a Gaussian distribution with standard deviation determined by the timing resolution of the system. . . . . . . . . . . . . . . . . . . . . 22 2.11 TOF SNR gain is approximately expressed by the FWHM of the TOF kernel function and object size. . . . . . . . . . . . . . . . . . . . . . 23 3.1 LOR sample interval variation as a function of sinogram radial index for the microPET Focus 220 scanner . . . . . . . . . . . . . . . . . . . . . 29 3.2 The variation ofδ (the tangent of oblique angleθ) as a function of radial index for oblique sinograms for the microPET Focus 220 scanner. Each curve shows the variation in δ for a fixed ring difference, with a maxi- mum ring difference (MRD) of 47 and a span of 3 . . . . . . . . . . . . 30 3.3 The 3D geometric projector P geom can be approximated using inverse rebinning projectors via the following three steps: (i) 2D projection, (ii) inverse rebinning (IRB), and (iii) post corrections. Three different vari- ations based on IRB are investigated here: Case A uses a 2D geometric projector with corrections for nonuniform sampling in s and δ, Case B uses the same 2D geometric projector but without the corrections, and Case C is the same as Case A but uses a 2D line integral projector in place of the 2D geometric projector. . . . . . . . . . . . . . . . . . . . 31 3.4 Block diagram of inverse rebinning operatorB δ IRB which computes the oblique sinograms at fixed oblique angle index δ from 2D direct sino- grams. The dotted line box indicates calculation for eachφ whereφ m is the largest index forφ. . . . . . . . . . . . . . . . . . . . . . . . . . . 34 ix 3.5 Illustration of oblique sinogram resampling procedure: (a) uniform radial and axial angular sampling pattern ins andδ after applying IRB to com- pute oblique sinograms from the 2D direct sinograms; (b) radial and axial angular sampling pattern after IAC; (c) radial and axial angular sampling pattern after ARS and IAC. . . . . . . . . . . . . . . . . . . . 36 3.6 Sinograms for a ring difference of 38 computed from the digital 3D Hoffman brain phantom using two different transaxial voxel sizes: (left) 256x256x95 .4mm voxel size, (right) 256x256x95 0.8mm voxel size. (a) The reference sinograms from the 3D geometric projector; (b) differ- ence images between reference and IRB projector Case A in figure 3.3; (c) difference images for Case B; (d) difference images for Case C. The normalized RMS error in (3.7) is indicated below each difference image. 41 3.7 Radial resolution (FWHM) in three different planes for the 3D geomet- ric and inverse rebinning projectors for 100 iterations of MAP and β = 0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.8 Tangential resolution (FWHM) at three different planes for the 3D geo- metric and inverse rebinning projectors for 100 iterations of MAP andβ = 0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.9 Axial resolution (FWHM) for three lines parallel to the central axis of the scanner for the 3D geometric and inverse rebinning projectors for 100 iterations of MAP andβ = 0.001. . . . . . . . . . . . . . . . . . . 48 3.10 Reconstructed transaxial sections through the brain phantom (c) using the 3D geometric (a) and inverse rebinning (b) projectors; (d) is the dif- ference between (a) and (b); (e) shows a profile along the white dashed line in (c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Sagittal (left), transaxial (middle) and coronal (right) views of recon- structed images on the microPET Focus 220 of an FDG study of a macaque monkey using (a) MAP with the 3D geometric projector; (b) MAP with IRB projector, Case A; (c) MAP with IRB projector, Case B; (d) FORE + 2D OSEM; (e) 3DRP with a ramp filter. . . . . . . . . . . . 52 4.1 All possible mappings between 3D TOF PET data and other lower dimen- sional data sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2 (a) Top and (b) side view of a cylindrical 3D PET scanner. For each line of response (LOR), the object is multiplied by the TOF kernelh and integrated along the line to form the TOF data. . . . . . . . . . . . . . 58 x 4.3 (a) The variabless andt constitute an orthogonal coordinate rotated by φ fromx-y coordinate. (b) Therefore, by taking Fourier transform with respect to s and t, their frequency coordinate is also rotated by φ. The strip weighted by H in (b) shows the extra information produced by having TOF information. When TOF information is not available, we have information only in the line whereω t = 0. . . . . . . . . . . . . . 61 4.4 The frequency response of the conventional 2D PET sinogram (Radon transform) has bow-tie shaped frequency response where the bow-tie boundary is determined by the radius of non zero region of object func- tionf. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.5 Example of the magnitude response of the series in (4.11). The figures shows the modified bow-tie region bounded by hyperbolic lines where the outside of the bow-tie region is nearly zero. . . . . . . . . . . . . . 64 4.6 2D TOF sinograms and their magnitudes after taking 1D Fourier trans- form int; The horizontal axis is the sinogram angleφ∈ [0,2π] and the vertical axis is the radial variables. . . . . . . . . . . . . . . . . . . . . 65 4.7 2D Fourier transformP(ω s ,ω φ ;ω t ) ofP(s,φ;ω t ) in (4.11) for four dif- ferent ω t frequency indices. The horizontal axis represents ω φ and the vertical axisω s . Also shown in (d) is the hyperbolic curve bounding the effective support of the bow-tie for TOF data in the 5th frequency bin. . 66 4.8 Bow tile filtering to noisy 2D TOF sinograms in (4.5); (a) magnitude response of 2D TOF sinogram after FT to TOF bin direction at 3-th TOF frequency index from zero frequency, (b) magnitude plot of 2D FT of (a), (c) Bow-tie filter in 2D frequency domain, (d) Radial filter in 2D frequency domain, (e) magnitude plot of Bow-tie filtered sinogram of (a), (f) magnitude plot of radial filtered sinogram of (a) . . . . . . . . . 67 4.9 Graphical interpretation of the exact inverse Fourier rebinning mappings from 2D non TOF data to 3D non TOF data. Inverse Fourier rebinning is equivalent to finding the coordinates (ω ′ s ,φ ′ ) of a point (for example, pointB) on a trajectory (for example, the line joining the pointsA and B) mapped from the 3D sinogram. In the approximate inverse rebinning, |CB| and∠ABC are approximated by|AB| and−δω z /ω s , respectively. The inside region of the dotted line circle is not covered by 3D non TOF PET data℘(ω s ,φ,ω z ,δ;0). So there is no rebinning mapping available for the 2D non TOF PET data℘(ω ′ s ,φ ′ ,ω z ,0;0) in this area from the 3D non TOF PET data℘(ω s ,φ,ω z ,δ;0). . . . . . . . . . . . . . . . . . . . 72 xi 4.10 Graphical interpretation of the exact inverse Fourier rebinning mappings from 2D TOF to 3D TOF PET data. The exact inverse Fourier rebinning is equivalent to finding the coordinate information of the line segment CB (|CB| and∠ABC) for the line segment AB. . . . . . . . . . . . . . 75 5.1 Comparison of rebinned TOF sinograms with corresponding directly computed TOF sinograms (TOF timing resolution 500 psec); Columns show TOF bin indices from 4 to 7; rows show the calculated direct sino- gram (Reference), the rebinned sinogram (TOF-FOREX), and their dif- ference (Difference). . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2 Comparison of rebinned TOF sinograms with corresponding directly computed TOF sinograms (TOF timing resolution 500 psec) shown in Figure 5.2; (a) radial profile through the sinograms for 4th TOF bin at 40-th angle; (b) radial profile through the sinograms for 7th TOF bin at 40-th angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3 Comparison of (top) noisy 2D direct sinogram with (bottom) noise reduc- tion resulting from rebinning of oblique data into direct planes for (a) 500 psec, and (b) 250 psec resolution. . . . . . . . . . . . . . . . . . . 88 5.4 3D reconstructions from 200 million counts for the central axial slice: (a) true image; (b) 3D non-TOF reconstruction: FOREX followed by direct 2D Fourier reconstruction; (c) 3D TOF reconstruction using TOF- FOREX followed by 2D TOF direct Fourier reconstruction with TOF timing resolutions of 500 psec: (d) as for (c) but with 250 psec resolution. 89 5.5 Graphical interpretation of the exact and approximate inverse Fourier rebinning mappings from 2D non TOF data to 3D non TOF data. Inverse Fourier rebinning is equivalent to finding the coordinates (ω ′ s ,φ ′ ) of a point (for example, pointB) on a trajectory (for example, the line join- ing the pointsA andB) mapped from the 3D sinogram. In the approx- imate inverse rebinning, |CB| and∠ABC are approximated by |AB| and−δω z /ω s , respectively. . . . . . . . . . . . . . . . . . . . . . . . . 92 5.6 Graphical interpretation of the exact and approximate inverse Fourier rebinning mappings from 3D non TOF data to 3D TOF data. The exact inverse Fourier rebinning is equivalent to finding the coordinate infor- mation of the line segment CB (|CB| and∠ABC) for the line segment AB. By the approximation in (5.9),|CB| and∠ABC are approximated by|DB| and∠ABD, respectively. . . . . . . . . . . . . . . . . . . . . 95 xii 5.7 Averaged mapping errors of approximate rebinnings for different ring differences and for a fixedφ = 10. The error value is normalized by the image frequency sample interval (Δω x ) in the transverse plane where Δω x = 1/2/256 mm −1 in our simulation. (See Table 5.3.3 for parame- ters). Notice the error is very large whenω s is small, so special consider- ation will be required in this region, similarly to the FORE implementation. 97 5.8 Averaged mapping errors of approximate rebinnings for different ring differences and for a fixedφ = 10. The error value is normalized by the image frequency sample interval (Δω x ) in the transverse plane where Δω x = 1/2/256 mm −1 in our simulation. (See Table 5.3.3 for parame- ters). Notice the error is very large whenω s is small, so special consider- ation will be required in this region, similarly to the FORE implementation. 98 5.9 A comparison of the mean and variance of rebinned sinograms, obtained by FORET, and non TOF sinograms for two ring differences: (a) Mean of rebinned sinograms for RD (ring difference) = 0, (b) Variance profiles of rebinned and non TOF sinograms at 160-th angle (RD=0). . . . . . . 102 5.10 A comparison of the mean and variance of rebinned sinograms, obtained by FORET, and non TOF sinograms for two ring differences: (a) Mean profiles of rebinned and non TOF sinograms at 160-th angle (RD=54), (b) Variance profiles of rebinned and non TOF sinograms at 160-th angle (RD=54). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.11 Resolution (FWHM) versus pixel variance plot for four different loca- tions. (a) pixel location (128,128); (b) pixel location (128,160) . . . . . 105 5.12 Resolution (FWHM) versus pixel variance plot for four different loca- tions. (a) pixel location (128,220); (b) pixel location (110,110) . . . . . 106 5.13 A comparison of 3D TOF data reconstruction by ‘FORET+MAP’ (top row) and ‘non TOF+MAP’ (bottom row) in the transverse view (first col- umn), coronal view (second column) and sagittal view (third column). . 107 xiii Abstract Positron emission tomography (PET) is a functional biomedical imaging technique that provides in vivo information about physiological processes within the body by recon- structing the 3D distribution of a radiolabelled tracer. It is challenging to estimate the tracer distribution accurately since photon-limited PET data have high statistical vari- ance. Statistical iterative reconstruction methods have shown superior image quality and better quantitative results compared to analytical reconstruction methods. However, the iterative methods involve huge computational loads, particularly for fully 3D PET, consequently limiting their routine use in practice. Clinical 3D time-of-flight (TOF) PET scanners, which become available due to recent developments of fast scintillators, involve an even higher computational cost for image reconstruction because of extra TOF information. A successful pragmatic approach that reduces computation cost for fully 3D PET is to use Fourier rebinning (FORE), reducing 3D PET data to 2D data, combined with a 2D reconstruction method such as the ordered subsets expectation- maximization (OSEM) algorithm. The combination of FORE and 2D OSEM is now routinely used in both clinical and small animal imaging systems. However, data dimen- sionality reduction using FORE can result in a loss in resolution or noise performance. In the first part of this dissertation, we develop fast projectors that map a 3D image into 3D PET data by using the exact inverse (Fourier) rebinning operator. The inverse rebin- ning operator is used for calculating 3D PET data from 2D PET data, and combined xiv with a 2D projector to construct a fast 3D projection operator. The inverse rebinning operator is implemented cost-effectively using the fast Fourier transform (FFT), so that the computation time is reduced by an order of magnitude. We utilize this fast projector in the context of fully 3D PET maximum a posteriori (MAP) reconstruction methods. We reduced the reconstruction time substantially while retaining the resolution recovery properties. In the second part of this dissertation, we develop fast 3D image recon- struction methods for TOF PET. We derive a unified framework based on a generalized projection slice theorem. In this framework, we derive mapping equations between dif- ferent data sets such as 3D TOF, 2D TOF, 3D non TOF and 2D non TOF PET data. Using the mapping equations, we develop a TOF Fourier rebinning method which rebins 3D TOF PET data to 2D TOF data. We also develop other exact and approximate rebinning methods where 3D TOF PET data are rebinned to 2D or 3D non TOF data. Our sim- ulation studies showed that the reconstructed images after TOF Fourier rebinning give superior resolution-variance properties compared to the non TOF case where the TOF information is ignored. xv Chapter 1 Introduction 1.1 Overview Positron emission tomography (PET) [100] is a functional imaging modality that pro- vides in vivo information about targeted physiological processes. In PET, a radiotracer is injected into a patient or an animal, which is a radiolabelled molecule selected or designed to provide a spatial or spatial-temporal map of a physiological process in the subject. For example, 18 FDG, an analog of glucose, is used to provide a measure of glucose metabolism [103]. Other tracers measure neuroreceptor density [46, 134], regional blood flow or volume [48, 121], rates of gene expression [12], and cardiac per- fusion [114]. It is challenging to estimate the tracer distribution accurately because PET data usu- ally have high statistical variance due to the limited dosage and scan time. Statistical iterative reconstruction methods [49, 76, 107, 126] have shown superior image quality and better quantitative results [26,30,44,87,105] compared to analytical reconstruction methods [23, 25, 71, 95, 119]. However, the statistical iterative methods involve huge computational loads, particularly for fully 3D PET [20], consequently limiting their routine use in practice. This high computational cost mainly comes from projection operations which map a 3D image to the large PET data [6]. The computational issues become even more severe with the recently available 3D time-of-flight (TOF) PET sys- tems [19, 118] because of the increased data size due to the extra TOF information. 1 The first part of this thesis is focused on reducing the high computational cost involved in fully 3D PET iterative image reconstruction. A successful approach in 3D PET image reconstruction is to use Fourier rebinning (FORE) [29, 32], reducing 3D PET data into 2D stacked sinograms, combined with a 2D iterative image reconstruc- tion method such as the ordered subsets expectation-maximization (OSEM) [22]. The combination, FORE+2D OSEM, is now in routine use in both clinical [75] and small animal imaging systems [74]. While this pragmatic solution results in relatively short reconstruction times, the data dimensionality reduction that occurs during the initial rebinning step induces a loss in resolution or noise performance compared to fully 3D iterative reconstruction methods. Here we retain the fully 3D iterative reconstruction framework, but develop a fast projector using inverse Fourier rebinning [32] to reduce the computational cost. Inverse Fourier rebinning is used to calculate 3D PET data from stacked 2D sinograms, so we combine the inverse rebinning operator with a 2D projector, which maps the 3D image to stacked 2D sinograms, to to compute a fully 3D forward projection. The inverse rebinning operator is implemented cost-effectively using the FFT, and the computational cost is reduced by approximately an order of mag- nitude. We use this fast projector in the context of a fully 3D maximum a posteriori (MAP) PET reconstruction method [53,87,108]. Substantial reconstruction time reduc- tions were achieved while retaining the resolution recovery properties of the original MAP method. The second main contribution in this thesis is to develop a method which can reduce the computation cost in 3D time-of-flight (TOF) PET data reconstruction. It has long been known that the TOF information can be used to improve image quality in PET. However, due to the lack of high performance detectors, TOF PET had not been practical until the recently developed fast detectors such as LSO [77, 85, 90] and 2 LaBr 3 [66, 117] became available. The TOF timing resolution in current TOF PET sys- tems is now approximately 500 ps corresponding to a spatial resolution along each line of response (LOR) of 70 mm at full-width-at-half-maximum (FWHM) [19, 117]. It has been reported that this timing resolution can be beneficial in clinical studies where the size of the subject is relatively large [52, 89, 92, 124, 127]. However, fully 3D TOF PET iterative reconstruction is prohibitively expensive due to the increased data size results from the extra TOF information. We were therefore motivated develop a rebinning algo- rithm, similar to FORE in 3D PET, which can effectively reduce the 3D TOF PET data to a lower dimensional representation while still retains the benefit of the TOF infor- mation. First, we develop a unified framework based on a generalized projection slice theorem. In this framework, we derive various data mappings between 3D TOF PET data and its lower dimensional versions such as 2D TOF PET, 2D non TOF PET and 3D non TOF PET data. Using the mapping equations, we develop an exact rebinning method to rebin 3D TOF PET data to 2D TOF PET data. Then, we develop other novel rebinning methods where the 3D TOF PET data are rebinned to 2D or 3D non TOF PET data. Through simulation studies, the rebinng methods were evaluated, and showed to have the better performance than the non TOF PET case where the TOF information is ignored. Although we use the Fourier rebinning concept here for data reduction pur- poses, developing a fast projector for fully 3D TOF PET iterative reconstruction will be a natural extension of this approach since we expect fully 3D TOF PET iterative recon- struction will produce superior image quality compared to rebinning-based TOF PET data reconstructions. 1.2 Contribution and Dissertation Organization The original contributions of this thesis can be summarized as follows. 3 • We designed a fast projector based on Fourier rebinning. When incorporated into a fully 3D MAP reconstruction, the reconstruction time was significantly reduced without significant resolution loss [24]. • We developed a unified framework for deriving mappings between 3D TOF, 3D non TOF, 2D TOF and 2D non TOF PET data sets based on a generalized projec- tion slice theorem. The framework enables us to interpret the rebinning equations geometrically including the original Fourier rebinning equation, and can be used as a tool for implementing data reduction (rebinning) algorithms, or designing fast projectors using inverse rebinnings [17, 18]. • We developed TOF Fourier rebinning algorithms for fast 3D TOF PET data recon- struction based on the unified framework. We first developed an exact 3D TOF Fourier rebinning algorithm where 3D TOF PET data are rebinned to 2D TOF PET data. Next, we developed other rebinning algorithms where the 3D TOF PET data are rebinned to 2D and 3D non TOF PET data [17, 18]. This dissertation is organized as follows. In Chapter 2, as background, we first review the basic physical principle and data structure of PET. The iterative PET image reconstruction methods and Fourier rebinning concepts are reviewed next, followed by a review of time-of-flight (TOF) PET. The first contribution developing a fast projector using an exact inverse rebinning operator, is described in Chapter 3. The fast projector is incorporated into fully 3D MAP reconstruction code to replace the original matrix based geometric projector for speed up. The reconstruction performance using the projector was evaluated using simulations and real data. In Chapter 4, we describe the unified framework that derives the mappings between different PET data types (3D TOF, 2D TOF, 3D non TOF and 2D non TOF data), as shown in Figure 4.1. All the mappings were derived based on a generalized projection slice theorem derived from an analytical data 4 model for 3D TOF PET data. In Chapter 5, fast 3D TOF PET reconstruction methods based on TOF PET data rebinnings are described. Finally, in Chapter 6, we provide a concluding summary, and future work is discussed. 5 Chapter 2 Background 2.1 Physical Principle of PET Positron emission tomography (PET) requires that a radiotracer be injected into the sub- ject. The radiotracer is a radioactive molecule labeled by an isotope emitting positrons during its decay. The emitted positron is annihilated with a nearby electron, producing two photons traveling nearly colinearly in opposite directions as shown in Figure 2.1. The emitted photon pairs are detected by the PET detectors surrounding the subject. Figure 2.2 shows a PET detector pair A and B detecting a coincident photon pair. PET neutrino Positron-electron annihilation γ γ An unstable neutron deficit isotope photon p n n p p p n p n n p n p n ν neutron proton Proton Decay Istope stablized ∼180 ◦ Figure 2.1: Physical principle of PET imaging. An unstable neutron deficit isotope undergoes a proton decay producing a neutron, a positron and a neutrino. The two prod- ucts of the decay, positron and neutrino, are emitted from the isotope, and the emitted positron is annihilated with a nearby electron. Due to the positron-electron annihilation, two 511 keV photons are produced, and emitted in opposite directions of approximately 180 ◦ . 6 γ γ A B Δ 1 <T Δ 2 >T Pulse generator T : Coincidence timing window Coincident event Figure 2.2: The emitted photon pair from the positron-electron annihilation in Figure 2.1 is detected by PET detector pair A and B. The photon pair is called a coincident event when their detected time difference is smaller than the preset coincidence timing win- dow,T . detector hardware can detect the approximate arrival time difference, and determine a coincident event when their detected time difference is smaller than a preset timing win- dow. The PET data described in the next section is the collection of all those detected photon pairs determined to be in coincidence between any pair of detectors. 2.2 3D PET Data Structure 3D PET data is histogrammed [31] so that each data element, or bin, indicating a detector pair index, contains the total number of events detected within a tube connecting the detector pair as shown in figure 2.3. The tube in the figure is called a line of response (LOR). In the transverse view of the scanner, the detectors are paired in the way shown in Figure 2.4, so that the data collected through the parallel lines can be regarded as parallel 7 Line (tube) of Response (LOR) Detector element (Crystal) Figure 2.3: A line ( or tube) of response (LOR) in 3D PET acquisition mode between two detectors in two different rings. Each histogrammed projection data bin contains the number of detected coincident events within the line of response. projections. When the parallel projection data are arranged for all angles from 0 to180 ◦ , the data set is called a sinogram. Since a 3D PET scanner has multiple rings, as shown in Figure 2.5(a), we can have multiple sinograms acquired between different scanner ring pairs. So there are two types of sinograms in 3D PET data. As shown in Figure 2.5(a), one is the 2D direct sinograms, which are acquired from detector pairs in the same ring, and the other one is the oblique sinograms, which are acquired from detector pairs in two different rings. If there areN detector rings, there can beN 2 sinograms. Therefore, the 3D PET data size increases rapidly with the number of scanner rings [6]. To reduce the large size of 3D PET data, several adjacent lines of response are often combined in the axial direction as shown in Figure 2.5(b). Part of the data are also some- times discarded. Two terminologies, span and maximum ring difference (MRD), are used in this context [31, 40]. The MRD defines maximum limit of the allowed ring dif- ference so that the data acquired between two rings where their ring difference exceeds 8 Parallel projections Object Figure 2.4: The detectors are paired in the way to give parallel LORs for a fixed angle to form a parallel projection data. the MRD are discarded. The span determines the number of axial LORs which will be combined together as illustrated in Figure 2.5(b). A convenient illustration tool for these data reduction methods is the Michelogram illustrated in Figure 2.6 as an example. In Figure 2.6, the asterisk (*) mark in each square box stands for one sinogram between two rings where their ring indices are marked in horizontal and vertical axes. Then, to reduce the PET data size, the sinograms (asterisks) connected by the solid line are combined, which is called axial mashing or variable axial rebinning [15, 41, 109]. Here, the span is the sum of the numbers of combined LORs in odd and even planes, that is, 5 (2+3) in this example. The dotted line divides the 3D PET data into segments, each of which has the same ring difference marked Δ. An MRD of 12 was used in this example, so the sinograms marked as circled asterisks are discarded. In this example, by using span of 5 and MRD of 12, the total number of sinograms is reduced to 111 from 256 (16 2 ). Once all the data are arranged in the 3D PET sinogram format, the 9 2D direct sinograms Oblique sinograms Combining together: (Axial mashing) (a) (b) Oblique angle Figure 2.5: (a) The 3D PET sinogram is the combination of 2D direct sinograms and all oblique sinograms. Notice that the number of oblique sinograms for a fixed oblique angle is less than the number of 2D direct sinograms due to the finite axial aperture of the scanner. (b) Axial mashing is performed to reduce the 3D PET data size. The parameter ‘span’ determines the number LORs which will be combined. data can be parametrized by m(s,φ,z,Δ) where s and φ are the sinogram radial and angle indices, respectively; Δ is the ring difference index (or segment index); z is the axial plane index indicating each solid line in Figure 2.6 for each segment (or fixed ring difference Δ). As a special case, when Δ is zero, the data m are the stacked 2D direct sinograms. The number of planes in each segment in Figure 2.6 decreases as the ring difference increases because of the finite axial aperture of the scanner. Table 2.1 shows the 3D PET data size with typical MRD/span values for four Siemens microPET and clinical PET scanners. As shown in the table, the data size can be approximately 1 Gbytes for the High Resolution Research Tomograph (HRRT) [135]. 10 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Δ =−10 Δ =−5 Δ = 0 Δ = 5 Δ = 10 1 6 11 1 6 11 16 16 Figure 2.6: Example of the Michelogram which is useful for illustrating 3D PET data reduction by span and maximum ring difference (MRD). The Michelogram is shown for a 16-ring scanner, where the span is 5 (2 + 3), and MRD is 12 (13-1) so that the sinograms marked as circled asterisks are discarded. 2.3 Iterative PET Image Reconstruction PET image reconstruction methods can be classified into two categories. One is ana- lytical methods, and the other is iterative reconstruction methods. In analytical meth- ods [23, 25, 71, 95, 119], the line integral based data model is used, and image recon- struction is performed using analytical inversion formula, so the calculations are typi- cally very fast. However, the methods do not effectively account for the signal depen- dent noise, and the simple line integral model does not effectively account for several 11 microPET microPET Biograph ECAT Parameters F220 [129] Inveon [91] TPTV [57] HRRT [135] Sinogram rays (s) 288 128 336 256 Number of angles (φ) 252 160 336 288 Number of rings 48 80 55 104 MRD/span 47/3 79/3 54/5 67/3 Number of sinograms (z,δ) 1567 4319 1275 6367 Total size (MBytes) 434* 337* 275** 895** Table 2.1: Sinogram dimension and 3D PET data size for four Siemens PET scanners. Note that * indicates a sinogram bin size of 4 byte integers, and ** is a sinogram bin size of 2 byte integers. complicating factors such as positron range, noncolinearity of the photon pair, vari- ations in detector-pair sensitivity along the LORs, and the spatially variant detector response [76, 100, 107]. In iterative methods, a statistical noise model of PET data can be used, and all the complicating factors can be accounted for in the system model. As a noise model for PET data, the Poisson distribution is used since the photon counting statistics in the measured data can be well modeled as Poisson. Using the system model, and the noise model, the image reconstruction problem is translated into an estimation problem. Here we focus on the following log-likelihood function L(x|y) = X i [−(ˆ y i +b i )+y i log(ˆ y i +b i )] (2.1) wherey i , ˆ y i andb i are elements of the measured data vectory, the mean vector ˆ y ofy and the estimated background noise vectorb, respectively. x is the mean radioactivity (image) vector which is related to ˆ y by the system matrixP as ˆ y =Px+b. 12 The MAP reconstruction method [109], one example of the statistical approach, is to find a maximum a posteriori (MAP) solution as the maximizer of the following objective function: ˆ x(y) = argmax x≥0 L(x|y)−βφ(x) (2.2) whereφ(x) is a penalty function encouraging local smoothness of the estimated image, andβ is the smoothing or hyper parameter. This problem can be solved using an iterative optimization technique. As a special case of (2.2), when the hyper parameterβ is zero, the problem becomes maximum likelihood estimation. In that case, the following closed form update equation can be obtained using expectation-maximization (EM) algorithm [68, 125] x (k+1) j = x (k) j P i p ij X i p ij y i P l p il x (k) l +b i (2.3) wherex (k) j is the image element estimated at the k-th iteration, andp ij is an element of the system matrixP . Equation (2.3) requires one forward projection (multiplication by P : the third summation term in (2.3)) and one back projection (multiplication by the transpose ofP : the second summation term in (2.3)). Most iterative algorithms need the forward and backprojections for each iteration, which are computationally expensive, particularly for 3D PET, due to the large data size as shown in Table 2.1. Several approaches have been proposed to address the computation issue in the fully 3D iterative reconstruction methods. The first approach is simply to utilize powerful computing devices such as parallel computing machines or graphics cards [13,102,122]. The second approach is to use a part of data at each iteration to reduce the computation. 13 2D direct sinograms 3D PET data Redundant part Fourier Rebinning 2D stacked direct sinograms in the same dimension Slice-by-slice 2D reconstructions Figure 2.7: Fast 3D PET data reconstruction using Fourier rebinning. 3D PET data is rebinned to 2D stacked direct sinograms so only a slice-by-slice 2D reconstruction is required after Fourier rebinning. Fourier rebinning is an algorithm to reduce the dimensionality of 3D PET data. Hudson and Larkin in [54] originally proposed the ordered subset EM (OS-EM) algo- rithm. Since then, other fast subsetting algorithms have been proposed [11,16,38]. How- ever, these methods do not guarantee global convergence to an optimal solution. This problem was remedied by using a relaxed stepsize parameter [1,11,38,63,64,97,98], or incremental EM type algorithms [2,45,50,55,56,61,99]. A third approach for reducing computational cost is to design a fast projector [27,51,60,80]. This approach is effective due to the fact that most of the computation time for fully 3D PET iterative reconstruc- tion is consumed by forward and backprojections. The method in Chapter 3 belongs to this category. As an alternative way of addressing this computation issue, we briefly review in the next section a hybrid method where the 3D PET datam(s,φ,z,Δ) are first rebinned to 2D stacked sinogramsm(s,φ,z,Δ = 0) using Fourier rebinning, and then a slice-by- slice 2D iterative reconstruction is applied to the reduced data (See Figure 2.7). The method has been very successful due to its reasonable image quality and very good computational efficiency. 14 2.4 Fourier Rebinning for Fast 3D PET Image Recon- struction The 3D PET data consist of stacked 2D direct sinograms and oblique sinograms. How- ever, in principle, all the oblique sinograms are redundant in the sense that the 3D image can be reconstructed slice-by-slice only using 2D direct sinogram data. Here, the 2D direct sinograms are the subset of the 3D PET data m(s,φ,z,Δ) for which Δ is zero. However, the oblique sinograms help to improve the signal-to-noise-ratio (SNR) in the reconstructed images. Defrise et al proposed the Fourier rebinning (FORE) method [29] which can reduce the 3D PET data to 2D stacked sinograms while retaining most of the SNR advantage of 3D PET data. After the Fourier rebinning is applied to 3D PET data, the reconstruction problem is reduced to a 2D slice-by-slice image reconstruction as shown in Figure 2.7. Since the data reduction step using FORE is done very quickly using the fast Fourier transform (FFT), the method has significant computational effi- ciency. The Fourier rebinning algorithm is derived based on an ideal line integral model [32]. Using the notation in [29, 32], the measured 3D PET sinogram data m(s,φ,z,Δ) is expressed as a line integral of an object functionf(x,y,z) through an LOR as follows m(s,φ,z,Δ) = Z ∞ −∞ dtf(scosφ−tu x ,ssinφ+tu y ,z +tu z ) (2.4) whereu = (u x ,u y ,u z ) is a unit vector through the LOR shown in figure. 2.8, expressed as u =   −sinφ,cosφ, Δ 2 √ R 2 −s 2 q 1+ Δ 2 4(R 2 −s 2 )   where s and φ are the radial and angular coordinate parameters, respectively; z is the axial coordinate of the midpoint of each LOR; andΔ is the ring difference, as illustrated 15 x y θ φ t (a) (b) δ = tanθ = Δ 2 √ R 2 −s 2 z LOR s Δ R 2 √ R 2 −s 2 Figure 2.8: (a) Top and (b) side view of a cylindrical 3D PET scanner. in Figure 2.8. Using this data model, the derivation of the Fourier rebinning algorithm is obtained from the following slightly modified data parametrizationp(s,φ,z,δ) [29,32]. p(s,φ,z,δ) = m(s,φ,z,Δ = 2δ √ R 2 −s 2 ) √ 1+δ 2 = Z ∞ −∞ dtf(scosφ−tsinφ,ssinφ+tcosφ,z +tδ) (2.5) whereδ is the tangent of the oblique angleθ (see Figure 2.8). We will callp “pseudo” 3D sinogram since it is slightly different from the parametrization for the measured sinogram m. The pseudo sinogram can be considered to be a version of the measured datam scaled by 1/ √ 1+δ 2 whenR >> s since thenδ ≈ Δ/2R. This is a reasonable assumption when the image field-of-view (FOV) is small, but when the image FOV is very large, the assumption may not be accurate. The implications of this approximation will be addressed on our work in the next chapter. Forδ = 0,p(s,φ,z,0) is equivalent to m(s,φ,z,0) and represents the stacked 2D sinograms for the 2D transaxial slices. From the parametrization, the aim of the Fourier rebinning algorithm is to find a mapping between oblique sinogram datap(s,φ,z,δ), whereδ6= 0, and the stacked 2D sinogram 16 data p(s,φ,z,0) so that the oblique sinogram data are converted and added up to gen- erate the corresponding stacked 2D sinogram data. Through this processing, we can reduce the 3D PET data to 2D stacked sinogram data, while retaining the SNR benefit of 3D PET data by fully utilizing the oblique sinograms. Historically, the first mapping was proposed in [29] based on the frequency distance principle [39] as follows P(ω s ,k,z,δ)≈P(ω s ,k,z− kδ ω s ,0), (2.6) which is an approximate mapping. Here, the functionP is the 2D Fourier transform ofp in (2.5) with respect to (s,φ),ω s is the frequency index fors, andk is the Fourier series index in the angular direction (the Fourier transform is discrete with respect to angleφ because p is a periodic function in φ). The data reduction method using this mapping is called Fourier rebinning (FORE). Whereas FORE is approximate, Defrise et al also derived the following exact rebinning equation by further taking Fourier transform in thez direction [32]. ℘(ω,k,ω z ,δ) = exp −ikarctan δω z ω ℘   ω s 1+ δω z ω 2 ,k,ω z ,0   (2.7) where ℘(ω,k,ω z ,δ) is the 3D Fourier transform (FT) ofp(s,φ,z,δ) with respect to s, φ and z at a fixed δ; ω z is the axial frequency variable; This exact equation is used in the exact rebinning method FOREX [70]. It requires taking Fourier transforms through the axial variable z. However, due to the limited axial aperture of a PET scanner, the complete data set is not available, in other words, the data is truncated in the axial direc- tion for oblique sinograms as shown in Figure 2.5(a). Therefore, we need to estimate the truncated data before taking the Fourier transform. Liu et al used inverse rebinning 17 based on (2.7), to estimate the truncated data from the 2D direct sinograms in measured 3D PET data. Due to this extra step and the extra Fourier transform taken in the z direction, the computation time for the exact rebinning method is larger compared to the previously mentioned FORE algorithm [70], and estimation of the truncated data may cause some artifacts [10]. To avoid the missing data problem, Defrise and Liu proposed another exact method FOREJ based on John’s equation [33]. This method does not suf- fer from the truncated data problem in FOREX, but still gives comparable accuracy to FOREX while reducing computation time. All of the work in this thesis is related to the Fourier rebinning concept. In Chap- ter 3, we use inverse Fourier rebinning mapping for designing a fully 3D fast projector since the inverse rebinning can be implemented cost-effectively using the FFT. In Chap- ter 4 and 5, we extend Fourier rebinning to 3D time-of-flight PET data reconstruction. We develop rebinning algorithms to reduce the large 3D TOF PET data to its lower dimensional versions. 2.5 Time-of-Flight (TOF) PET In PET, the time-of-flight (TOF) is the time difference between detection of the two photons produced by the positron annihilation. High resolution measurement of this time difference would allow one to determine the precise location at which the anni- hilation occurred. Although the idea of utilizing this TOF information appeared in the 1960’s [5, 14], lack of fast scintillators has limited building practical TOF PET systems until the recent development of fast detectors such as LSO [77, 85, 90] and LaBr 3 [66, 117]. Table 2.2 summarizes some important properties of various scintillation materials [21, 72, 78, 84, 101, 131]. Here, the density of the scintillator is directly related to the 18 Density Light output Decay Scintillators (g/cm 3 ) (NaI = 100) time (ns) Nal(TI) 3.67 100 230 CsF 4.64 5 5 BaF 2 4.88 28 0.8 BGO 7.13 15 300 LSO 7.40 70 40 LaBr 3 5.30 160 25 Table 2.2: Property comparison of various PET scintillators. stopping power of high energy photons, and therefore is an important factor determining the necessary crystal thickness which affects the intrinsic spatial resolution of the PET scanner. The high light output enables us to couple more crystals to one photomultiplier tube (PMT) so we can reduce the number of PMTs required. For example, with BGO, the block detector uses up to 16 crystals per PMT, but over 144 crystals can be coupled to a PMT in the LSO case [78]. Finally, the decay time determines coincidence resolving time, and is the key factor enabling TOF PET system feasibility. The first choice of scintillator in the 1970’s was NaI(TI) because of its high light out- put, but it is very difficult to manufacture since it is hygroscopic [78]. In the early 1980’s, CsF and BaF 2 were the detectors used for TOF PET systems [79, 93, 128] because of their short decay time. However, they were not dense enough to give good spatial reso- lution. So BGO detector-based 3D PET system [20] were dominant by the early 990’s. The performance of systems using BGO detectors was limited by dead time, randoms, and scatter because of the scintillator’s long decay time. The more recently developed LSO detector [90] has the advantages of fast decay time, high density and fairly good light output. The most important benefit of systems using this scintillator is that this can also generate TOF information [19, 89] without introducing disadvantages such as the low density of CsF or BaF 2 . A commercial TOF PET scanner using LYSO scintillator, 19 a slightly modified version of LSO, is now available, and achieves better image quality by utilizing the TOF information [118]. LaBr 3 is also used for TOF PET systems since it has very high light output and short decay time although it has a lower density than LSO [66, 117]. Most of the TOF PET systems developed in 1980’s were 2D systems where each transverse plane is reconstructed using mostly 2D analytical reconstructions based on line integral model [124, 127]. In principle, localization of the point of annihilation, that is, image reconstruction, can be direct if we can measure the time difference accu- rately. However, there is measurement uncertainty in the TOF information which has been narrowed down to around 70 mm at full-width-at-half-maximum (FWHM) with the current technology [19, 117]. It is reported that this timing resolution can still help to increase SNR of the reconstructed image when incorporated into the image reconstruction [52, 89, 92, 124, 127]. As indicated in [89], and described analytically in [124, 127, 132], the SNR improvement in TOF PET image reconstruction is deter- mined by the TOF timing resolution. Conceptually, as Figure 2.9 shows, when there is no TOF information available, we just assume the positron annihilation happened some- where along the line of response, so the error in the measured data is distributed along the whole LOR when backprojected. When TOF information is available, the uncer- tainty region is reduced by the weighted backprojection using the TOF kernel function, so the the error propagation is limited proportionally to the width of the TOF kernel function. The TOF PET data can be modeled as a line integral weighted by the kernel function through a line of response as shown in Figure 2.10 [124,133]. Based on this data model, the image reconstruction problem is to recover the object functionf from the measured 20 No information about where the activity happened Measuring certainty probability distribution: constant Line of Response (LOR) (a) Conventional case Actual positions where coincident events happened Measuring certainty probability distribution Line of Response (LOR) (b) TOF PET case Figure 2.9: Time-of-flight (TOF) information can be utilized for better image recon- struction. In (a), when backprojected, the uncertainty in the measured data is distributed along the whole LOR, but in (b), uncertainty is reduced by weighted backprojection with the TOF kernel function. data p as illustrated in Figure 2.10. In [133], this is done by taking the Fourier trans- form of the data model under the assumption that the TOF kernel function is spatially invariant: P(ν,φ) =F(ν)H(ν·d) (2.8) whereν is the 2D frequency vector for spatial coordinatesu in Figure 2.10. This means that, in principle, we can recoverF at a fixed angleφ, which the TOF data in other angles are redundant. In order to fully use the TOF data, for maximum noise suppression, a filter W(d,ν) is applied and all the data is summed with respect to φ to generate following equation: Z dφW(d,ν)P(ν,φ) =F(ν) Z dφW(d,ν)H(ν·d) (2.9) 21 TOF Kernel function h True annihilation position Detector Detector f(u) φ u=(x,y) p(u,φ)= Z f(u+ld)h(l)dl d = (cosφ,sinφ) Line of Response (LOR) Figure 2.10: TOF PET data can be modelled using a line integral weighted by a TOF kernel function which is a probability density function of the measurement error asso- ciated with the time difference of two detected photons at the detector pair. It is often assumed spatially invariant, and approximated as a Gaussian distribution with standard deviation determined by the timing resolution of the system. Then,F(ν) can be recovered by: F(ν) = Z dφG(d,ν)P(ν,φ) (2.10) where G(d,ν) = W(d,ν) R dφW(d,ν)H(ν·d) . (2.11) The problem findingF reduces to finding the filterG that optimizes in the reconstruc- tion. Watson suggested several filters in [133], and noted that (2.10) is equivalent to the “Most likely position (MLP)” case indicated in [124] when G = 1, and “Confidence weighting (CW)” whenW = G, which gives reasonably good results. This choice was also suggested in [127] as an optimized weighting function used in the backprojection through an LOR. Notice the different parametrization using the 2D space variableu and 22 Detector Detector Δt D Imaging object TOF Kernel function h Figure 2.11: TOF SNR gain is approximately expressed by the FWHM of the TOF kernel function and object size. angleφ in Figure 2.10, in contrast to the Radon transform parametrization for conven- tional sinogram data used in other work [28]. To quantify the benefit produced by using the TOF information, the variance ratio between TOF (V TOF ) and non TOF (V nonTOF ) was derived analytically by an assuming infinitely large uniform image source [127]: V TOF V nonTOF = 0.62 D Δt (2.12) where D is the object diameter, and Δt is the FWHM of the TOF kernel function as shown in Figure 2.11. A similar equation is also derived in [14], and recently confirmed by Fisher information matrix based analysis in [132]. This analysis tells us that, if we study a 30cm diameter patient, we can achieve approximately a 3 times variance reduction when Δt = 7cm. All of the work described here is for 2D TOF PET. There is very little literature dealing with 3D TOF analytical reconstruction, Mallon and Granget [81] extend the fully 3D PET analytical reconstruction [36] to 3D TOF PET, but with very preliminary 23 results. There have been several approaches to reconstructing 3D TOF PET data itera- tively using list mode data [47, 117], histogrammed data [19], or by backprojecting the TOF PET data to image space [82]. However, we also expect a heavy computational cost in fully 3D TOF PET iterative reconstruction since the data size will be even larger than 3D PET data due to the extra TOF information. Defrise extended the rebinning concept to the 3D TOF PET data [28, 35]. Here we also develop rebinning methods for data dimensionality reduction using mapping equations derived on a new unified frame- work. We can also use the mapping equations to develop a fast projector for fully 3D TOF PET iterative reconstruction as we did for fully 3D PET iterative reconstruction. 24 Chapter 3 Fast Iterative 3D PET Image Reconstruction 3.1 Introduction Iterative reconstruction methods based on a Poisson noise model have been widely demonstrated to provide superior image quality compared to analytic reconstruction methods [22, 44, 59]. However, fully three dimensional (3D) iterative reconstruction methods can become prohibitively expensive since they need multiple passes going through the large data, each of which requires a minimum of one forward and one back projection. Consequently, the availability of fast projection methods would facilitate more routine use of fully 3D iterative reconstruction methods. Here, as our first work, we address the problem of iterating between the fully 3D data and the 3D image volume while using Fourier rebinning [32] concept to reduce the computation cost. Several fast projection algorithms have been published for the purpose. Matej et al in [80, 86] proposed Fourier-based projectors in which the fast Fourier transform (FFT) algorithm is used in combination with the Fourier slice theorem to transform the computation from the spatial to the frequency domain. The major limitation to this approach is that it is not well suited to use with algorithms that iterate over subsets of the data. This other method, which achieves a similar speed up, is based on a hierarchical decomposition that exploits, so called, bow-tie property that the angular bandwidth of the sinogram is proportional to the image size [112]. This decomposition 25 was used to develop fast forward and back projection operators [7–9] that can be applied to fully 3D PET, but to date published fully 3D applications are limited to cone beam tomography. A third approach to development of fast projectors is to use one of the Fourier rebinning methods. Here we use an exact inverse rebinning method similar to that described in [32, 70]. As we show below, it is sufficient for our purposes to define an inverse rebinning operator from which we develop both the forward and back pro- jection operators. In contrast to exact rebinning, exact inverse rebinning can be imple- mented relatively efficiently. The inverse rebinning results assume an idealized line integral model rather than a more accurate geometric model based on solid angle sen- sitivity considerations [108]. However, the projectors based on inverse rebinning can still be combined with our factored system model [109], which includes a blur kernel in sinogram space that models shift-variant detector response. In this way we can largely retain the resolution recovery properties of our accurate 3D system model. While our approach differs in its use of inverse rebinning, a number of other related studies have also demonstrated the particular importance of detector response modeling in resolution recovery [3, 4, 44, 73, 74, 104, 111, 120, 123, 137]. We evaluate the fast projector pair in comparison to our 3D geometric factored system model which has been widely used in small animal PET systems [26]. We also investigate the impact of modifying the inverse rebinning procedure to correct for nonuniform radial sampling and nonconstant oblique angles in measured sinograms from a cylindrical PET scanner. Finally, we illustrate performance with simulated data and an in vivo 18 FDG monkey brain scan. 26 3.2 Background 3.2.1 Factored System Model We focus on a maximum a posteriori (MAP) reconstruction based on the Poisson model [76, 100] briefly reviewed in the previous chapter. The major cost of the reconstruction comes from forward and back projection operations. The projection or system matrix P models the sensitivity of each detector pair to positron emissions from each voxel. In fully 3D systems this matrix is huge (> 10 14 elements for the microPET Focus 220 small animal scanner [129]) but can be stored efficiently through the use of geometrical symmetries and sparse matrix structures in combination with the factorization [88,108]: P =P det.sens P det.blur P attn P geom P positron (3.1) whereP det.sens is the detector-pair sensitivity matrix; P det.blur is a local sinogram blur- ring operator that models detector response effects resulting from photon pair nonco- linearity, crystal penetration, inter crystal scatter, and the effects of the block detector geometry; P attn is a diagonal matrix containing the attenuation factors for each LOR; P geom is the geometric sensitivity matrix representing the probability that an emission from each voxel is detected by the detector pair associated with each line of response (LOR) in the absence of attenuation and noncolinearity effects; and P positron is a local image-space blurring operator that models positron range effects [110]. Symmetries and sparse structures reduceP geom to a manageable size, e.g.∼ 500MB for the fully 3D microPET Focus 220 [129]. However, the geometric forward projection (i.e. multiplication by P geom ) and back projection (multiplication by P T geom ) are time consuming, requiring a mapping between each sinogram element and all those voxels that intersect the corresponding LOR. Here we propose an alternative form of the 3D 27 geometric projectorP geom and its transpose based on the exact inverse rebinning oper- ator. By retaining other factored system models as are, we can largely preserve the resolution recovery properties of the original system. 3.2.2 Two Complicating Factors Encountering When Using Fourier Rebinning Here, we want to utilize the exact inverse rebinning to design a fast projector which can replace P geom in (3.1). The inverse rebinning equation (2.7) gives an exact mapping between 2D direct sinogram and 3D sinogram. However, there are two complicating factors that should be considered when applying the rebinning methods to measured sinogram data in cylindrical PET scanners [29, 32, 33, 70]. For the explanation, same notations from previous chapter will be continuously used. First, in cylindrical PET scanners with discrete detectors, the LORs between detector pairs do not correspond to uniform samples of m(s,φ,z,Δ) in radial (s) direction. Sample spacing decreases towards the edge of the transaxial field of view as a result of the curvature of the detector ring. Consequently, “arc correction” should be applied to the measured data to achieve uniform radial sampling. To indicate the extent of this sampling nonuniformity, we show the intersample interval in mm as a function of radial index in figure 3.1 for the microPET Focus220 scanner. A further complication is that for fixed (φ,z,Δ), the LORs for different s also have different δ values since δ = Δ/2 √ R 2 −s 2 . Consequently, the oblique angle θ = arctanδ of an LOR for an oblique sinogram increases as we move away from the center of the field-of-view (FOV) in the radial direction. The oblique angle variation will be small for a small FOV because we can assumeδ ≃ Δ/2R whens≪ R. In this case,p(s,φ,z,δ) can be approximated as the 3D sinogram datam(s,φ,z,Δ) scaled by the constant √ 1+δ 2 [32]. However, for larger fields of view, the variation is significant 28 160 180 200 220 240 260 280 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 Sinogram radial index from the center(145−th sample) Radial sampling interval(mm) Figure 3.1: LOR sample interval variation as a function of sinogram radial index for the microPET Focus 220 scanner as shown in figure 3.2 for the full field of view of the microPET Focus220 scanner. In this case we need to account for the variation in δ with respect to s by interpolating between oblique sinograms as described in Chapter 3.3.3. 3.3 Fast 3D Projector Using Exact Inverse Rebinning The projection operatorP geom in (3.1) can be approximately factored into the product of three operators: (i) a 2D projector that maps the 3D image data into a stack of 2D sino- gramsp(s,φ,z,0); (ii) the inverse rebinning operator that maps these 2D sinograms into the full 3D pseudo sinogramsp(s,φ,z,δ); and (iii) a post correction step that accounts for differences between pseudo sinograms p(s,φ,z,δ) and the measured 3D sinogram datam(s,φ,z,Δ). 29 160 180 200 220 240 260 280 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 δ = tanθ Sinogram radial index from the center(145−th sample) Figure 3.2: The variation of δ (the tangent of oblique angle θ) as a function of radial index for oblique sinograms for the microPET Focus 220 scanner. Each curve shows the variation inδ for a fixed ring difference, with a maximum ring difference (MRD) of 47 and a span of 3 We will investigate three different forward projectors based on exact inverse rebin- ning, as illustrated in figure 3.3. Case A uses a 2D geometric projector (P2DG) to map the 3D image into a set of stacked 2D sinograms. This method includes arc correction (AC) and inverse arc correction (IAC) to map between the nonuniform sampling ins for the geometric projector and the uniform sampling required when applying the Fourier rebinning equations using FFTs. Case A also uses axial resampling (ARS) to interpo- late between oblique pseudosinograms to account for variations in δ as a function ofs in the measured sinograms. The final normalization (N) step accounts for differences in geometric sensitivity between the measured and pseudo sinograms as a function of radial position s and oblique angle θ. Case B is the same as Case A except that we 30 P2DG AC IRB IAC ARS N P2DG IRB N P2DL IRB IAC ARS N (a) Case A (b) Case B (c) Case C (i) 2D Projection (iii) Post Corrections (ii) IRB Figure 3.3: The 3D geometric projectorP geom can be approximated using inverse rebin- ning projectors via the following three steps: (i) 2D projection, (ii) inverse rebinning (IRB), and (iii) post corrections. Three different variations based on IRB are inves- tigated here: Case A uses a 2D geometric projector with corrections for nonuniform sampling ins andδ, Case B uses the same 2D geometric projector but without the cor- rections, and Case C is the same as Case A but uses a 2D line integral projector in place of the 2D geometric projector. ignore nonuniform sampling ins andδ. We would expect this case to produce reason- able results for a FOV that is small relative to the radius of the detector ring. Case C replaces the 2D geometric projector with a line integral projector so that we directly compute uniform samples ofp(s,φ,z,0) in the first stage and therefore do not need the initial arc correction step. We then apply the same postcorrections as for Case A. We now describe each of the three stages in figure 3.3 in detail. 3.3.1 2D Projection: 3D Image to 2D Direct Sinograms PET images are typically reconstructed as 2N − 1 transaxial sections for an N-ring scanner. The 2D direct sinograms for the N odd numbered planes are formed directly from LORs corresponding to a single detector ring. The 2D sinograms for theN−1 even numbered planes are formed by adding LORs from adjacent rings. Combining adjacent sinograms in this way corresponds to using a span of 3 and using only the direct 2D 31 sinograms is a maximum ring difference (MRD) of zero. All of the fully 3D results presented below use minimal axial rebinning with a span of 3 and a maximum ring difference ofN−1. For the 48 ring Focus220 scanner referred to in Figure 3.1 and 3.2, this reduces the number of sinograms from 2,304 to 1,567 and the number of different oblique angles from N − 1 = 47 to N/3 = 16. We note that our fully 3D geometric form ofP geom [108] models this rebinning process by summing the appropriate rows of the full matrix. To compute the stacked 2D sinograms from the 3D image we can use our P geom geometric matrix computed for the 2D case, i.e. span of 3 and MRD of 1. The 2D geometric calculation includes considerations of geometric response based on the solid angles subtended by each voxel at each detector pair as described in [108]. Further- more, since the calculations are based on the detector ring geometry, the resulting 2D sinograms are nonuniformly sampled in s and must be resampled before applying the inverse rebinning operator in the next step. We use a cubic B-spline interpolation to perform arc correction for Case A, but also investigate performance in the Case B where this effect is ignored. A third alternative, Case C, is to directly use a 2D line integral based projector to compute uniformly sampled 2D sinograms so that arc correction is not required.Furthermore, use of the true line integral rather than the 2D geometric pro- jection may lead to improved performance since the resulting sinograms should more closely obey the consistency conditions for true line integrals [96] than those computed using the geometric projector. For the 2D line integral projector, the ‘radon’ function in Matlab was used. 32 3.3.2 Inverse Rebinning: 2D Direct Sinograms to Pseudo 3D Sino- gram For inverse rebinning to map from 2D to 3D sinograms, we use a modified version of (2.7). Using the shift property of the FT and taking the inverse FT with respect to the angular frequency indexk, the exponential term on the right hand side drops out leaving the 2D FT relationship: P(ω,φ,ω z ,δ) =P   ω s 1+ δω z ω 2 ,φ−arctan δω z ω ,ω z ,0   (3.2) where P(ω,φ,ω z ,δ) is the 2D FT of p(s,φ,z,δ) with respect to s and z. The key advantage of this form is that it facilitates the use of inverse rebinning in combination with ordered-subset methods. The additional cost associated with working withφ rather thank is that we need to perform 2D interpolation in(φ,ω) rather than 1D interpolation inω as would be the case when using (2.7). However this cost is offset by the reduced cost of the 2D FFT relative to the 3D FFT. A further potential advantage is that there is no need to resample in angle in order to produce a highly factorizable number of samples inφ for efficient computation of the FFT. Based on (3.2), the inverse rebinning operator, denoted B δ IRB , can be written in a factored matrix form as: B δ IRB ≡S δ X φ B φ F −1 2D I φ,δ ! F. (3.3) The operator B δ IRB calculates oblique sinograms p(s,φ,z,δ) at a fixed δ. The block matrixF computes the 2D FFT with respect to (s,z) for all anglesφ to yield the com- plete P(ω,φ,ω z ,0) data in (3.2). The FFTs P(ω,φ,ω z ,δ) of the oblique sinograms at fixed δ are obtained from P(ω ′ ,φ ′ ,ω z ,0) using (3.2) where ω ′ = q 1+ δωz ω 2 and 33 I 1,δ F −1 2D I 2,δ I φ m ,δ X ... B 1 B 2 B φ m F −1 2D F −1 2D F p(s,φ,z,0) √ 1+δ 2 p(s,φ,z,δ) S δ Figure 3.4: Block diagram of inverse rebinning operator B δ IRB which computes the oblique sinograms at fixed oblique angle indexδ from 2D direct sinograms. The dotted line box indicates calculation for eachφ whereφ m is the largest index forφ. φ ′ = φ− arctan δωz ω . The matrices I φ,δ perform the interpolation with respect to (ω ′ ,φ ′ ) to compute P(ω,φ,ω z ,δ) from P(ω ′ ,φ ′ ,ω z ,0). The 2D inverse FFT operator F −1 2D with respect to (ω,ω z ) yieldsp(s,φ,z,δ). The concatenation operatorB φ is used after the inverse FFT so that the summation overφ combines the components of the 3D sinogram for all values of φ at a fixed value of δ. The diagonal matrixS δ represents a scaling by √ 1+δ 2 to map from the pseudo to measured sinogram according to (2.5). A block diagram for the inverse rebinning operator is shown in figure 3.4. The decomposition with respect to φ in the dotted line boxes indicates that the method can be readily applied to ordered subset algorithms [1, 54]. Since the inverse rebinning operator for eachδ and each dotted line box in figure 3.4 can be computed independently, the inverse rebinning operation can also be easily parallelized. Note that the 2D FFT operationF in figure 3.4 does not need to be computed for eachδ. Once it is calculated, it can be used for allδ. Finally, we note that the direct planes are passed directly to the output in figure 3.4 without the need to recompute them using inverse rebinning. 34 We used bilinear interpolation via precomputed lookup table forI φ,δ in (3.3), which produced reasonably accurate results as reported in Chapter 3.4. Higher order inter- polants may improve accuracy but would impact overall runtime. The tables for the bilinear interpolation and its adjoint (for back projection) have relatively small storage requirements, e.g.∼10Mbytes for the microPET Focus 220 scanner. The computational cost of inverse rebinning is dominated by the 2D FFTs and 2D interpolation. The cost is almost linear in the number of oblique angles as determined by the maximum ring difference and the span [31, 32, 40]. The computational complexity C can be approximated as C =N ob {N v (O s logN s +O z logN z )+2LN v +N v } +N v (O s logN s +O z logN z ) (3.4) where N v = N s N φ N z ; N s , N φ , and N z are the number of sample points through the radial, angular and axial directions respectively; N ob is the total number of oblique angles; O s and O z are the scaling factors for FFT complexity. L is the length of the 2D interpolation kernel with L = 4 for bilinear interpolation. Since we interpolate values in the complex frequency domain, 2L operations are required. 3.3.3 Post Corrections: Pseudo 3D Sinogram to 3D Sinogram The final step after calculating the pseudo 3D sinogram data using inverse rebinning is to map them to the measured sinogram space. This involves resampling or interpolation with respect to the radial variable s and the tangent δ of the oblique angle, followed by normalization to account for differences in geometric sensitivity between the line integral and our fully 3D geometric model. In Case B in figure 3.3 we ignore the need for resampling in s and δ and simply apply the normalization step. In the other two cases, A and C, we first apply the two resampling procedures. Both are performed using 35 δ s δ s IAC δ s ARS p(·,φ,z,·) p(·,φ,z,·) m(·,φ,z,·) Extra oblique samples for ARS (a) (b) (c) Axial angular sampling interval δ max Figure 3.5: Illustration of oblique sinogram resampling procedure: (a) uniform radial and axial angular sampling pattern in s and δ after applying IRB to compute oblique sinograms from the 2D direct sinograms; (b) radial and axial angular sampling pattern after IAC; (c) radial and axial angular sampling pattern after ARS and IAC. a 1D cubic B-spline interpolation procedures, first in s then in δ, as we now describe and illustrated in figure 3.5. After inverse rebinning, we have samples of the functionm(s,φ,z,Δ) uniformly in s and δ = Δ/2 √ R 2 −s 2 as illustrated in figure 3.5a. We use the same 1D cubic B spline interpolant for each oblique sinogram to resamplem(s,φ,z,Δ) nonuniformly in s. The resulting sample locations match the actual sample locations for the measured LORs as illustrated in figure 3.5b. As with the earlier 2D interpolants, we compute and store a look-up table to define the interpolation procedure. We then resample in oblique angle so that we have m(s,φ,z,Δ) uniformly sam- pled in Δ = 2δ √ R 2 −s 2 rather than in δ. We refer to this procedure as axial angular resampling (ARS). Since δ increases withs for fixed Δ, the oblique angles for large s are larger than their nominal values at the center of the field of view, as shown in figure 3.2. Consequently we need to compute additional oblique sinograms using IRB for val- ues ofδ larger than the valueδ max , which corresponds to the maximum ring difference. While it is also possible to increase the accuracy of the interpolation by sampling the sinograms more densely inδ, i.e. by computing sinograms at more oblique angles using IRB, the cost will increase almost linearly with this denser sampling. Consequently, we maintained the same sampling interval in δ but computed sufficient additional oblique 36 sinograms for δ > δ max to interpolate the largest value ofδ required for the maximum ring difference, as illustrated in figure 3.5a and c. We note that the number of oblique sinograms is actually a function of θ, since the maximum ring difference occurs only between the top and bottom most ring of the scanner. Consequently the additional num- ber of sinograms required for axial angular resampling is less than might be inferred from Figure 3.2 and 3.5c. Finally, for all three cases in figure 3.3, we apply normalization at the last step to account for reduced geometric sensitivity as the oblique angle increases. This drop off in sensitivity is caused by a reduced solid angle. While it is modeled explicitly in our calculation of the 3D geometric projector P geom [108], it is not considered in the inverse rebinning projector, since inverse rebinning is based on an ideal line integral model. The normalization factors were computed using a uniform cylindrical phantom covering the entire FOV . After computing forward projections using both the inverse rebinning projector and the 3D geometric projector, we obtain the normalization factor as a function of oblique angleθ and radial displacements by taking ratios between the two sinograms. 3.3.4 IRB Projector Pair: From 3D Image to 3D Sinogram and Back By combining all three steps described above, we can construct the IRB projector which maps a 3D image into 3D sinogram data as illustrated in figure 3.3. All three cases can be described by the factored matrix: y =N X δ ′ A δ ′ 2 T δ ′C δ ′ 3 X δ A δ 1 C 2 B δ IRB !! C 1 P 2D x (3.5) 37 where x is the 3D image vector; y is the vector form of the 3D sinogram data; B δ IRB is the inverse rebinning operator (IRB) in (3.3); P 2D is the 2D (geometric or line inte- gral) projector; C 1 is the arc correction operator (AC); C 2 is the inverse arc correction operator (IAC); C δ ′ 3 is the axial angular resampling operator (ARS); T δ ′ is a trunca- tion matrix which zeroes the unobserved or missing data in the oblique sinograms as described below;A δ 1 andA δ ′ 2 are the concatenation operators which, when summed over δ and δ ′ , assemble all oblique sinograms into a single vector; N is the diagonal nor- malization matrix; and δ ′ is the index which corresponds to Δ in m(s,φ,z,Δ). The operators C i (i = 1,2,3) become identity matrices for those cases in figure 3.3 where that resampling operator is omitted. The truncation matrixT δ ′ is a diagonal rectangular matrix with one’s on the diagonal which has the effect of truncating the input vector to remove those elements of the sinogram computed using IRB but which are not observed in the measured data because of the finite axial extent of the scanner. The back projection operator is obtained simply by transposing the whole operator in (3.5): ˆ x =P T 2D C T 1 X δ ′ X δ B δ IRB T C 2 T A δ 1 T ! C δ ′ 3 T T T δ ′A δ ′ 2 T ! Ny. (3.6) It is interesting to note that the backprojector in (3.6), which is the adjoint of the inverse rebinning procedure, is not the equivalent Fourier rebinning procedure. In exact rebin- ning we would need to perform reprojection to compute the missing oblique data, while in (3.6) the missing data remain zero through the transpose of the truncation operatorT δ ′ which simply adds zeroes to produce arrays of the correct dimension before applying further operations. An equivalent observation was noted by Matej et al. with respect to the use of the Fourier slice theorem for fast iterative reconstruction [80]. 38 As shown in figure 3.4 and 3.5, the inverse rebinning operator and all corrections are performed independently with respect to the sinogram angular variableφ, so that we can apply ordered subset methods [1, 54] to this projector for further speed up. However, before inverse rebinning can be applied, we must first compute the stacked 2D projec- tions. Therefore in any ordered subset method all 2D projections must be recomputed at each subiteration. While we anticipate that it is possible to achieve speedup using subset methods, the number of subsets into which the data is decomposed should not be too large, because of the fixed cost of computing these 2D projections at each iteration. Note that the current 2D blur kernel in (3.1) cannot be used for an ordered subset based method when the subsetizing is performed in the azimuthal angle direction. Instead, other approximated blur kernels [4, 104] which are not a function of azimuthal angle can be used. 3.4 Results 3.4.1 Computation Cost and Accuracy We now compare the new inverse rebinning projectors with the 3D geometric projector P geom in terms of accuracy and run time. All calculations were based on the microPET Focus 220 small animal scanner [129] whose specifications are listed in Table 3.1. To generate the fully 3D sinogram data we used a maximum ring difference of 47 and a span of 3 for a total of 1,567 sinograms. To study the effects of image size on projector accuracy we computed 3D sinograms from a256×256×95 digital 3D Hoffman brain phantom image scaled to two different transaxial voxel sizes: 1) 0.8mm (large FOV) and 2) 0.4mm (small FOV). The larger phantom covers the full transaxial FOV of the scanner. figure 3.6 shows sinograms for the two cases for a ring difference of 38. We also show the difference images between 39 Parameter Value Ring Diameter (mm) 262 Detectors per ring 504 Number of rings 48 Rays (LORs) per angle 288 FOV (radial) (mm) 190 Detector crystal size (mm 3 ) 1.5×1.5×10 Table 3.1: MicroPET Focus 220 Scanner Specifications sinograms obtained from the 3D geometric projector and IRB projectors. All three projectors (Cases A,B and C in figure 3.3) produce very small errors for the small FOV (left column) but the error increases for the larger FOV (right column), particularly towards the edge of FOV . Case A shows the most uniform error around the center of the sinogram and yields the smallest normalized root-mean-squared (RMS) error defined as RMSE = ky r −y p k ky r k ×100(%) (3.7) wherek·k is the Euclidean norm, andy r andy p are the sinograms computed usingP geom and the IRB projectors, respectively. For the larger FOV , Case B shows the largest error, and Case C, which uses the line integral based 2D projector, shows a more structured difference image compared to Case A, although the normalized RMS values differ only slightly. Run times for geometric forward and back projection are shown in Table 3.2 for the 3D geometric and IRB projectors. Note that these are the run times for geometric projection only, and exclude the other factors in the factored matrix in (3.1). We used an AMD Opteron 870 workstation with four 2.0GHz dual core CPUs and 32Gbytes RAM. The projection times are measured using all eight processors by multithreading. We obtained more than five times speed up using multithreading compared to their single threaded versions. 40 50 100 150 200 -5 0 5 10 -5 0 5 10 -5 0 5 10 (a) (c) (d) (b) RMSE = 0.77% RMSE = 2.36% RMSE = 0.78% RMSE = 2.52% RMSE = 0.95% RMSE = 2.40% Figure 3.6: Sinograms for a ring difference of 38 computed from the digital 3D Hoffman brain phantom using two different transaxial voxel sizes: (left) 256x256x95 .4mm voxel size, (right) 256x256x95 0.8mm voxel size. (a) The reference sinograms from the 3D geometric projector; (b) difference images between reference and IRB projector Case A in figure 3.3; (c) difference images for Case B; (d) difference images for Case C. The normalized RMS error in (3.7) is indicated below each difference image. 41 Three different images were used for the comparison: Image I of size128×128×95 with 0.4mm transaxial voxel size, Image II of size256×256×95 with a transaxial voxel size of 0.4mm, and Image III of size256×256×95 with a transaxial voxel size of 0.8mm. We show the run times for the inverse rebinning projectors in two parts: (i) the cost for inverse rebinning (IRB) (including corrections) and (ii) the cost of 2D projection for computing the 2D direct sinograms. We do not include the run time for Case C since its time is approximately the same as that for Case A, differing only in the details of the 2D projection. In determining the run time, we truncated the sinograms to the maximum extent of the reconstructed field of view. We then determined the cost of inverse rebinning, which is a function only of the sinogram size but not the image size, and of 2D geometric projection, which depends on both the 2D sinogram size and the size of the image vol- ume and number of voxels. For this reason we see that the biggest savings are achieved when images are reconstructed on a larger transaxial array of 256 rather than 128 voxels. Serendipitously, we therefore achieve the largest fractional speed up in those cases that are currently most time consuming. In comparing Cases A and B we see that the cost of performing all corrections is contained entirely in the inverse rebinning step and does not affect the 2D projection step and increases the overall cost by 50% to 75%. 3.4.2 Application to MAP Reconstruction: Resolution Properties We evaluated the performance of the inverse rebinning projectors in MAP reconstruc- tion of simulated point sources, digital brain phantom images, and in vivo FDG monkey brain data. Images were reconstructed by maximizing (2.2) using a preconditioned con- jugate gradient (PCG) method with a quadratic penalty function and factored system model as described in [108]. Spatially variant weighting of the quadratic penalty terms by the estimated diagonal elements of the Fisher information matrix was used to achieve 42 Projector Image I Image II Image III 3D geometric 14.33s 99.93s 173.15s IRB Run time 4.16s 8.15s 12.75s Forward projector: (3.33 + 0.83) (4.45 + 3.70) (7.65+5.10) projection Case A Speed-up 3.44 12.26 13.58 IRB Run time 3.18s 6.75s 8.53s projector: (2.35 + 0.83) (3.05 + 3.70) (3.43+5.10) Case B Speed-up 4.51 14.80 20.30 3D geometric 16.50s 114.15s 202.77s IRB Run time 4.68s 9.60s 14.31s Back projector: (3.65+ 1.03) (5.05 + 4.55) (8.04 + 6.27) projection Case A Speed-up 3.53 11.89 14.07 IRB Run time 3.58s 7.75s 10.56s projector: (2.55 + 1.03) (3.20 + 4.55) (4.29+6.27) Case B Speed-up 4.61 14.73 19.20 Table 3.2: Comparison of computation times for forward and back projections with different projectors. For the IRB projectors, the run times for inverse rebinning and 2D projection are also shown separately. See text for details. approximate count-independence in reconstructed image resolution [43, 67, 106]. The digital brain phantom and in vivo monkey scans were precorrected for randoms and we therefore used the shifted-Poisson modification [108, 136] of the Poisson likelihood in (2.1) for these data. We compared the performance of the new inverse rebinning pro- jectors with the 3D geometric projector P geom , in each case replacing P geom in (3.1) with the inverse rebinning approximation. In these studies we do not include positron range modeling [110] in the factorization in (3.1). The detector response modelP det.blur used a nonstationary 2D blur kernel applied to each sinogram. The blur kernels were computed based on Monte Carlo simulation as described in [108] and account for block- detector effects as well as intercrystal scatter and penetration and photon-pair noncolin- earity [113]. All reconstructions used an initial estimate computed using 2 iterations of 3D OSEM with 6 subsets. To evaluate the resolution properties of the IRB projectors, we simulated point sources in a uniform background (size: 256×256×95, transaxial voxel size: 0.8mm), 43 with single voxel point sources at an intensity ratio of 40:1 relative to background. The non-zero background was used to avoid artificially small full-width-half-maximum (FWHM) resolution as a result of the nonnegativity constraint. Since we were inter- ested in comparing resolution properties, noise was not added to this data. It has been established elsewhere [42, 106] that reconstructions of noiseless data provide a good approximation to the mean of reconstructions from Poisson data but avoid the need for large scale Monte Carlo studies. In each case shown, we used 100 iterations of MAP reconstruction with a smoothing parameter value ofβ = 0.001. Figure 3.7 and 3.8 show radial and tangential resolution in reconstructed point sources for all three IRB projectors and the 3D geometric projector. The resolutions were measured and plotted in three different planes (9-th, 25-th and 49-th plane from the top of the image volume) as a function of radial distance from the image center. Since the data were noiseless, one might expect the case where the projector model exactly match that used to simulate the data, i.e. the 3D geometric P geom , to produce near perfect resolution. However, in figure 3.7 we see a clear fall off in radial reso- lution with distance from the center of the field of view. We found that running more iterations did result in further improvements in resolution, indicating that PCG had not fully converged after 100 iterations. However the same overall trend remained, ie. the resolution decreased away from the center of the field of view. Furthermore, since a very small value of β was used, instabilities occur at higher iterations when using the inverse rebinning operators because of illconditioning coupled with their mismatch with theP geom operator. In practice, such a small value ofβ would rarely be used. Increasing β would result both in faster convergence and in more stable solutions. It would also lead to reduced differences between the projectors. Consequently the results presented here could be interpreted as a worse case analysis of the differences between projectors. 44 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (a) 9-th plane 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (b) 25-th plane 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (c) 49-th plane Figure 3.7: Radial resolution (FWHM) in three different planes for the 3D geometric and inverse rebinning projectors for 100 iterations of MAP andβ = 0.001. 45 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (a) 9-th plane 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (b) 25-th plane 0 10 20 30 40 50 60 70 80 90 100 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (c) 49-th plane Figure 3.8: Tangential resolution (FWHM) at three different planes for the 3D geometric and inverse rebinning projectors for 100 iterations of MAP andβ = 0.001. 46 The curves in figure 3.7 show similar drop off in radial resolution with distance from the center for all four projectors. To confirm that the loss in radial resolution is a result of our modeling of detector response, we repeated the analysis without including P det.blur in generation of the sinograms or their reconstructions. We found in this case that resolution was approximately shift invariant. In figure 3.8 tangential resolution is approximately uniform with distance from the center and there are clearer differences between the projectors, with P geom performing best. The projector using the 2D line integral projector (Case C) performs worst. Rather surprisingly, Cases A and B show similar performance indicating that including arc correction and axial resampling has little impact on either tangential or radial resolution in the transaxial plane. In fact, performance drops slightly when including these corrections although these differences between Cases A and B are probably insignificant in practice. Figure 3.9 shows axial resolution for lines parallel to the central axis at three differ- ent displacements plotted as a function of the axial distance from the central transaxial plane. Differences between the projectors increase as we move away from the central axis with highest resolution resulting from the 3D geometric projector. In this case, the corrections included in Case A have a clear impact, resulting in noticeably better resolu- tion than when they are omitted as in Case B. The axial resolution for Case B improves as we move towards the boundary as figure 3.9 shows. This possibly counter intuitive observation may be due to the fact that the central planes include a higher fraction of sinograms with large oblique angles than do the outer planes. Since the approximations implicit in Case B become less accurate with increasing oblique angle, the resulting errors are worse towards the axial center of field of view. Conversely, since oblique angle effects are modeled more accurately in Case A, there is no significant variation in axial resolution in the axial direction. Overall, the 2D line integral projector of Case C performs worst in axial resolution. 47 0 5 10 15 20 25 30 35 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (a) 60mm from the center axis 0 5 10 15 20 25 30 35 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (b) 72mm from the center axis 0 5 10 15 20 25 30 35 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Distance from the center (mm) FWHM (mm) 3D Geometric Projector IRB Projector (Case A) IRB Projector (Case B) IRB Projector (Case C) (c) 82mm from the center axis Figure 3.9: Axial resolution (FWHM) for three lines parallel to the central axis of the scanner for the 3D geometric and inverse rebinning projectors for 100 iterations of MAP andβ = 0.001. 48 We also investigated the effect of arc correction without also including axial resam- pling for correcting oblique angles, as is common practice when using Fourier rebinning. We found that this actually resulted in worse performance, in terms of resolution, than Case B where no corrections are included. Consequently, for this application, it is better not to include any corrections than to include only arc correction. To summarize these resolution results, we see that Case A outperforms Case C in tangential and axial resolution. Since they differ only in the form of their 2D projector (solid angle based vs. line integral) and the computation costs are similar, we dropped Case C from further consideration. In comparing Cases A and B, which differ in whether or not they include arc and oblique angle corrections, we see little difference in the transaxial plane but substantial differences in axial resolution, particularly as we move towards the transaxial edge of the field of view. This suggests that Case B may be adequate for smaller objects (< 60mm radius for this scanner) while Case A would result in improved performance for objects taking up a larger fraction of the field of view. 3.4.3 Application to MAP Reconstruction: Phantom and in vivo Studies To investigate the performance of IRB projectors in a structured phantom, we generated a 3D sinogram (MRD=47, span=3) from a 256× 256× 95 digital 3D Hoffman brain phantom image with the transaxial voxel size 0.8mm so that the phantom filled the full field of view. Poisson noise was added to simulate 2.0× 10 9 total counts with a uniform 10% randoms background. A separate randoms process was then simulated with the same mean and subtracted to simulate the process of randoms precorrection. We reconstructed images using 30 iterations of PCG with a hyperparameter β = 0.1, following 2 iterations of 3D OSEM with 6 subsets. figure 3.10 shows a transaxial section 49 (a) 3D Geometric (b) IRB projector (Case A) (c) True (d) Difference 50 100 150 200 250 0 5 10 15 20 25 True 3D Geometric Projector IRB Projector (Case A) (e) Trans-axial Profile Figure 3.10: Reconstructed transaxial sections through the brain phantom (c) using the 3D geometric (a) and inverse rebinning (b) projectors; (d) is the difference between (a) and (b); (e) shows a profile along the white dashed line in (c) 50 and profile of the images reconstructed using the 3D geometric projector and the IRB projector in Case A. The difference image confirms the observation in the resolution study that the images are very similar near the center and differ most towards the edge of the field of view. Finally we compare reconstructions for an in vivo study using a 5.5kg cynomolo- gus macaque injected with 6.9 mCi FDG. The data were collected from the microPET Focus 220 scanner. A total of 8.2× 10 8 counts were detected. The image was recon- structed using a 256x256x95 array with transaxial voxel size of 0.8mm. In figure 3.11 we show sagittal, coronal and transaxial sections through reconstructions of the head obtained using MAP reconstruction with 30 iterations of PCG with a hyperparameter β = 0.1, following 2 iterations of 3D OSEM with 6 subsets. We compare images for Cases A and B to that using the 3D geometric projector. For comparison we also include a reconstruction using Fourier Rebinning (FORE) plus 4 iterations of 2D OSEM. This version of 2D OSEM used an on-the-fly, rather than stored projector/backprojector pair, calculated using line integrals with linear interpolation. Finally, we show a reconstruc- tion using the analytic 3D reprojection (3DRP) method [65] with a ramp filter. Throughout the field of view, the MAP images, figure 3.11a,b and c, are more similar than the alternative FORE+2D OSEM and 3DRP reconstructions. Since the snout of the monkey is close to the edge of field of view of the scanner, as can be seen in the sagittal view, it appears that the IRB projectors can achieve very similar performance to the 3D geometric projector in our earlier work [108, 109]. 3.4.4 Application to MAP Reconstruction: Reconstruction Time Comparison To compare the total reconstruction times, MAP reconstruction code was run in different three image sizes as done in Table 3.2 in the same computer. We reconstructed images 51 (a) (b) (c) (d) (e) Full FOV Figure 3.11: Sagittal (left), transaxial (middle) and coronal (right) views of recon- structed images on the microPET Focus 220 of an FDG study of a macaque monkey using (a) MAP with the 3D geometric projector; (b) MAP with IRB projector, Case A; (c) MAP with IRB projector, Case B; (d) FORE + 2D OSEM; (e) 3DRP with a ramp filter. 52 Image I Image II Image III Reference 27.7 I 125.6 240.2 IRB projector: Time 13.9 22.2 35.4 Case A Speed-up 1.99 5.66 6.79 IRB projector: Time 12.4 20.8 29.8 Case B Speed-up 2.23 6.04 8.06 Table 3.3: Comparison of the MAP reconstruction time (Unit: Minutes). using 20 iterations of PCG with a hyperparameterβ = 0.1, following 2 iterations of 3D OSEM with 6 subsets. The reconstruction code is also programmed in the way to utilize all processors by multithreading. Table 3.3 shows the comparison result. ’Reference’ indicates the our original MAP code which has the original 3D geometric projector, and others show the MAP codes using fast projectors, Case A and B. As the table shows, the speed up factor increases as the image FOV does as expected from the result in Table 3.2, and 2 to 8 times speed up was achieved. 3.5 Conclusions We have developed an approach for combining inverse Fourier rebinning with iterative 3D PET reconstruction as a means of reducing reconstruction time while retaining the advantages of fully 3D iterative reconstruction. For accurate use of inverse rebinning in a practical detector-ring based PET scanner, we applied arc correction and axial angu- lar resampling and showed that these corrections improve sinogram and reconstructed image quality, although the effects are most pronounced towards the edge of the field of view. With the full corrections we achieve an order of magnitude speed up using IRB compared to 3D geometric projectors. When implemented as part of a MAP reconstruc- tion algorithm, small reductions in resolution were observed, but in vivo results indicate 53 small differences as a function of the projector when these are compared with alter- native reconstructions using 3DRP and FORE + 2D OSEM. This small image quality difference may indicate that considering the spatially varying detector response model is more important than having the depth dependant geometric sensitivity model in the 3D geometric projector. 54 Chapter 4 Unified Framework for Fourier Rebinning 4.1 Introduction Fourier rebinning concept plays a useful role not only for data reduction, but also for developing a fast projector, in fully 3D PET image reconstructions. Here, we extend the similar concept to 3D TOF PET reconstructions. We derive a generalized projection slice theorem which provides the relationship between the 3D object function and the 3D TOF PET data. The theorem is general in the sense that it also applies to other data sets such as 3D non TOF, 2D TOF and 2D non TOF PET data. Based on the theorem, we develop a unified framework which is used for deriving all the mappings between different data sets as shown in Figure 4.1. We first derive the generalized projection slice theorem based on an analytical data model of 3D TOF PET data. The time of flight information is included in the line integral based data model as a convolution kernel function through lines of response. Using the analytical model, we derive the generalized projection slice theorem in the frequency domain. The projection slice theorem gives the relationship between a common object function and all the different data sets, so we can find mappings between different data sets in Figure 4.1 by using the Fourier transform of the common object function as an intermediate quantity which links two different data sets. The mapping equations are 55 A B C D E 3D TOF Data 3D non TOF Data 2D non TOF Data 2D TOF Data Figure 4.1: All possible mappings between 3D TOF PET data and other lower dimen- sional data sets all exact, and can be utilized for designing fast projectors or data rebinnings for fast 3D TOF PET data reconstructions. 4.2 Analytical Data Model and Generalized Projection Slice Theorem Here, we extend the conventional projection slice theorem to the one applied to 3D TOF PET data. The derivation is started from a line integral based analytical data model. Based on the model, we derive projection slice theorem for 3D TOF PET data which shows the relationship between the object function and 3D TOF PET data in the fre- quency domain. The theorem is general in the sense that it also includes the previous non-TOF PET case. 56 4.2.1 Line Integral Model with Spatially Invariant TOF Kernel Function Equation (2.5) in Chapter 2 is the analytical model for 3D non TOF PET data. Here, we add a TOF kernel function to this line integral model to represent 3D TOF PET data. Following the notations in [28], the 3D TOF PET datap for a cylindrical scanner geometry, illustrated in Figure 4.2, can be written as a line integral weighted by the TOF kernelh: p(s,φ,z,δ;t) = √ 1+δ 2 Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z +lδ)h(t−l √ 1+δ 2 )dl (4.1) wheref represents the 3D object,s andφ are the radial and angular coordinates, respec- tively,z is the axial midpoint of each LOR,δ is the tangent of the oblique angleθ, andt is the TOF variable. Notice that the scaling term √ 1+δ 2 is also included here compared to (2.5). The TOF kernel h, often modeled as a Gaussian function [19], is assumed to be spatially invariant, so that the integral in (4.1) is a convolution of the object function f with the kernel function through an LOR. As a special case, we can write the stacked 2D TOF direct sinograms by settingδ, the tangent of the oblique angle, to zero in (4.1) as follows: p(s,φ,z,0;t) = = Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z)h(t−l)dl. (4.2) 57 x y θ φ t (a) (b) δ = tanθ z LOR s R TOF kernel h Figure 4.2: (a) Top and (b) side view of a cylindrical 3D PET scanner. For each line of response (LOR), the object is multiplied by the TOF kernel h and integrated along the line to form the TOF data. The absence of time of flight information is equivalent to settingh = 1 in which casep in (4.1) reduces to the 3D non-TOF data: p nonTOF (s,φ,z,δ) = √ 1+δ 2 Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z+lδ)dl, (4.3) which is equivalent to (2.5) except the constant scaling term √ 1+δ 2 . The same thing applies to (4.2) to obtain 2D non-TOF data model. Therefore, the model in (4.1) includes all data sets, 3D TOF, 3D non TOF, 2D TOF and 2D non-TOF PET data. 4.2.2 Taking Fourier Transform in TOF Bin Direction The key difference of our approach in handling the 3D TOF PET data from others [28, 133] is taking Fourier transform of the 3D TOF PET datap in (4.1) with respect to the 58 TOF variablet. Under the assumption thath is spatially invariant, the Fourier transform of (4.1) in the TOF variablet becomes: P(s,φ,z,δ;ω t ) = Z ∞ −∞ p(s,φ,z,δ;t)e −jωtt dt = √ 1+δ 2 H(ω t ) Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z +lδ)e −jωtl √ 1+δ 2 dl (4.4) whereω t is the frequency variable corresponding tot andH is the 1D FT ofh. A stacked 2D version of the above equation, that is, the 1D Fourier Transform ofp(s,φ,z,0;t) in (4.2) with respect tot, is given by P(s,φ,z,0;ω t ) =H(ω t ) Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z)e −jωtl dl. (4.5) Taking the Fourier transform with respect to the TOF variable may cause aliasing since the sampling frequency in the TOF variable is low. However, the bandwidth along the LOR is limited through convolution with the TOF kernel which limits the effects of aliasing. The Fourier transform with respect to the TOF variablet takes the TOF kernel term H out of the integrals in (4.4) and (4.5) so P 3D in (4.4) is similar to the non-TOF p nonTOF in (4.3) except for the scaling factor H(ω t ) and the phase modulation term exp(−jω t l √ 1+δ 2 ). This similarity between the non-TOF datap nonTOF and the Fourier transformP of the TOF datap, with respect tot, enables us to extend the properties of non-TOF data to the TOF case. The generalized projection slice theorem is the first one which is extended from the non TOF PET data case to the TOF PET data case. The property, so called bow-tie property [112], for the Radon transform (equivalently, 2D 59 non TOF sinogram based on the analytical model) can be extended to 2D TOF PET data. The bow-tie property and more detail description of the 2D version of the projec- tion slice theorem is described in the following. 4.2.3 Analytical Properties of 2D TOF PET Data Two properties of 2D TOF PET data will be described before explaining the general- ized projection slice theorem. These properties can be also extended to 3D TOF PET data. All derivations are started from (4.5). However, we will omit the axial variable z in the equation since we are interested in just one 2D TOF PET sinogram plane, so 2D TOF sinogram will be parametrized as p(s,φ;t) and the non TOF 2D sinogram as p nonTOF (s,φ) in here. Projection Slice Theorem By omittingz in (4.5), the equation for a 2D TOF PET sinogram is expressed as follows: P(s,φ;ω t ) =H(ω t ) Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ)e −jωtl dl. (4.6) Now, we want to take further Fourier transform ins direction as follows: ˜ ℘(ω s ,φ;ω t ) =H(ω t ) Z ∞ −∞ Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ)e −j(ωtl+ωss) dlds =H(ω t ) Z ∞ −∞ Z ∞ −∞ f(x,y)e −i((ωscosφ−ωtsinφ)x+(ωssinφ+ωtcosφ)y) dxdy =H(ω t )F(ω s cosφ−ω t sinφ,ω s sinφ+ω t cosφ) (4.7) whereF is the 2D Fourier transform of 2D objection functionf(x,y). Compared to the conventional case of the 2D Radon transform, we have the extra coordinate information 60 x y p(s,φ;t) t y x t s φ s f(x,y) φ LOR (a) Space domain ω x ω y ω s F(ω x ,ω y ) φ ω t H(ω t ) ˜ ℘(ω s ,φ;ω t ) (b) Frequency domain Figure 4.3: (a) The variables s and t constitute an orthogonal coordinate rotated by φ fromx-y coordinate. (b) Therefore, by taking Fourier transform with respect tos andt, their frequency coordinate is also rotated by φ. The strip weighted by H in (b) shows the extra information produced by having TOF information. When TOF information is not available, we have information only in the line whereω t = 0. t andω t as shown in Figure 4.3. In the figure,t ands variables constitute an orthogonal coordinate rotated by φ from the x-y coordinate in the space domain. Therefore, the same rotation occurs in the frequency domain with respect to the coordinates constituted by their frequency variables as illustrated in Figure 4.3(b). The conventional non-TOF case is when ω t is zero, that is, the sum in the TOF bin direction, then it tells that ˜ ℘(ω s ,φ;0) for fixedφ corresponds to the solid line in Figure 4.3(b) indicated asω s axis in the sinogram domain. However, in the TOF case, we have the extra TOF information which extends the line to a strip weighted by H, the Fourier transform of TOF kernel functionh, through TOF bin direction. In principle, when we have a complete data set for all φ from 0 to 180 ◦ , the increased information due to having the TOF information would be redundant, but will help to increase SNR in real situations. The projection slice theorem in (4.7) is used in the next chapter for implementing a direct Fourier reconstruction of 2D TOF PET data. 61 x φ y R M s p nonTOF (s,φ) f(x,y) LOR (a) Space domain ω φ ω s R M 1 (b) Frequency domain Figure 4.4: The frequency response of the conventional 2D PET sinogram (Radon trans- form) has bow-tie shaped frequency response where the bow-tie boundary is determined by the radius of non zero region of object functionf. Bow-tie Property of 2D TOF PET Sinogram Rattey and Lindgren showed that the 2D Fourier transform of a 2D non-TOF PET sino- gram has almost all of its energy confined to a characteristic bow-tie shaped region, where the size of the bow-tie is determined by the radius of the object field-of-view [112] as illustrated in Figure 4.4. This is called bow-tie property. The property is derived by analyzing the impulse response of the 2D Radon transform data in the frequency domain. Here, we apply the similar technique to examine the frequency response of a 2D TOF PET sinogram, and show that the 2D TOF PET data also have the similar property. To derive the property, consider P(s,φ;ω t ) in (4.6), the 2D TOF sinogram after taking the Fourier transform int. First, we take the 1D Fourier transform ofP(s,φ;ω t ) 62 with respect to s for the case. And a 2D impulse at (x 0 ,y 0 ) is used for the object as f(x,y) =δ 2 (x−x 0 ,y−y 0 ), to see the impulse response: ˜ ℘(ω s ,φ;ω t ) = Z ∞ −∞ P(s,φ;ω t )e −iωss ds =H(ω t ) Z ∞ −∞ Z ∞ −∞ f(x,y)e −i(ωt(ycosφ−xsinφ)+ωs(xcosφ+ysinφ)) dxdy =H(ω t ) Z ∞ −∞ Z ∞ −∞ δ 2 (x−x 0 ,y−y 0 )e −i(ωt(ycosφ−xsinφ)+ωs(xcosφ+ysinφ)) dxdy =H(ω t )e −ir 0 √ ω 2 t +ω 2 s cos(φ−φ 0 +φ 1 ) (4.8) wherer 0 = p x 2 0 +y 2 0 ,φ 0 = arctan(y 0 /x 0 ),φ 1 = arctan(ω t /ω s ). We make the change of variablescosφ−lsinφ = x andssinφ+lcosφ = y in (4.8). Next we take the 1D Fourier transform of℘(ω s ,φ;ω t ) with respect toφ: P(ω s ,ω φ ;ω t ) = Z ∞ −∞ ˜ ℘(ω s ,φ;ω t )e −iω φ φ dφ (4.9) = H(ω t ) Z ∞ −∞ e −ir 0 √ ω 2 s +ω 2 t cos(φ−φ 0 +φ 1 ) e −iω φ φ dφ. (4.10) We note that, since ˜ ℘ is a2π periodic function inφ, it can be directly represented through its Fourier series. However, here we follow the notation of [112] and use the general- ized Fourier transform of the periodic function in (4.10), which results in the weighted impulse train inω φ : P(ω s ,ω φ ;ω t ) =H(ω t ) ∞ X m=−∞ 2πe −im(π/2+φ ′ 0 ) J m r 0 q ω 2 s +ω 2 t δ(ω φ −m) (4.11) whereφ ′ 0 =φ 0 −φ 1 andJ m is themth order Bessel function of the first kind: J m (r) = 1 2π Z π −π e i(mφ−rsinφ) dφ. 63 Figure 4.5: Example of the magnitude response of the series in (4.11). The figures shows the modified bow-tie region bounded by hyperbolic lines where the outside of the bow-tie region is nearly zero. As observed by Rattey and Lindgren (1981) for the non TOF case (equivalent toω t = 0 in (4.11)), the properties of the Bessel function are such that P(ω s ,ω φ ;ω t ) in (4.11) is approximately zero for|ω φ | > |r 0 p ω 2 s +ω 2 t | + 1 and the energy of the function is approximately contained within the region B(ω t ) = (ω s ,ω φ ) :|ω φ |≤ r 0 q ω 2 s +ω 2 t +1 (4.12) This region is bounded by a hyperbolic curve in the (ω s ,ω φ ) plane forω t > 0 as shown in Figure 4.5. When ω t = 0,B becomes the well known bow-tie region for non-TOF data [112]. 64 (a) center TOF bin (b) 2nd TOF bin from the center bin (c) After 1D Fourier transform in t at ω t frequency index zero (d) After 1D Fourier transform in t at ω t frequency index 2 Figure 4.6: 2D TOF sinograms and their magnitudes after taking 1D Fourier transform in t; The horizontal axis is the sinogram angle φ ∈ [0,2π] and the vertical axis is the radial variables. Validation of the analysis was conducted using simulations based on the parameters listed in Table 5.1. Figure 4.6 shows examples of TOF data for a direct plane (δ = 0). Figures 4.6(a) and (b) show the raw TOF data p 3D in (4.1) for different values of the TOF variablet. Figures 4.6(c) and (d) show the magnitude of the TOF data after taking the 1-D Fourier transform int, as defined in (4.4), for different values ofω t . Note that P forω t = 0, shown in Figure 4.6(c), is equivalent to the non-TOF data. Figure 4.7 shows examples of the magnitude of the Fourier transformP(ω s ,ω φ ;ω t ) of P(s,φ;ω t ) for a direct plane (δ = 0) at 250 ps TOF resolution. Each block in the figure shows the magnitude for a different frequency index in ω t . As ω t increases, the separation between the two branches of the hyperbola increases. The solid hyperbolic curve in Figure 4.7(d) shows the boundaries of the bow-tie region calculated from the relationship in (4.12). Note that the calculated values ofP(ω s ,ω φ ;ω t ) are indeed largely contained in these hyperbolically confined regions. 65 (a) 2nd frequency (DFT) index (b) 3rd frequency index (c) 4th frequency index (d) 5th frequency index Figure 4.7: 2D Fourier transformP(ω s ,ω φ ;ω t ) ofP(s,φ;ω t ) in (4.11) for four different ω t frequency indices. The horizontal axis represents ω φ and the vertical axis ω s . Also shown in (d) is the hyperbolic curve bounding the effective support of the bow-tie for TOF data in the 5th frequency bin. This property could be used for reduction of noise or scatter components in the data, or for sinogram restoration for limited angle tomography [62]. Here, we will show a simple way of noise reduction filtering to the 2D TOF PET noisy sinogram data. The filter is a modified bow-tie shape filter in the frequency domain, where each filter is designed differently for each TOF frequency index. Figure 4.8 is for showing the bow tie filters and the filtered 2D TOF sinograms of a noisy 2D TOF sinogram. To the 2D TOF sinogram, first we take FT through TOF bin direction to get the signal in (4.6) as shown in Figure 4.8(a) for its magnitude response at a fixedω t (3rd DFT index from zero). Figure 4.8(b) shows the magnitude plot of 2D FT which isP(ω s ,ω φ ;ω t ) in (4.10). Figure 4.8(c) shows the bow-tile filter in the 2D FT domain where the center is the origin of the radialω s (vertical axis) and angularω φ (horizontal axis) frequency axes. Notice that we also cascaded a low-pass filter in the radial direction to the bow-tie filter, and here we call the filting process by two cascased filters the bow-tie filtering. The low pass filter in the radial direction is shown in Figure 4.8(d). To see the bow-tie filtering 66 (a) (b) (c) (d) (e) (f) Figure 4.8: Bow tile filtering to noisy 2D TOF sinograms in (4.5); (a) magnitude response of 2D TOF sinogram after FT to TOF bin direction at 3-th TOF frequency index from zero frequency, (b) magnitude plot of 2D FT of (a), (c) Bow-tie filter in 2D frequency domain, (d) Radial filter in 2D frequency domain, (e) magnitude plot of Bow-tie filtered sinogram of (a), (f) magnitude plot of radial filtered sinogram of (a) effect, we applied two filters shown in Figure 4.8(c) and (d) to the frequency response in Figure 4.8(b). Figure 4.8(e) and (f) are the magnitude plots after the filterings using two filters in Figure 4.8(c) and (d), respectively. As you can see the bow-tie filtering well suppresses noise. 67 4.2.4 Generalized Projection Slice Theorem We derived the projection slice theorem for 2D TOF PET data in (4.7). Now, we extend the theorem to 3D TOF PET data which will be called a generalized projection slice theorem. The theorem is derived based on the analytical data model in (4.4). We take the 2D Fourier transform ofP(s,φ,z,δ;ω t ) in (4.4) with respect tos andz as follows: ℘(ω s ,φ,ω z ,δ;ω t ) = Z ∞ −∞ Z ∞ −∞ P 3D (s,φ,z,δ;ω t )e −i(ωss+ωzz) dsdz = √ 1+δ 2 H(ω t )˜ ℘(ω s ,φ,ω z ,δ;ω t ) (4.13) where ˜ ℘(ω s ,φ,ω z ,δ;ω t ) = Z ∞ −∞ Z ∞ −∞ Z ∞ −∞ f(scosφ−lsinφ,ssinφ+lcosφ,z +lδ) ·e −iωtl √ 1+δ 2 dle −i(ωss+ωzz) dsdz (4.14) = Z ∞ −∞ Z ∞ −∞ F z (x,y,ω z )e −ix(ωscosφ−(ωt √ 1+δ 2 −δωz)sinφ) ·e −iy(ωssinφ+(ωt √ 1+δ 2 −δωz)cosφ) dxdy (4.15) = F 3D (˜ ω x , ˜ ω y ,ω z ). (4.16) where ˜ ω x =ω s cosφ−(ω t √ 1+δ 2 −δω z )sinφ, ˜ ω y =ω s sinφ+(ω t √ 1+δ 2 −δω z )cosφ, F 3D is the 3D Fourier transform off, andF z (x,y,ω z ) in (4.15) is the 1D Fourier trans- form off(x,y,z) with respect to z. The change of variable, scosφ−lsinφ = x and ssinφ +lcosφ = y, and the shifting property of Fourier transform are used to obtain (4.15). One can view (4.16) as an extension to 3D TOF PET data, which is equivalent to (3) in [133], although here the result is derived directly in terms of the parametrization 68 used for data collected in a cylindrical scanner. Then, the final form of the generalized projection slice theorem for 3D TOF PET data can be written as follows: ℘(ω s ,φ,ω z ,δ;ω t ) = √ 1+δ 2 H(ω t )F 3D (ω s cosφ−χsinφ,ω s sinφ+χcosφ,ω z ) (4.17) where χ =ω t √ 1+δ 2 −δω z . (4.18) The theorem Equation (4.17) is general in the sense that it also applies to other cases: 3D non TOF PET case (when ω t = 0), 2D TOF case (when δ = 0) and 2D non TOF case (when ω t = 0 and δ = 0). Note that ω t = 0 corresponds to the DC component in the TOF variable direction, which is the sum through the TOF bin variable t and therefore represents non TOF data. This theorem can be used as the basis for developing a direct Fourier reconstruction method of 3D TOF data as done for 3D non TOF PET data reconstruction in [83]. Notice that the last argument of the functionF 3D is onlyω z in (4.17). Therefore, after taking Fourier transform in thez direction, the theorem, for each fixedω z , becomes very similar to its 2D version in (4.7). There are differences only in scaling term √ 1+δ 2 and χ in the arguments in (4.17). So we can extend the bow-tie property in (4.12) to this 3D TOF PET data after taking Fourier transform inz. The 3D TOF direct Fourier reconstruction method can be also implemented slice by slice for each ω z . In order to take the Fourier transform in thez direction, the truncated part of the 3D TOF oblique sinograms needs to be estimated first. This missing data problem will be also addressed in the next chapter when we implement the exact rebinning where the Fourier transform inz direction is required. 69 4.3 All Exact Mappings Between Different PET Data Sets Once we have obtained the generalized projection slice theorem, the problem finding all mappings between different data sets is straightforward. The projection slice theorem provides us with the relationship between a common object functionf and all different data sets in the frequency domain, so the Fourier transform of object functionf is used as an intermediate medium to derive a mapping equation between two different data sets. Now, we derive all the mappings in Figure 4.1 one by one. 4.3.1 Mapping D: 3D non TOF and 2D non TOF PET Data This mapping is the already known mapping derived previously in [32, 70]. Here we derive in a different way. From (4.17), the 3D non TOF PET data can be written by settingω t to zero as follows: ℘(ω s ,φ,ω z ,δ;0) = √ 1+δ 2 H(0)F 3D (ω s cosφ+δω z sinφ,ω s sinφ−δω z cosφ,ω z ). (4.19) This is the projection slice theorem for 3D non TOF PET data, and was introduced in [70] when the exact mapping equation is derived. The 2D non TOF PET data, stacked 2D non TOF sinograms, can be expressed as follows by settingω t andδ to zero: ℘(ω s ,φ,ω z ,0;0) =H(0)F 3D (ω s cosφ,ω s sinφ,ω z ). (4.20) Then, the following mapping equation, ℘(ω s ,φ,ω z ,δ;0) = √ 1+δ 2 ℘(ω ′ s ,φ ′ ,ω ′ z ,0;0), (4.21) 70 can be obtained by equating both arguments of F 3D function in (4.19) and (4.20) as follows: ω ′ s cosφ ′ = ω s cosφ+δω z sinφ (4.22) ω ′ s sinφ ′ = ω s sinφ−δω z cosφ (4.23) ω ′ z = ω z . (4.24) That is, finding the exact mapping equation of (4.21) is equivalent to solve the equations (4.22)-(4.24). The equations can be solved by rewriting (4.22) and (4.23) as follows: ω ′ s cosφ ′ = q ω 2 s +(δω z ) 2 cos(φ+β) =ω s s 1+ δω z ω s 2 cos(φ+β) (4.25) ω ′ s sinφ ′ = q ω 2 s +(δω z ) 2 sin(φ+β) =ω s s 1+ δω z ω s 2 sin(φ+β) (4.26) where β = arctan(−δω z /ω s ). From the equations, we can obtain the solution, ω ′ s = ω s q 1+(δω z /ω s ) 2 andφ ′ =φ+β. Then, (4.21) can be written as follows: ℘(ω s ,φ,ω z ,δ;0) = √ 1+δ 2 ℘(ω s s 1+ δω z ω s 2 ,φ−arctan δω z ω s ,ω z ,0;0). (4.27) From (4.27), we can obtain the exact equation (2.7) derived in [32] by further taking the Fourier transform inφ direction. Note that the scaling constant √ 1+δ 2 in (4.27) is different from previous (2.7) because the data model in (4.1) includes this scaling term while it was considered later as a scaling constant in (2.5) in [32]. The two equations (4.19) and (4.20) show that once the PET data are Fourier trans- formed in z direction, the relationship between PET data and the Fourier transform of object functionf is made in 2D (ω x ,ω y ) plane where each coordinate indicates the first 71 B φ ′ =φ+β ω ′ s = ω s v u u u u t1+   δω z ω s   2 β = arctan   −δω z ω s   φ −δω z ω x ω y ω s = 0 A C |AB| =|ω s | |CB| =|ω ′ s | 6 ABC =β Figure 4.9: Graphical interpretation of the exact inverse Fourier rebinning mappings from 2D non TOF data to 3D non TOF data. Inverse Fourier rebinning is equivalent to finding the coordinates (ω ′ s ,φ ′ ) of a point (for example, point B) on a trajectory (for example, the line joining the points A and B) mapped from the 3D sinogram. In the approximate inverse rebinning, |CB| and∠ABC are approximated by |AB| and −δω z /ω s , respectively. The inside region of the dotted line circle is not covered by 3D non TOF PET data℘(ω s ,φ,ω z ,δ;0). So there is no rebinning mapping available for the 2D non TOF PET data ℘(ω ′ s ,φ ′ ,ω z ,0;0) in this area from the 3D non TOF PET data ℘(ω s ,φ,ω z ,δ;0). and second argument of F 3D function. Then, the exact mapping (4.27) can be graphi- cally interpreted in the (ω x ,ω y ) plane for fixed ω z as Figure 4.9 shows. In Figure 4.9, the line joining the points A and B is a trajectory (ω x ,ω y ) corresponding to 3D non TOF PET data ℘(ω s ,φ,ω z ,δ;0) parametrized by ω s for fixed φ, ω z and δ, through the relationships ω x = ω s cosφ + δω z sinφ and ω y = ω s sinφ − δω z cosφ, as given in equation (4.19). Then, exact inverse Fourier rebinning finds the coordinates (ω ′ s ,φ ′ ) of (ω x ,ω y ) (for example, point B in Figure 4.9) through the relationships ω x = ω ′ s cosφ ′ andω y = ω ′ s sinφ ′ for 2D non TOF data℘(ω ′ s ,φ ′ ,ω z ,0;0) in (4.20), and vice versa for 72 rebinning from 3D non TOF PET data to 2D non TOF data. Note that the inside region of the the dotted line circle in Figure 4.9 is the area where 3D PET data℘(ω s ,φ,ω z ,δ;0) are not available, so there is no rebinning mapping available for the 2D non TOF PET data ℘(ω ′ s ,φ ′ ,ω z ,0;0) in this area from the 3D non TOF PET data ℘(ω s ,φ,ω z ,δ;0). That is, the rebinning is defined only when|ω s | >|δω z | which is also indicated in (15) in [32]. 4.3.2 Mapping B: 3D TOF and 2D TOF PET Data The 2D TOF PET data, the stack of 2D TOF PET sinograms, can be written as follows by settingδ is zero in (4.17): ℘(ω s ,φ,ω z ,0;ω t ) =H(ω t )F 3D (ω s cosφ−ω t sinφ,ω s sinφ+ω t cosφ,ω z ). (4.28) Then, the same method in Chapter 4.3.1 can be applied to find the following exact mapping equation between 3D TOF PET data, ℘(ω s ,φ,ω z ,δ;ω t ), and 2D TOF PET data,℘(ω s ,φ,ω z ,0;ω t ): ℘(ω s ,φ,ω z ,δ;ω t ) = √ 1+δ 2 ℘(ω ′ s ,φ ′ ,ω ′ z ,0;ω t ). (4.29) It can be obtained by solving the following mapping equations: ω ′ s cosφ ′ −ω ′ t sinφ ′ = ω s cosφ−(ω t √ 1+δ 2 −δω z )sinφ (4.30) ω ′ s sinφ ′ +ω ′ t cosφ ′ = ω s sinφ+(ω t √ 1+δ 2 −δω z )cosφ (4.31) ω ′ z = ω z . (4.32) 73 Those equations can be solved by rewriting (4.30) and (4.31) as follows: q ω ′ s 2 +ω ′ t 2 cos(φ ′ +α ′ ) = p ω 2 s +χ 2 cos(φ+β) (4.33) q ω ′ s 2 +ω ′ t 2 sin(φ ′ +α ′ ) = p ω 2 s +χ 2 sin(φ+β) (4.34) whereχ =ω t √ 1+δ 2 −δω z ,α ′ = arctan(ω ′ t /ω ′ s ),β = arctan(χ/ω s ). To find an exact mapping, we need to find a transformation between{ω ′ s ,φ ′ ,ω ′ z ,ω ′ t } and{ω s ,φ,ω z ,ω t } that satisfies (4.30)-(4.32). Note that there are more than one trans- formations that satisfy these three equations since we have more unknown variables. Here, we propose the following mapping satisfying (4.30)-(4.32) using (4.33) and (4.34): ω ′ s =ω s q 1+(χ 2 −ω 2 t )/ω 2 s , φ ′ =φ+β−α ′ , ω ′ z =ω z , ω ′ t =ω t . (4.35) Here, we fixedω ′ t toω t as forω ′ z to avoid the interpolation across the TOF bin direction. Those mappings in (4.35) give inverse rebinning mapping in (4.29). Figure 4.10 is a modified version of Figure 4.9 for explaining this mapping. Here, point B corresponds to the 3D TOF PET data ℘(ω s ,φ,ω z ,δ;ω t ) in the (ω x ,ω y ) plane. The inverse rebinning is equivalent to finding the 2D TOF data℘(ω ′ s ,φ ′ ,ω z ,0;ω t ) where ω ′ s denotes the line segment from point C to B andφ ′ is the sum ofφ and∠ABC. Note that the mapping forω s is only valid forω 2 s +χ 2 ≥ω 2 t to guarantee that the term inside of the square root be non-negative. This also can be interpreted graphically by modifying the Figure 4.10 with the case where|χ|<|ω t | as done in Figure 4.9. Using the inverse of the transformation in (4.35), we can also obtain the following exact rebinning equation: ℘ 2D (ω ′ s ,φ ′ ,ω ′ z ;ω ′ t ) = 1 √ 1+δ 2 ℘ 3D (˜ ω ′ s ,φ ′ +α ′ −β ′ ,ω ′ z ,δ;ω ′ t ) (4.36) 74 ω x ω y α B ω t α φ φ ′ =φ+α ω ′ s = ω s v u u u u t1+ χ 2 −ω 2 t ω 2 s χ |AB| =|ω s | |CB| =|ω ′ s | A C 6 ABC =α Figure 4.10: Graphical interpretation of the exact inverse Fourier rebinning mappings from 2D TOF to 3D TOF PET data. The exact inverse Fourier rebinning is equivalent to finding the coordinate information of the line segment CB (|CB| and∠ABC) for the line segment AB. where χ ′ =ω ′ t √ 1+δ 2 −δω ′ z , ˜ ω ′ s =ω ′ s q 1+(ω ′ t 2 −χ ′2 )/ω ′ s 2 , β ′ = arctan(χ ′ /˜ ω ′ s ). Similar comment can be made for the validity of this mapping. The mapping is valid only whenω ′ s 2 +ω ′ t 2 ≥χ ′ 2 . All the previously proposed TOF rebinning methods [28, 35, 94] belong to this cat- egory. They all rebin 3D TOF PET data to 2D TOF PET data. Defrise et. al proposed exact rebinning equations in [35] based on a consistency condition expressed by partial 75 differential equations. The equation (4.36) is also exact, and does not require partial derivatives. In the next two subsections, we show the cases where 3D TOF PET data can be rebinned to 2D or 3D non TOF sinograms so that the TOF information is utilized implic- itly so that the rebinned data are not a function of TOF information anymore. 4.3.3 Mapping C: 3D TOF and 2D non TOF PET Data For this mapping, we need to find a mapping satisfying the following equation: ℘(ω s ,φ,ω z ,δ;ω t ) = √ 1+δ 2 {H (ω t )/H(0)}℘(ω ′ s ,φ ′ ,ω ′ z ,0;0). (4.37) Using (4.17) and (4.20), the following mappings need to be satisfied to have the mapping in (4.37): ω ′ s cosφ ′ = ω s cosφ−χsinφ (4.38) ω ′ s sinφ ′ = ω s sinφ+χcosφ (4.39) ω ′ z = ω z , (4.40) where, again, χ = ω t √ 1+δ 2 −δω z . The mapping equation becomes identical to the mapping D if we replacingχ term with−δω z in (4.38)-(4.39). Therefore, we have the following exact mapping by simply changing the−δω z term in (4.26) withχ: ω ′ s =ω s q 1+(χ/ω s ) 2 , φ ′ =φ+β, ω ′ z =ω z , ω ′ t = 0, (4.41) whereβ = arctan(χ/ω s ). 76 Substituting these transformations into (4.37) we obtain the final form of the inverse rebinning relationship: ℘(ω s ,φ,ω z ,δ;ω t ) = √ 1+δ 2 {H(ω t )/H(0)}℘ ω s q 1+(χ/ω s ) 2 ,φ+β,ω z ,0;0 . (4.42) This mapping can be used for rebinning 3D TOF PET data to 2D non TOF PET data, which gives the most data reduction factor. The graphical interpretation in Figure 4.9 can be applied to this mapping by usingχ in place of−δω z . 4.3.4 Mapping A: 3D TOF and 3D non TOF PET Data The mapping A in Figure 4.1 can be expressed as follows: ℘(ω s ,φ,ω z ,δ;ω t ) ={H (ω t )/H(0)}℘(ω ′ s ,φ ′ ,ω ′ z ,δ;0). (4.43) Using (4.17) and (4.19), the following mappings need to be satisfied to have (4.43): ω ′ s cosφ ′ +δω ′ z sinφ ′ = ω s cosφ−χsinφ (4.44) ω ′ s sinφ ′ −δω ′ z cosφ ′ = ω s sinφ+χcosφ (4.45) ω ′ z = ω z , (4.46) whereχ =ω t √ 1+δ 2 −δω z . This equations can be solving by re-writhe (4.44) and (4.45) as follows: q ω ′ s 2 +(−δω z ) 2 cos(φ ′ +α ′ ) = p ω 2 s +χ 2 cos(φ+β) (4.47) q ω ′ s 2 +(−δω z ) 2 sin(φ ′ +α ′ ) = p ω 2 s +χ 2 sin(φ+β) (4.48) 77 whereχ =ω t √ 1+δ 2 −δω z ,α ′ = arctan(−δω z /ω ′ s ),β = arctan(χ/ω s ). By solving above equations forω ′ s ,φ ′ andω z , we obtain the following exact inverse rebinning mapping from 3D non TOF data to 3D TOF data: ω ′ s =ω s s 1+ χ 2 −(δω z ) 2 ω 2 s φ ′ =φ+arctan χ ω s +arctan δω z ω ′ s . (4.49) Note that the mapping E is one special case of this mapping whenδ is zero. Figure 4.10 can be applied to this case by using−δω z in place ofω t . 4.4 Summary We derived all exact mappings between 3D TOF PET data or 3D non TOF PET data and their lower dimensional versions using the generalized projection slice theorem. The theorem was derived based on an analytical data model where the object function is convolved with a spatially invariant TOF kernel function through lines of response. After taking Fourier transform inz direction, all the mappings are derived in 2D(ω x ,ω y ) plane of the function F 3D . That is, all the mapping equations are derived by equating the first and second arguments ofF 3D function in (4.17). So the key mapping equations for each case can be summarized as in Table 4.1. The mapping D is the original Fourier rebinning mapping, and the mapping B was the mapping between 3D TOF PET data and 2D TOF PET data, which can be used for a TOF PET rebinning which was also proposed by others [28, 35, 94]. Other mappings are between TOF PET data and non TOF PET data which can be used for rebinning TOF PET data to non TOF PET data, so that previously developed non TOF PET data reconstruction methods can be used as is 78 Mapping Involved data sets Mapping equations 3D TOF data A ℘(ω s ,φ,ω z ,δ;ω t ) ω x =ω s cosφ−χsinφ =ω ′ s cosφ ′ +δω z sinφ ′ 3D non TOF data ω y =ω s sinφ+χcosφ =ω ′ s sinφ ′ −δω z cosφ ′ ℘(ω ′ s ,φ ′ ,ω z ,δ;0) 3D TOF data B ℘(ω s ,φ,ω z ,δ;ω t ) ω x =ω s cosφ−χsinφ =ω ′ s cosφ ′ −ω t sinφ ′ 2D TOF data ω y =ω s sinφ+χcosφ =ω ′ s sinφ ′ +ω t cosφ ′ ℘(ω ′ s ,φ ′ ,ω z ,0;ω t ) 3D TOF data C ℘(ω s ,φ,ω z ,δ;ω t ) ω x =ω s cosφ−χsinφ =ω ′ s cosφ ′ 2D non TOF data ω y =ω s sinφ+χcosφ =ω ′ s sinφ ′ ℘(ω ′ s ,φ ′ ,ω z ,0;0) 3D non TOF data D ℘(ω s ,φ,ω z ,δ;0) ω x =ω s cosφ+δω z sinφ =ω ′ s cosφ ′ 2D non TOF data ω y =ω s sinφ−δω z cosφ =ω ′ s sinφ ′ ℘(ω ′ s ,φ ′ ,ω z ,0;0) 2D TOF data E ℘(ω s ,φ,ω z ,0;ω t ) ω x =ω s cosφ−ω t sinφ =ω ′ s cosφ ′ 2D non TOF data ω y =ω s sinφ+ω t cosφ =ω ′ s sinφ ′ ℘(ω ′ s ,φ ′ ,ω z ,0;0) Table 4.1: Summary of all mapping equations derived from the generalized projection slice theorem while preserving SNR benefits by implicitly utilizing the TOF information through the rebinning. All the mappings require to take Fourier transform in axial direction (z). This causes missing data problem [32,70] when implementing a rebinning algorithm. Some approx- imation was applied for Fourier rebinning (FORE) which can avoid this problem [32]. Exact rebinning method was also proposed where the inverse rebinning equation is used to estimate the missing part before applying the rebinning [70]. In Chapter 5, we also 79 address this issue, and show both exact and approximate TOF PET data rebinning meth- ods for a fast 3D TOF PET reconstruction. Finally, notice that our first work developing a fast projector described in Chapter 3 was based on the mapping B. 80 Chapter 5 Fourier Rebinning of 3D TOF PET Data for Fast Reconstruction 5.1 Introduction Fully 3D time-of-flight (TOF) PET scanners offer the potential for previously unachiev- able signal to noise ratio in clinical PET. However, fully 3D TOF PET image reconstruc- tion using accurate system and noise models is a challenging task due to the huge data size. So the rebinning techniques using the mapping equations derived in the previous chapter will be very useful in the 3D TOF PET data reconstruction. Previously, the sin- gle slice rebinning (SSRB-TOF) [94] was proposed where the oblique TOF sinograms are combined together to direct plane to form a set of stacked 2D TOF sinograms in a similar manner to SSRB [37] applied to 3D non TOF data. As an alternative to SSRB- TOF, Defrise and his colleagues proposed an approximate Fourier rebinning method, where the rebinning is performed in the frequency domain. This approximate approach shows improved performance over SSRB-TOF [28]. An exact rebinning equation was derived based on a consistency condition expressed by a partial differential equation in the continuous data domain [35], where rebinning is performed with respect to the axial variables. This result motivated the development of an approximate discrete axial rebinning method. In this method a cost function based on a bias and variance tradeoff is used to find optimal pre-computable rebinning coefficients. Using these coefficients, 81 a weighted average of the axial lines of response is taken to estimate an appropriate line of response in a 2D direct plane. All the rebinning methods reviewed above rebin 3D TOF data to 2D TOF data and specifically retain the TOF information in the rebinned data. We also present an exact rebinning method based on the mapping B in Figure 4.1 in Chapter 4. Then, we pro- pose some different rebinning methods where the 3D TOF data are rebinned to 2D or 3D non TOF data. Mapping A between 3D TOF data and 3D non TOF data, and map- ping C between 3D TOF data and 2D non TOF data will be used for those rebinnings. Mapping E is a special case of mapping A for 2D TOF PET data. Note that, compared to mapping B, the rebinnings A and C use the TOF information implicitly during the rebinning process so that the result of rebinning is not a function of the TOF informa- tion. However, since the mappings themselves do make use of the extra TOF informa- tion to convert them to the corresponding frequency components in non TOF PET data, the rebinned non TOF data are able to retain SNR advantages relative to count-matched data acquired without TOF information, or equivalently rebinnings computed by simply summing sinograms over the TOF bins. All mapping equations require to take Fourier transform in the axial direction. However, the oblique data are axially truncated due to the finite axial aperture of the scanner. To address the missing data problem, we apply approximate mapping equations as in [32]. The proposed rebinning methods are evaluated by simulations. This one time data reduction can be implemented in the frequency domain using the fast Fourier transform (FFT) algorithm, and once the 3D TOF PET data is reduced to their lower dimensional form, the data can be reconstructed in a much faster way com- pared to fully 3D TOF PET data reconstruction. Here, by fully 3D TOF PET data recon- struction, we mean an iterative 3D TOF PET data reconstruction where the iterations are performed between image and the fully 3D TOF PET data. The work developing a fast 82 projector in Chapter 3 can be also applied for reducing the reconstruction time of the fully 3D iterative TOF PET data reconstruction. 5.2 Rebinning of 3D Time of Flight (TOF) Data to 2D TOF Data Here we describe an exact rebinning method where 3D TOF PET data are rebinned to 2D TOF PET data. The exact rebinning method based on the mapping B in Chapter 4.3.2 requires both inverse rebinning and rebinning mappings. The inverse rebinning is required to estimate the missing part of oblique sinograms in axial direction for taking Fourier transform. After rebinning 3D TOF PET data to 2D TOF PET data, a slice-by- slice 2D TOF PET sinogram reconstruction is applied. 5.2.1 Method An exact inverse rebinning equation The inverse rebinning is used for estimating missing part in axial direction of oblique sinograms. Here, both mapping B and C, in Figure 4.1, can be used for estimating the missing part. The mapping C will have computational advantage compared to the mapping B since it can estimate the 3D TOF PET data from 2D non TOF PET data, and the non TOF 2D PET data can be obtained simply by summing the TOF sinograms in TOF bin direction. So the following mapping C is used for the inverse rebinning: ℘(ω s ,φ,ω z ,δ;ω t ) = √ 1+δ 2 {H (ω t )/H(0)}℘(ω s s 1+ χ ω s 2 ,φ+arctan χ ω s ,ω ′ z ,0;0) (5.1) 83 whereχ =ω t √ 1+δ 2 −δω z . An exact rebinning equation For the rebinning, we rebin the 3D TOF PET data to 2D TOF PET data, so we need to use the mapping B. As derived in Chapter 4.3.2, we use the following exact rebinning equation: ℘(ω ′ s ,φ ′ ,ω ′ z ,0;ω ′ t ) = 1 √ 1+δ 2 ℘(˜ ω ′ s ,φ ′ +α ′ −β ′ ,ω ′ z ,δ;ω ′ t ) (5.2) where χ ′ =ω ′ t √ 1+δ 2 −δω ′ z , ˜ ω ′ s =ω ′ s q 1+(ω ′ t 2 −χ ′2 )/ω ′ s 2 , β ′ = arctan(χ ′ /˜ ω ′ s ). This exact rebinning equation is defined only forχ ′2 ≤ω ′ s 2 +ω ′ t 2 . This constraint is because of the similar hole property of TOF PET data as indicated in Figure 4.9. Implementation To implement an exact rebinning mapping, we first perform inverse rebinning to esti- mate the missing oblique TOF sinograms to compute the Fourier transform in z. The direct plane sinograms are used as an input for the estimation process. After construct- ing a complete set of oblique sinogram data, the rebinning is performed based on (5.2) to map each of frequency component of the oblique sinogram data into the corresponding one in the 2D TOF PET data. A normalization step is required to account for variable contributions to each frequency bin in℘(ω ′ s ,φ ′ ,ω ′ z ,0;ω ′ t ) from the oblique sinograms in ℘(ω ′ s ,φ ′ ,ω ′ z ,δ;ω ′ t ) for all oblique angle variable δ as done in [70]. After the 3D TOF 84 Parameter Value Ring diameter (mm) 421 Detectors per ring 672 Number of rings 41 Rays (LORs) per angle 336 Maximum ring difference (MRD) 27 Span 3 TOF resolution (psec) 250 (or 500) Number of TOF bins 15 (or 7) Image size 256× 256× 81 Transverse voxel size (mm) 2 Table 5.1: Simulated 3D TOF cylindrical PET system parameters sinogram is rebinned to stacked 2D direct TOF sinograms, any 2D TOF image recon- struction method can be applied to each 2D direct TOF sinogram. We use a 2D direct Fourier reconstruction method based on the TOF projection slice theorem in (4.17). The exact inverse rebinning equation (5.1) and the exact rebinning equation (5.2) are reduced to the exact rebinning equations used in FOREX [70] whenω ′ t = ω t = 0. We therefore refer the TOF PET data rebinning method to TOF-FOREX. 5.2.2 Simulations We performed simulation studies to evaluate the exact 3D TOF PET data rebinning. We simulated a 3D TOF cylindrical PET system whose parameters are shown in Table 5.1. A 3D TOF projector was implemented using Siddon’s ray tracing algorithm [116] with a TOF kernel weight as in [47]. A Gaussian function was used for the TOF kernel. Two different TOF resolutions, 500 and 250 psecs, were tested. These values were used for both the full width at half maximum (FWHM) of the TOF kernel and for the TOF bin interval, leading to 7 and 15 TOF bins, respectively over the field of view. The 3D Hoffman digital brain phantom was used as the object for data generation; to better 85 illustrate the advantages of TOF data, we scaled the brain volume isotropically to the maximum size that can be contained within the transaxial field of view. 0 100 200 300 400 500 600 Reference TOF-FOREX Difference 4-th TOF bin 5-th TOF bin 6-th TOF bin 7-th TOF bin Figure 5.1: Comparison of rebinned TOF sinograms with corresponding directly com- puted TOF sinograms (TOF timing resolution 500 psec); Columns show TOF bin indices from 4 to 7; rows show the calculated direct sinogram (Reference), the rebinned sino- gram (TOF-FOREX), and their difference (Difference). To investigate a practical implementation of the exact rebinning algorithm, TOF- FOREX, we generated noiseless 3D TOF sinograms and rebinned ones using the exact rebinning described in Chapter 5.2.1. We then compared the rebinned stacked 2D TOF sinograms with those directly calculated from the phantom image as shown in Figure 5.1 and 5.2. The rebinned values agree closely with the directly calculated values. The nor- malized root-mean-squared (RMS) error between the two 2D direct sinograms was less than 6% for both 500 and 250 psec timing resolutions. The differences are probably due 86 50 100 150 200 250 0 50 100 150 200 250 300 350 400 450 500 Reference TOF−FOREX (a) 50 100 150 200 250 0 50 100 150 200 250 300 350 400 450 500 Reference TOF−FOREX (b) Figure 5.2: Comparison of rebinned TOF sinograms with corresponding directly com- puted TOF sinograms (TOF timing resolution 500 psec) shown in Figure 5.2; (a) radial profile through the sinograms for 4th TOF bin at 40-th angle; (b) radial profile through the sinograms for 7th TOF bin at 40-th angle 87 0.00 5.00 0.00 1.32 (a) TOF resolution: 500 psec 0.00 6.00 0.00 1.91 (b) TOF resolution: 250 psec Figure 5.3: Comparison of (top) noisy 2D direct sinogram with (bottom) noise reduction resulting from rebinning of oblique data into direct planes for (a) 500 psec, and (b) 250 psec resolution. to a combination of numerical differences when rebinning compared to direct calcula- tion, and the impact of the small degree of aliasing reported above. Next we generated noisy 3D TOF sinograms based on a total of 200 million counts for TOF timing resolutions of 500 and 250 psecs. Figure 5.3 shows the SNR improve- ment resulting from rebinning of all oblique data into a set of stacked 2D sinograms. After rebinning, we reconstructed 2D images from the noisy rebinned data using 2D direct Fourier reconstruction based on the projection slice theorem in (4.7). We used a gridding method for 2D interpolation based on a Kaiser-Bessel convolution kernel [58] to obtain frequency samples in the Cartesian grid before taking the 2D inverse Fourier Transform to reconstruct each slice. For comparison, a 3D non-TOF sinogram with the same total number of counts was also generated and FOREX applied. We compared the reconstructed images at the axial center plane in Figure 5.4. As expected, we see 88 −0.5 0 0.5 1 1.5 2 2.5 (a) (b) (c) (d) Figure 5.4: 3D reconstructions from 200 million counts for the central axial slice: (a) true image; (b) 3D non-TOF reconstruction: FOREX followed by direct 2D Fourier reconstruction; (c) 3D TOF reconstruction using TOF-FOREX followed by 2D TOF direct Fourier reconstruction with TOF timing resolutions of 500 psec: (d) as for (c) but with 250 psec resolution. the SNR gains resulting from the use of TOF data compared to non-TOF data with a matched number of counts, and further improvements when increasing timing resolu- tion from 500 to 250 psec. Although TOF-FOREX is exact, it requires more computation compared to the approximated TOF-FORE in [28] due to the inverse rebinning step required to estimate missing planes. Consequently, in the following, we propose more practical versions of 89 rebinning methods, where 3D TOF PET data are rebinned to 2D or 3D non TOF PET data, and approximated to avoid the inverse rebinning step. 5.3 Rebinning of 3D TOF Data to 2D and 3D non TOF Data Previously, we described an exact 3D TOF PET data rebinning method where 3D TOF PET data are rebinning to 2D TOF PET data similar to [28,35,94]. But, here we want to propose other interesting rebinning methods where the 3D TOF PET data are rebinned to non TOF PET data using the mapping A and C in Figure 4.1. The TOF information is used implicitly during the rebinning process so that the result of rebinning is not a function of the TOF information. However, since the mappings themselves utilizes the TOF information, the rebinned non TOF data are able to retain SNR advantages relative to count-matched data acquired without TOF information. We already showed exact mapping equations for rebinning 3D TOF data to 3D or 2D non TOF data in the previous chapter, and both equations require Fourier transform in the axial direction. However, the oblique sinogram data are axially truncated due to the finite axial aperture of the scanner. To address the missing data problem, we apply some approximations to change the exact mapping equations into the ones which do not require the Fourier transform in the axial direction as in [32]. First, we investigate the accuracy of the approximate rebinning methods by compar- ing the error induced by the approimations for two mapping A and C. While the mapping C produces more data reduction factor, it shows larger approximation error compared to the approximation applied to the mapping C. So the approximate rebinning method using the mapping A was implemented, and evaluated by simulations. Although a fully 3D non TOF PET data reconstruction is required after rebinning using the mapping A, 90 the compution can be significantly reduced by using fast projectors [24, 27, 51, 80] or dedicated hardware [13, 102, 122]. We also conduct 2D Monte Carlo simulations to examine the noise properties of the rebinning. Finally, we perform 3D TOF PET recon- struction from noisy data using the rebinning from 3D TOF to 3D non TOF and compare results with those from using non TOF 3D data. 5.3.1 Non TOF Fourier Rebinning and Approximation Since the approximate TOF rebinning methods described below use similar approxima- tions to those for the non TOF case, we first review the exact and approximate Fourier rebinning for non TOF data [29, 32]. The exact inverse Fourier rebinning equation [32] can be written as ℘(ω s ,φ,ω z ,δ;0) = √ 1+δ 2 ℘(ω s s 1+ δω z ω s 2 ,φ−arctan δω z ω s ,ω z ,0;0), (5.3) which was derived using the Mapping D in Table 4.1 in the previous chapter. The exact mapping equation requires taking Fourier transform inz direction. How- ever, the oblique sinograms in 3D PET data are not axially complete since the data are truncated due to the finite axial aperture of the scanner. So the missing part needs to be estimated before applying the rebinning [70]. Alternatively, one can use the first order Taylor series truncation with respect toδω z /ω s , ω ′ s =ω s s 1+ δω z ω s 2 ≈ω s , φ ′ =φ−arctan δω z ω s ≈φ− δω z ω s , (5.4) to obtain the Fourier rebinning (FORE) equation [29, 32] P(ω s ,k,z,δ;0)≈ √ 1+δ 2 P(ω s ,k,z− δk ω s ,0;0) (5.5) 91 B φ ′ = φ+α ω ′ s = ω s v u u u u t1+   δω z ω s   2 α = arctan   −δω z ω s   φ −δω z ω x ω y ω s = 0 A C |AB| =|ω s | |CB| =|ω ′ s | 6 ABC =α |ω ′ s |≈|AB| α≈ −δω z ω s Figure 5.5: Graphical interpretation of the exact and approximate inverse Fourier rebin- ning mappings from 2D non TOF data to 3D non TOF data. Inverse Fourier rebinning is equivalent to finding the coordinates (ω ′ s ,φ ′ ) of a point (for example, point B) on a trajectory (for example, the line joining the points A and B) mapped from the 3D sinogram. In the approximate inverse rebinning,|CB| and∠ABC are approximated by |AB| and−δω z /ω s , respectively. where k is the Fourier series index corresponding to φ andP is the 2D Fourier trans- form of p(s,φ,z,δ;0) in s and φ. Then the approximation in (5.5) does not require the Fourier transform in the z direction. The approximation error becomes larger as |δω z /ω s | increases. To reduce the side effect caused by the approximation errors, Defrise et al used oblique sinograms only for smallδ when|ω s | is small [32]. Figure 5.5 is a modified version of Figure 4.9 by adding the remarks for the approx- imation. 92 5.3.2 TOF Rebinning Mappings and Approximations We now consider mappings A and C in Table 4.1, which rebin 3D TOF data to 3D and 2D non TOF data, respectively. To address the missing data problem, we propose approximations, similar to the non TOF case, which do not require the Fourier transform inz. Approximate mapping between 3D TOF and 2D non TOF data (Mapping C) Noting that the mapping equations for the mapping C in Table 4.1 are the same as those for the mapping D of non TOF data except for the−δω z andχ terms, we use a first order Taylor series truncation with respect toχ/ω s , similar to (5.4), ω ′ s =ω s s 1+ χ ω s 2 ≈ω s , ˜ ω ′ s , φ ′ =φ+arctan χ ω s ≈φ+ χ ω s , ˜ φ ′ . (5.6) This approximation leads to the following approximate inverse rebinning equation: P(ω s ,k,z,δ;ω t )≈ √ 1+δ 2 {H(ω t )/H(0)}e ik „ ω t √ 1+δ 2 ωs « P(ω s ,k,z− δk ω s ,0;0), (5.7) which does not require a Fourier transform in thez direction. Similarly, an approximate rebinning equation can be written as P(ω s ,k,z,0;0)≈ 1 √ 1+δ 2 {H(0)/H(ω t )}e −ik „ ω t √ 1+δ 2 ωs « P(ω s ,k,z + δk ω s ,δ;ω t ). (5.8) Figure 5.5 also applies to this case if we useχ in place of−δω z . The errors due to the approximations in (5.6) are comparable to those from (5.4) for FORE, as we show in Chapter 5.3.2. 93 Approximate mapping between 3D TOF data and 3D non TOF data (Mapping A) The exact equations mapping between 3D TOF and 3D non TOF data are given in (4.49). To remove the dependency onω z , we make the following approximation, ω ′ s =ω s s 1+ χ 2 −(δω z ) 2 ω 2 s =ω s v u u t 1+ ω t √ 1+δ 2 ω s ! 2 −2 ω t √ 1+δ 2 ω s ! δω z ω s ≈ω s r 1+ ω t √ 1+δ 2 /ω s 2 , ˜ ω ′ s φ ′ =φ+arctan ω t √ 1+δ 2 −δω z ω s ! +arctan δω z ω ′ s ≈φ+arctan ω t √ 1+δ 2 /ω s , ˜ φ ′ , (5.9) which can be seen as the zeroth order Taylor series truncation with respect to δω z /ω s . This approximation yields the following approximate inverse rebinning equation: ℘(ω s ,φ,z,δ;ω t )≈ H(ω t ) H(0) ℘ ω s r 1+ ω t √ 1+δ 2 /ω s 2 ,φ+arctan ω t √ 1+δ 2 /ω s ,z,δ;0 ! . (5.10) The corresponding approximate rebinning equation can be written as ℘(ω s ,φ,z,δ;0)≈ H(0) H(ω t ) ℘ ω s r 1− ω t √ 1+δ 2 /ω s 2 ,φ−arctan ω t √ 1+δ 2 /ω ′ s ,z,δ;ω t ! , (5.11) whereω ′ s =ω s q 1− ω t √ 1+δ 2 /ω s 2 . 94 ω x ω y α ω t √ 1+δ 2 B −δω z α φ φ ′ = φ+α ω ′ s = ω s v u u u u t1+ χ 2 −(δω z ) 2 ω 2 s χ |AB| =|ω s | |CB| =|ω ′ s | A C D |ω ′ s |≈|DB| 6 ABC =α 6 ABD =α approx α≈ α approx α approx = arctan    ω t √ 1+δ 2 ω s    |DB| =|ω s | v u u u u u t1+    ω t √ 1+δ 2 ω s    2 α approx Figure 5.6: Graphical interpretation of the exact and approximate inverse Fourier rebin- ning mappings from 3D non TOF data to 3D TOF data. The exact inverse Fourier rebinning is equivalent to finding the coordinate information of the line segment CB (|CB| and∠ABC) for the line segment AB. By the approximation in (5.9),|CB| and ∠ABC are approximated by|DB| and∠ABD, respectively. The approximation error will increase with|δω z /ω s |, and also with|ω t | due to the term ω t √ 1+δ 2 /ω s ignored when approximatingω s ′ in (5.9). Note that the approxi- mations in (5.10) and (5.11) are exact for direct sinograms (δ = 0), and the approximate rebinning method uses fixedδ andz. In other words, the non TOF sinogram for a given (z,δ) is computed only from the TOF data for the same axial planez with the oblique angleδ. Figure 5.6 is illustrating the approximation. The point B corresponds to the 3D TOF PET data℘(ω s ,φ,ω z ,δ;ω t ) in the(ω x ,ω y ) plane. The inverse rebinning is equivalent to finding the 3D non TOF data℘(ω ′ s ,φ ′ ,ω z ,δ;0) whereω ′ s denotes the line segment from 95 mapping C mapping A ω exact x ω ′ s cosφ ′ ω ′ s cosφ ′ +δω z sinφ ′ ω exact y ω ′ s sinφ ′ ω ′ s sinφ ′ −δω z cosφ ′ ω approx x ˜ ω ′ s cos ˜ φ ′ ˜ ω ′ s cos ˜ φ ′ +δω z sin ˜ φ ′ ω approx y ˜ ω ′ s sin ˜ φ ′ ˜ ω ′ s sin ˜ φ ′ −δω z cos ˜ φ ′ Table 5.2: Exact and approximate points in (ω x ,ω y ) plane for the inverse rebinning for mappings C and A. point C to B and φ ′ is the sum of φ and∠ABC. The approximation in (5.9) is to take point D in place of point C in Figure 4.10. That is, ω ′ s and φ ′ are approximated by the coordinate information (|DB| andα approx ) of the line segment from point D to B. Comparison of Approximation Errors We have described two approximate rebinning mappings (5.8) and (5.11). The inverse rebinning equations are used to estimate at each frequency the 3D TOF data ℘(ω s ,φ,z,δ;ω t ) from 2D non TOF data ℘(ω s ,φ,z,0;0) or 3D non TOF data ℘(ω s ,φ,z,δ;0). To quantify the errors from the approximations in (5.6) and (5.9), we calculated the Euclidean distance between the exact and approximate points used for estimating the℘(ω s ,φ,ω z ,δ;ω t ) data in the (ω x ,ω y ) plane as shown in Figure 4.9 and Figure 4.10. The error measure is defined as E = q (ω exact x −ω approx x ) 2 + ω exact y −ω approx y 2 (5.12) whereω exact x ,ω exact y ,ω approx x andω approx y are defined in Table 5.3.2 for the mappings C and A. Figure 5.7 and 5.8 show the approximation errors averaged overω z andω t at four dif- ferent ring differences as a function of transaxial radial frequencyω s for a fixedφ. The 96 0 20 40 60 80 100 120 140 0 1 2 3 4 5 6 ω s frequency index Error Mapping A Mapping C Mapping D (FORE) (a) Ring difference: 10 0 20 40 60 80 100 120 140 0 1 2 3 4 5 6 7 8 9 ω s frequency index Error Mapping A Mapping C Mapping D (FORE) (b) Ring difference: 25 Figure 5.7: Averaged mapping errors of approximate rebinnings for different ring dif- ferences and for a fixedφ = 10. The error value is normalized by the image frequency sample interval (Δω x ) in the transverse plane where Δω x = 1/2/256 mm −1 in our sim- ulation. (See Table 5.3.3 for parameters). Notice the error is very large when ω s is small, so special consideration will be required in this region, similarly to the FORE implementation. 97 0 20 40 60 80 100 120 140 0 2 4 6 8 10 12 ω s frequency index Error Mapping A Mapping C Mapping D (FORE) (a) Ring difference: 35 0 20 40 60 80 100 120 140 0 2 4 6 8 10 12 14 16 ω s frequency index Error Mapping A Mapping C Mapping D (FORE) (b) Ring difference: 50 Figure 5.8: Averaged mapping errors of approximate rebinnings for different ring dif- ferences and for a fixedφ = 10. The error value is normalized by the image frequency sample interval (Δω x ) in the transverse plane where Δω x = 1/2/256 mm −1 in our sim- ulation. (See Table 5.3.3 for parameters). Notice the error is very large when ω s is small, so special consideration will be required in this region, similarly to the FORE implementation. 98 FORE case (mapping D) is also included in the figure for comparison. All the approxi- mate mappings tend to have larger errors as the oblique angle (arctanδ) increases. The approximate mapping A from 3D TOF data to 3D non TOF data shows smaller errors than the other approximate mappings. The approximation errors for the mapping C from 3D TOF data to 2D non TOF data are comparable to that for the conventional Fourier rebinning (FORE). The figure shows that the approximation errors are large when ω s is small. Therefore, as in the FORE implementation, special treatment is required for smallω s , as described in the next section. 5.3.3 Results We performed simulation studies in the geometry of Siemens Biograph PET/CT True Point TrueV scanner [57], a state-of-the-art clinical scanner with planned TOF capa- bilities. The scanner geometry and simulation parameters are shown in Table 5.3.3. A Gaussian TOF kernel function was used for TOF data generation and 3D projections computed based on line integrals as in [47]. A TOF resolution of 500 ps was used and the data were sampled with the sampling period of 250 ps, leading to 13 TOF bins over the field of view. The NCAT torso phantom [115] was used as a 3D object. Attenuation, randoms, scatters and detector efficiencies were not considered. Here we evaluate the performance of the approximate rebinning method based on (5.11), which rebins 3D TOF data to 3D non TOF data, since this has less approximation error than the rebinning of 3D TOF to 2D non TOF data in (5.8), as shown in Figure 5.7 and 5.8. For each ω t , we estimate 3D non TOF data℘(ω s ,φ,z,δ;0) from the 3D TOF data using the rebinning equation (5.11) and calculate an average over all ω t ’s. As mentioned previously, the approximation error is large when |ω s | is small and |ω t | is large. To reduce such errors as in [32], when|ω s | is less than a thresholdS TH , we only 99 Parameter Value Ring radius (mm) 421 Detectors per ring 672 Number of rings 55 Rays (LORs) per angle 336 Maximum ring difference (MRD) 54 Span 5 TOF resolution (ps) 500 Number of TOF bins 13 Image size 256× 256× 109 Transverse voxel size (mm) 2 Scanner axial FOV (mm) 216 Table 5.3: Simulated 3D TOF cylindrical PET system parameters use the data℘(ω s ,φ,z,δ;ω t ) for|ω t | smaller than a threshold T TH . In our simulation, S TH = 7 andT TH = 1 were used. We did not implement the other rebinning equation (5.8) yet. An extensive compar- ison of different rebinning strategies will be left as future work. Comparison of Rebinned Data and non TOF Data To evaluate the performance of the approximate rebinning algorithm, which we refer to as FORET (FOurier REbinning of Time of flight), we compare rebinned data and non TOF data by Monte Carlo simulation. Poisson noise was added to noiseless 3D TOF data to generate 500 noisy realizations. Each noisy data set was rebinned by FORET to obtain rebinned sinograms. The approximate rebinning algorithm was applied to TOF datap(s,φ,z,δ;t) independently for each(z,δ). Non TOF data was obtained by simply summing the 3D TOF data in TOF bin direction. The mean and variance of the rebinned data and non TOF data were calculated and compared. Results for planes at axial center for two different ring differences are shown in figure 5.9 and 5.10. Figure 5.9(a) shows the mean of the rebinned sinogram at zero ring difference. Since the approximate rebinning is exact when the ring difference is 100 zero, it is unbiased. Figure 5.10(a) shows the profiles of the mean for the maximum ring difference. Even at this maximum ring difference there appears to be no significant approximation error. Figure 5.9(b) and Figure 5.10(b) show the profiles of the variance of the rebinned and non TOF sinograms at two different ring differences. As the figures show, significant variance reduction was achieved by the rebinning when compared to simply summing over the TOF bins. This is a perhaps surprising result, since each non TOF sinogram for a fixed (z,δ) is computed from the same data for the two methods compared in these figures: in one case, by summing over the TOF bins, in the other by rebinning using (5.11). Monte Carlo Simulation for 2D TOF PET Data Reconstruction To conclude we compare reconstructed images from non TOF data, TOF data rebinned to non TOF, and unrebinned TOF data. For the Monte Carlo study we consider the 2D case. To perform the study, we used the central plane of the stacked 2D sinograms (zero ring difference) from the previous section. Images were reconstructed by the MAP method [109]. Three cases were compared. First, the noisy 2D TOF data were used for TOF MAP reconstruction (denoted by ‘TOF’). Second, we rebinned the TOF data to 2D non TOF data by FORET (denoted by ‘FORET’) and then used the rebinned data for MAP reconstruction. Lastly, the noisy TOF data was summed over the TOF bins to generate non TOF data (denoted by ‘non TOF’) and then used for MAP reconstruction. To quantify the resolution and noise properties of the reconstructed images, we cal- culated the full-width-at-half-maximum (FWHM) of the local impulse response [43] and the sample variance of noisy reconstructed images for different regularization parame- ters. Figure 5.11 and 5.12 show the results at 4 different pixel locations. The results show that ‘FORET’ gives a better resolution vs. noise trade-off compared to the non TOF case (‘non TOF’) and its performance is somewhat worse than that for the TOF 101 (a) 50 100 150 200 250 0 2 4 6 8 10 12 14 16 18 20 Variance Sinogram radial index non TOF FORET (b) Figure 5.9: A comparison of the mean and variance of rebinned sinograms, obtained by FORET, and non TOF sinograms for two ring differences: (a) Mean of rebinned sinograms for RD (ring difference) = 0, (b) Variance profiles of rebinned and non TOF sinograms at 160-th angle (RD=0). 102 50 100 150 200 250 0 2 4 6 8 10 12 14 16 18 20 Mean Sinogram radial index non TOF FORET (a) 50 100 150 200 250 0 2 4 6 8 10 12 14 16 18 20 Variance Sinogram radial index non TOF FORET (b) Figure 5.10: A comparison of the mean and variance of rebinned sinograms, obtained by FORET, and non TOF sinograms for two ring differences: (a) Mean profiles of rebinned and non TOF sinograms at 160-th angle (RD=54), (b) Variance profiles of rebinned and non TOF sinograms at 160-th angle (RD=54). 103 case (‘TOF’). Here, the MAP reconstruction method assumed a Poisson model for the ‘FORET’ case although in practice the Poisson statistics are destroyed by rebinning as discussed in [22, 69]. We expect that reconstruction using more accurate statistical models will enhance the performance for rebinned data. 3D TOF Data Reconstruction Finally, we show the performance of ‘FORET’ and ‘non TOF’ for 3D TOF PET recon- struction. A fully 3D MAP reconstruction method [109] was applied to 3D non TOF PET data rebinned from 3D TOF PET data. The total number of counts in the data was 20M. In order to match the resolution, the regularization parameters were adjusted so that the local impulse response at the image center has a FWHM of 3.4mm in the trans- verse place. Figure 5.13 shows the reconstructed images in three orthogonal views. As the figures show, the 3D MAP reconstruction with the approximate rebinning yields less noisy images compared to the non TOF 3D reconstruction case. 5.4 Conclusions Our all derivations were based on the analytical model, and taking Fourier transform in TOF bin direction. Clearly, current TOF timing resolution limits the corresponding spatial resolution along the LOR in the time of flight variable. While we were concerned that this would limit our ability to use Fourier transforms int, in practice the bandlimit- ing effect of the TOF kernel appears to ameliorate potential aliasing effects. This appears to be the case even at currently achievable timing resolutions of 500 psec [19,117,118]. We proposed first an exact TOF PET data rebinning method TOF-FOREX. It requires more computation compared to the approximated TOF-FORE in [28] due to the inverse rebinning step required to estimate missing planes. The issue was addressed 104 3.5 4 4.5 5 5.5 6 6.5 0 0.01 0.02 0.03 0.04 0.05 0.06 FWHM (mm) Pixel Variance Variance comparison at (128,128) pixel non TOF FORET TOF (a) 3.5 4 4.5 5 5.5 6 6.5 0 0.01 0.02 0.03 0.04 0.05 0.06 FWHM (mm) Pixel Variance Variance comparison at (128,160) pixel non TOF FORET TOF (b) Figure 5.11: Resolution (FWHM) versus pixel variance plot for four different locations. (a) pixel location (128,128); (b) pixel location (128,160) 105 3.5 4 4.5 5 5.5 6 6.5 0 0.01 0.02 0.03 0.04 0.05 0.06 FWHM (mm) Pixel Variance Variance comparison at (128,220) pixel non TOF FORET TOF (a) 3.5 4 4.5 5 5.5 6 6.5 0 0.01 0.02 0.03 0.04 0.05 0.06 FWHM (mm) Pixel Variance Variance comparison at (110,110) pixel non TOF FORET TOF (b) Figure 5.12: Resolution (FWHM) versus pixel variance plot for four different locations. (a) pixel location (128,220); (b) pixel location (110,110) 106 (a) (b) (c) (d) (e) (f) Figure 5.13: A comparison of 3D TOF data reconstruction by ‘FORET+MAP’ (top row) and ‘non TOF+MAP’ (bottom row) in the transverse view (first column), coronal view (second column) and sagittal view (third column). in our alternative rebinning methods where the 3D TOF data are rebinned to either 3D or 2D non TOF data. This method is different from previous rebinning meth- ods [28, 35, 94, 130] where the 3D TOF data are rebinned to the 2D TOF data. Our approximations allow practical rebinnings of TOF to non TOF data. Using simulated 3D TOF data, we evaluated the performance of an approximate rebinning from 3D TOF data to 3D non TOF data, which has less approximation error than the method rebinning 3D TOF data to 2D non TOF data. Although a fully 3D non TOF PET data reconstruc- tion is computationally more demanding than 2D non TOF PET data reconstruction, fast projectors [24, 27, 51, 80] or dedicated hardware [13, 102, 122] significantly reduce the cost. 107 Monte Carlo simulations showed that the approximate rebinning method produced improved SNR compared to data without TOF information, but with some loss in per- formance relative to full use of the TOF data without rebinning. 108 Chapter 6 Summary and Future Work 6.1 Summary Our work was to reduce the high computation load involved in fully 3D PET iterative image reconstruction and 3D TOF PET image reconstruction. As our first work, we developed a fast projector using the inverse (Fourier) rebinning calculating 3D PET data from 2D PET data. The inverse rebinning operator was combined with a 2D projector to construct a fast 3D projector which maps 3D image to 3D PET data. The inverse rebin- ning operator was implemented cost-effectively using FFT, so the computation time was reduced by an order of magnitude. We utilized the fast projector in the context of our fully 3D MAP iterative image reconstruction, and achieved significant reconstruc- tion time savings. A slight resolution loss was observed in the reconstruction using the fast projector, however, the image quality degradation was hardly visible in the images reconstructed from real in vivo data. Next, we extended the Fourier rebinning concept to 3D TOF PET image reconstruc- tion. We derived a projection slice theorem based on a line integral based data model using a spatially invariant TOF kernel function. In the data model, each TOF PET data bin is a convolution of the source activity function with the TOF kernel function through an LOR. Under the assumption of the spatial invariance for the TOF kernel function, we take Fourier transform in TOF variable to derive the projection slice theorem. Taking the Fourier transform in the TOF variable may cause aliasing since the sampling frequency in the TOF variable is very low. However, the bandwidth along the LOR is limited 109 through the convolution with the TOF kernel which limits the aliasing effects. The gen- eralized projection slice theorem derived in the frequency domain gives the relationship between the Fourier transform of 3D object function and the Fourier transform of 3D TOF PET data. Based on the theorem, a unified framework was built to derive mapping equations between different data sets, such as 3D TOF, 3D non TOF, 2D TOF and 2D non TOF PET data. In the framework, all the mapping equations can be interpreted geo- metrically in the frequency domain. The framework also includes the previous Fourier rebinning applied to non TOF 3D PET data. Using the mapping equations, we imple- mented an exact rebinning method where 3D TOF PET data are rebinned to 2D TOF PET data. We also proposed other rebinning techniques where 3D TOF PET data are rebinned to 2D or 3D non TOF PET data. Some approximations are taken to avoid the missing data problem. The rebinning methods were evaluated by simulations, and showed superior performance compared to the non TOF case where the TOF informa- tion is ignored. 6.2 Future Work The future work is listed as follows. • Real data study We need to evaluate the proposed TOF PET data rebinning meth- ods with real TOF PET data to evaluate the performance of our rebinning methods. For the real data study, we also need to consider many complicating factors, such as detector efficiency, attenuation, randoms and scatters. • Information theoretic analysis We want to further investigate optimal rebinning methods. An information theoretical analysis will help to find an optimal rebin- ning method, to evaluate how our current rebinning methods compared to the optimal rebinning. 110 • Developing efficient iterative reconstruction methods for rebinned data After rebinning is applied, the rebinned data do not follow the Poisson distribution an y more and the data become correlated. It will improve the image quality if the accurate statistical model is incorporated into reconstruction methods. • Comparison study of our TOF PET data rebinning algorithms It may be also worth conducting comparison studies for our TOF rebinning methods with other published approximated rebinning methods [28, 34, 94]. • Designing a fast projector for fully 3D TOF PET iterative reconstruction We will use the inverse rebinning equation to design a fast fully 3D TOF PET projector for fully 3D TOF PET iterative reconstruction as an extension of the non TOF PET data case in Chapter 3. 111 References [1] S. Ahn and J. A. Fessler, “Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms,” IEEE Trans. Med. Imag., 22(3):613–626, May 2003. [2] S. Ahn, J. A. Fessler, D. Blatt, and A. O. Hero, “Convergent incremental opti- mization transfer algorithms: Application to tomography,” IEEE Trans. Med. Imag., 25(3):283–96, March 2006. [3] A. M. Alessio, P. E. Kinahan, R. L. Harrison, and T. K. Lewellen, “Measured spatially variant system response for PET image reconstruction,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., volume 4, pages 1986–90, 2005. [4] A. M. Alessio, P. E. Kinahan, and T. K. Lewellen, “Modeling and incorporation of system response functions in 3-D whole body PET,” IEEE Trans. Med. Imag., 25:828–37, 2006. [5] H. O. Anger, “Survey of radioisotope cameras,” ISA Trans., 5:311–34, 1966. [6] E. Asma, D. W. Shattuck, and R. M. Leahy, “Lossless compression of list-mode 3D PET data,” IEEE Trans. Nucl. Sci., 50:9–16, 2003. [7] S. Basu and Y . Bresler, “O(N 2 logN) filtered backprojection reconstruction algo- rithm for tomography,” IEEE Trans. Image Processing, 9(10):1760–73, October 2000. [8] A. Boag and Y . Bresler, “A multilevel domain decomposition algorithm for fast O(N 2 logN) reprojection tomographic images,” IEEE Trans. Image Processing, 9:1573–81, September 2000. [9] S. Basu and Y . Bresler, “O(N 3 logN) backprojection algorithm for the 3-D Radon transform,” IEEE Trans. Med. Imag., 21(2):76–88, February 2002. 112 [10] F. B. Bouall` egue, J. Crouzet, C. Comtat, M. Fourcade, B. Mohammadi, and D. Mariano-Goulart, “Exact and approximate Fourier rebinning algorithms for the solution of the truncation problem in 3-D PET,” IEEE Trans. Med. Imag., 26:1001–9, 2007. [11] J. A. Browne and A. R. De Pierro, “A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography,” IEEE Trans. Med. Imag., 15(5):687–699, October 1996. [12] R. Blasberg, “PET imaging of gene expression,” E. Journal of Cancer, 38:2137– 46, 2002. [13] B. Bai and A. M. Smith, “Fast 3D iterative reconstruction of PET images using PC graphics hardware,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 2787–90, 2006. [14] T. F. Budinger, “Time-of-flight positron emission tomography: status relative to conventional PET,” J. Nuc. Med., 24:73–6, 1983. [15] T. Bruckbauer, K. Wienhard, S. B. Hansen, L. Eriksson, and G. Blomquist, “Eva- lutation of the ECAT EXCAT HR with ASCII for clinical routine 3D,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., volume 3, pages 1378–82, 1995. [16] C. L. Byrne, “Accelerating the EMML algorithm and related iterative algorithms by rescaled block-iterative methods,” IEEE Trans. Image Processing, 7(1):100– 109, January 1998. [17] S. Cho, S. Ahn, Q. Li, and L. M. Leahy, “Analytical properties of time-of-flight PET data,” Phys. Med. Biol., 53:2809–21, 2008. [18] S. Cho, S. Ahn, Q. Li, and L. M. Leahy, “Exact and approximate rebinning of PET data from time-of-flight to non time-of-flight,” Phys. Med. Biol., 53, 2008. submitted. [19] M. Conti, B. Bendrien, M. Casey, M. Chen, F. Kehren, C. Michel, and V . Panin, “First experimental results of time-of-flight reconstruction on an LSO PET scan- ner,” Phys. Med. Biol., 50:4507–26, 2005. [20] S. R. Cherry, M. Dahlbom, and E. J. Hoffman, “3D PET using a conventional multislice tomograph without septa,” J. Comput. Assist. Tomogr., 15(4):655–68, July 1991. [21] S. R. Cherry, “The 2006 Henry N. Wagner lecture: Of mice and men (and positrons)-Advances in PET imaging technology,” J. Nuc. Med., 47(11):1735– 45, November 2006. 113 [22] C. Comtat, P. E. Kinahan, M. Defrise, C. Michel, and D. W. Townsend, “Fast reconstruction of 3D PET data with accurate statistical modeling,” IEEE Trans. Nucl. Sci., 45(3):1083–9, June 1998. [23] R. Clack, “Towards a complete description of three-dimensional filtered back- projection,” Phys. Med. Biol., 37(3):645–60, 1992. [24] S. Cho, Q. Li, S. Ahn, B. Bai, and R. M. Leahy, “Iterative image reconstruction using inverse Fourier rebinning for fully 3-D PET,” IEEE Trans. Med. Imag., 26(5):745–56, May 2007. [25] J. G. Colsher, “Fully three dimensional positron emission tomography,” Phys. Med. Biol., 25(1):103–15, January 1980. [26] A. Chatziioannou, J. Qi, A. Annala, A. Moore, R. M. Leahy, and S. R. Cherry, “Comparison of 3D MAP and FBP algorithms for image reconstruction in microPET,” IEEE Trans. Med. Imag., 19(5):507–12, May 2000. [27] B. De Man and S. Basu, “Distance-driven projection and backprojection in three dimensions,” Phys. Med. Biol., 49:2463–75, 2004. [28] M. Defrise, M. E. Casey, C. Michel, and M. Conti, “Fourier rebinning of time- of-flight PET data,” Phys. Med. Biol., 50:2749–63, 2005. [29] M. Defrise, “A factorization method for the 3-D X-ray transform,” Inverse Prob., 11:983–94, April 1995. [30] M. Defrise, A. Geissbuhler, and D. W. Townsend, “A performance study of 3D reconstruction algorithms for positron emission tomography,” Phys. Med. Biol., 39(3):305–20, March 1994. [31] M. Defrise and P. Kinahan, “Data acquisition and image reconstruction for 3D PET,” In B Bendriem D W Townsend, editor, The theory and practice of 3D PET, pages 11–53. Kluwer, Netherlands, 1998. [32] M. Defrise, P. E. Kinahan, D. W. Townsend, C. Michel, M. Sibomana, and D. F. Newport, “Exact and approximate rebinning algorithms for 3-D PET data,” IEEE Trans. Med. Imag., 16:145–58, April 1997. [33] M. Defrise and X. Liu, “A fast rebinning algorithm for 3D positron emission tomography using John’s equation,” Inverse Prob., 15(4):1047–65, August 1999. [34] M. Defrise, V . Panin, C. Michel, and M. E. Casey, “Discrete axial rebinning for time-of-flight PET,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., 2006. 114 [35] M. Defrise, V . Panin, C. Michel, and M. E. Casey, “Continuous and discrete data rebinning in time-of-flight PET,” IEEE Trans. Med. Imag., 27, 2008. in press. [36] M. Defrise, D. W. Townsend, and R. Clack, “Three-dimensional image recon- struction from complete projections,” Phys. Med. Biol., 34(5):573–87, 1989. [37] M. E. Daube-Witherspoon and G. Muehllehner, “Treatment of axial data in three- dimensional PET,” J. Nuc. Med., 28:1717–24, 1987. [38] A. R. De Pierro and M. E. B. Yamagishi, “Fast EM-like methods for maxi- mum ‘a posteriori’ estimates in emission tomography,” IEEE Trans. Med. Imag., 20(4):280–288, April 2001. [39] P.R. Edholm, R. M. Lewitt, and B. Lindholm, “Novel properties of the Fourier decomposition of the sinogram,” In Proc. of the SPIE, Int. Workshop on Physics and Engineering of Computerized Multidimensional Imaging and Processing, pages 8–18, 1986. [40] F. H. Fahey, “Data acquisition in PET imaging,” Journal of Nucl. Med. Technol., 30:39–49, 2002. [41] T. H. Farquhar, A. Chatziioannou, and S. R. Cherry, “An evaluation of exact and approximate 3-D reconstruction algorithms for a high-resolution, small-animal PET scanner,” IEEE Trans. Med. Imag., 17:1073–80, 1998. [42] J. A. Fessler, “Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography,” IEEE Trans. Image Processing, 5(3):493–506, March 1996. [43] J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized- likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Processing, 5(9):1346–1358, September 1996. [44] T. Frese, N. C. Rouze, C. A. Bouman, K. Sauer, and G. D. Hutchins, “Quanti- tative comparison of FBP, EM, and Bayesian reconstruction algorithms for the IndyPET scanner,” IEEE Trans. Med. Imag., 22:258–76, 2003. [45] A. Gunawardana and W. Byrne, “Convergence of EM variants,” Technical Report CLSP Research Note No. 32, ECE Dept., Johns Hopkins University, February 1999. [46] E. S. Garnett, G. Firnau, and C. Nahmias, “Dopamine visualized in the basal ganglia of living man,” Nature, 305:137–8, 1983. 115 [47] C. J. Groiselle and S. J. Glick, “3D PET list-mode iterative reconstruction using time-of-flight information,” In Proc. IEEE Nuclear Science Symp. Medical Imag- ing Conf., pages 2633–8, 2004. [48] S. T. Grafton, J. C. Mazziotta, S. Presty, K. J. Frston, R. S. J. Frackowiak, and M. E. Phelps, “Functional anatomy of human procedural learning determined with regional cerebral blood flow and PET,” J. Neuroscience, 12:2542–8, 1992. [49] D. Gourion and D. Noll, “The inverse problem of emission tomography,” Inverse Prob., 18(5):1435–60, October 2002. [50] A. J. R. Gunawardana, “The information geometry of EM variants for speech and image processing,” PhD thesis, Johns Hopkins University, Baltimore, MD., 2001. [51] I. K. Hong, S. T. Chung, H. K. Kim, Y . B. Kim, Y . D. Son, and Z. H. Cho, “Ultra fast symmetry and SIMD-based projection-backprojection(SSP) algorithm for 3- D PET image reconstruction,” IEEE Trans. Med. Imag., 26:789–803, 2007. [52] R. L. Harrison, S. B. Gillipspie, A. M. Alessio, P. E. Kinahan, and T. K. Lewellen, “The effects of object size, attenuation, scatter and random conincidences on sig- nal to noise ratio in simulations of time-of-flight positron emission tomography,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 1900–4, 2005. [53] T. Hebert and R. Leahy, “A Bayesian reconstruction algorithm for emission tomography using a Markov random field prior,” In Proc. SPIE 1092, Med. Im. III: Im. Proc., pages 458–4662, 1989. [54] H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imag., 13(4):601–609, December 1994. [55] I. Hsiao, A. Rangarajan, and G. Gindi, “A new convergent MAP reconstruction algorithm for emission tomography using ordered subsets and separable surro- gates,” In Proc. IEEE Intl. Symp. Biomedical Imaging, pages 409–412, 2002. [56] I. T. Hsiao, A. Rangarajan, and G. Gindi, “A provably convergent OS-EM like reconstruction algorithm for emission tomography,” In Proc. SPIE 4684, Medical Imaging 2002: Image Proc., pages 10–19, 2002. [57] B. W. Jakoby, Y . Bercier, C. C. Watson, V . Rappoport, J. Young, B. Bendriem, and D. W. Townsend, “Physical performance and clinical workflow of a new LSO HI-REZ PET/CT scanner,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 3130–4, 2006. 116 [58] J. I. Jackson, C. H. Meyer, D. G. Nishimura, and A. Macovski, “Selection of a convolution function for Fourier inversion using gridding,” IEEE Trans. Med. Imag., 10(3):473–8, September 1991. [59] C. A. Johnson, J. Seidel, R. E. Carson, W. R. Gandler, M. V . Green, M. E. Daube- Witherspoon, and A. Sofer, “Evaluation of 3D reconstruction algorithms for a small animal PET camera,” IEEE Trans. Nucl. Sci., 44(3):1303–8, 1997. [60] D. J. Kadrmas, “Rotate-and-slant projector for fast LOR-based fully-3-D iterative PET reconstruction,” IEEE Trans. Med. Imag., 27:1071–83, 2008. [61] P. Khurd, I. T. Hsiao, A. Rangarajan, and G. Gindi, “A globally convergent regularized ordered-subset EM algorithm for list-mode reconstruction,” IEEE Trans. Nucl. Sci., 51(3):719–725, June 2004. [62] J. S. Karp, G. Muehllehner, and R. M. Lewitt, “Constrained Fourier space method for compensation of missing data in emission computed tomography,” IEEE Trans. Med. Imag., 7(1):21–5, March 1988. [63] H. Kudo, H. Nakazawa, and T. Saito, “Convergent block-iterative method for general convex cost functions,” In Proc. of the 1999 Int. Mtg. Fully 3D Im. Recon. in Rad. Nuc. Med., pages 247–250, 1999. [64] H. Kudo, H. Nakazawa, and T. Saito, “Block-gradient method for image recon- struction in emission tomography,” Trans. IEICE, J83-D-II(1):63–73, January 2000. In Japanese. [65] P. E. Kinahan and J. G. Rogers, “Analytic 3D image reconstruction using all detected events,” IEEE Trans. Nucl. Sci., 36(1):964–8, February 1989. [66] A. Kuhn, S. Surti, J.S. Karp, P.S. Raby, K.S. Shah, A.E. Perkins, and G. Muehllehner, “Design of a lanthanum bromide detector for time-of-flight PET,” IEEE Trans. Nucl. Sci., 51:2550–7, 2004. [67] Q. Li, E. Asma, J. Qi, and R. M. Leahy, “Accurate estimation of the Fisher information matrix for the PET image reconstruction problem,” IEEE Trans. Med. Imag., 23:1057–64, 2004. [68] K. Lange and R. Carson, “EM reconstruction algorithms for emission and trans- mission tomography,” J. Comput. Assist. Tomogr., 8(2):306–316, April 1984. [69] X. Liu, C. Comtat, C. Michel, P. Kinahan, M. Defrise, and D. Townsend, “Com- parison of 3D reconstruction with OSEM and FORE+OSEM for PET,” IEEE Trans. Med. Imag., 20(8):804–13, August 2001. 117 [70] X. Liu, M. Defrise, C. Michel, M. Sibomana, C. Comtat, P. Kinahan, and D. Townsend, “Exact rebinning methods for three-dimensional PET,” IEEE Trans. Med. Imag., 18(8):657–664, August 1999. [71] R. M. Lewitt, “Reconstruction algorithms: transform methods,” Proc. IEEE, 71(3):390–408, March 1983. [72] T. K. Lewellen, “Recent development in PET detector technology,” Phys. Med. Biol., 53:R287–R317, 2008. [73] Z. Liang, “Detector response restoration in image reconstruction of high resolu- tion positron emission tomography,” IEEE Trans. Med. Imag., 13:314–21, 1994. [74] K. Lee, P. E. Kinahan, J. A. Fessler, R. S. Miyaoka, M. Janes, , and T. K. Lewellen, “Pragmatic fully 3D image reconstruction for the MiCES mouse imag- ing PET scanner,” Phys. Med. Biol., 49:4563–78, 2004. [75] C. Lartizien, P. E. Kinahan, R. Swensson, C. Comtat, M. Lin, V . Villemagne, and R. Trebossen, “Evaluating image reconstruction methods for tumor detec- tion in 3-Dimensional whole-body PET oncology imaging,” Journal of Nuclear Medicine, 44(2):276–90, 2003. [76] R. M. Leahy and J. Qi, “Statistical approaches in quantitative positron emission tomography,” Statistics and Computing, 10(2):147–165, April 2000. [77] W. W. Moses and S. E. Derenzo, “Prospects for time-of-flight PET using LSO scintillator,” IEEE Trans. Nucl. Sci., 46:474–8, 1999. [78] C. L. Melcher, “Scintillation crystals for PET,” J. Nuc. Med., 41(6):1051–5, June 2000. [79] N. A. Mullani, D. C. Ficke, R. Hartz, J. Markham, and G. Wong, “System design of fast PET scanner utilizing time-of-flight,” IEEE Trans. Nucl. Sci., 28:104–7, 1981. [80] S. Matej, J. A. Fessler, and I. G. Kazantsev, “Iterative tomographic image recon- struction using Fourier-based forward and back-projectors,” IEEE Trans. Med. Imag., 23:401–12, 2004. [81] A. Mallon and P. Grangeat, “Three-dimensional PET reconstruction with time- of-flight measurement,” Phys. Med. Biol., 37:717–29, 1992. [82] S. Matej, S. Jayanthi, S. Surti, M. E. Daube-Witherspoon, G. Muehllenhner, and J. S. Karp, “Efficient 3D TOF PET reconstruction using view-grouped histo- images: DIRECT - Direct Image Reconstruction for TOF,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 1728–35, 2006. 118 [83] S. Matej and I. G. Kazantsev, “Fourier-based reconstruction for fully 3-D PET: Optimization of interpolation parameters,” IEEE Trans. Med. Imag., 25:845–54, 2006. [84] J. Muehllehner and J. S. Karp, “Positron emission tomography,” Phys. Med. Biol., 51:R117–37, 2006. [85] M. Moszynski, M. Kapusta, A. Nassalski, T. Szczesniak, D. Wolski, L. Eriksson, and C. L. Melcher, “New prospects for time-of-flight PET with LSO scintilla- tors,” IEEE Trans. Nucl. Sci., 53:2484–8, 2006. [86] S. Matej and R. M. Lewitt, “3-FRP: direct Fourier reconstruction with Fourier reprojection for fully 3-D PET,” IEEE Trans. Nucl. Sci., 48(4-2):1378–1385, August 2001. [87] E. ¨ U. Mumcuoˇ glu, R. M. Leahy, and S. R. Cherry, “Bayesian reconstruction of PET images: methodology and performance analysis,” Phys. Med. Biol., 41(9):1777–1807, September 1996. [88] E. Mumcuoglu, R. Leahy, S. Cherry, and E. Hoffman, “Accurate geometric and physical response modeling for statistical image reconstruction in high resolution PET scanners,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., volume 3, pages 1569–73, 1996. [89] W. W. Moses, “Time of flight in PET revisited,” IEEE Trans. Nucl. Sci., 50:1325– 30, 2003. [90] C. L. Melcher and J. S. Schweitzer, “Cerium-doped lutetium oxyorthosilicate: a fast, efficient new scintillator,” IEEE Trans. Nucl. Sci., 39:502–5, 1992. [91] R. A. Mintzer and S. B. Siegel, “Design and performance of a new pixelated- LSO/PSPMT gamma-ray detector for high resolution PET imaging,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 3418–22, 2007. [92] R. M. Manjeshwar, Y . Shao, and F. P. Jansen, “Image quality improvements with time-of-flight positron emission tomography for molecular imaging,” In Proc. IEEE International Conference on Acoustics, Speech, and Signal Process- ing, volume 5, pages 853–6, 2005. [93] B. Mazoyer, R. Trebossen, C. Schoukroun, B. Verrey, and A. Syrota, “Physical characteristics of TTV03, a new high spatial resolution time-of-flight positron tomograph,” IEEE Trans. Nucl. Sci., 37(2):778–82, 1990. [94] N. Mullani, W. Wong, R. Hartz, E. Philippe, and K. Yerian, “Sensitivity improve- ment of TOFPET by the utilization of the inter-slice coincidence,” IEEE Trans. Nucl. Sci., 29(1):479–83, October 1982. 119 [95] F. Natterer, “Fourier reconstruction in tomography,” Numerische Mathematik, 47:343–53, 1985. [96] F. Natterer, “The mathematics of computerized tomography,” Teubner-Wiley, Stuttgart, 1986. [97] A. Nedi´ c and D. Bertsekas, “Convergence rate of incremental subgradient algo- rithms,” In S. Uryasev and P. M. Pardalos, editors, Stochastic Optimization: Algo- rithms and Applications, pages 263–304. Kluwer, New York, 2000. [98] A. Nedi´ c and D. P. Bertsekas, “Incremental subgradient methods for nondiffer- entiable optimization,” SIAM J. Optim., 12(1):109–138, 2001. [99] R. Neal and G. E. Hinton, “A view of the EM algorithm that justifies incremental, sparse and other variants,” In M. I. Jordan, editor, Learning in Graphical Models, pages 255–268. Kluwer, Dordrencht, 1998. [100] J. M. Ollinger and J. A. Fessler, “Positron emission tomography,” IEEE Signal Processing Mag., 14(1):43–55, January 1997. [101] C. M. Pepin, P. Berard, A. Perrot, C. Pepin, D. Houde, R. Lecomte, C. L. Melcher, and H. Dautet, “Properties of LYSO and recent LSO scintillators for phoswich PET detectors,” IEEE Trans. Nucl. Sci., 51(3):789–95, 2004. [102] G. Pratx, G. Chinn, F. Habte, P. Olcott, and C. Levin, “Fully 3-D list-mode OSEM accelerated by graphics processing units,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 2196–202, 2006. [103] M. E. Phelps, S. C. Huang, E. J. Hoffman, C. Selin, L. Sokoloff, and D. E. Kuhl, “Tomographic measurement of local cerebral glucose metabolic rate in humans with (F-18)2-Fluoro-2-Deoxy-D-Glucose: Validation of method,” Ann. Neurol., 6:371–88, 1979. [104] V . Y . Panin, F. Kehren, C. Michel, and M. Casey, “Fully 3-D PET reconstruction with system matrix derived from point source measurements,” IEEE Trans. Med. Imag., 25:907–21, 2006. [105] J. Qi and R. H. Huesman, “Theoretical study of lesion detectability of MAP reconstruction using computer observers,” IEEE Trans. Med. Imag., 20(8):815– 22, August 2001. [106] J. Qi and R. M. Leahy, “Resolution and noise properties MAP reconstruction for fully 3D PET,” IEEE Trans. Med. Imag., 19(5):493–506, May 2000. [107] J. Qi and R. M. Leahy, “Iterative reconstruction techniques in emission computed tomography,” Phys. Med. Biol., 51:R541–78, 2006. 120 [108] J. Qi, R. M. Leahy, S. R. Cherry, A. Chatziioannou, and T. H. Farquhar, “High resolution 3D Bayesian image reconstruction using the microPET small-animal scanner,” Phys. Med. Biol., 43(4):1001–14, April 1998. [109] J. Qi, R. M. Leahy, C. Hsu, T. H. Farquhar, and S. R. Cherry, “Fully 3D Bayesian image reconstruction for the ECAT EXACT HR+,” IEEE Trans. Nucl. Sci., 45(3):1096–1103, June 1998. [110] A. Ruangma, B. Bai, J. S. Lewis, X. K. Sun, M. J. Welch, R. M. Leahy, and R. Laforest, “Three-dimensional maximum a posteriori (MAP) imaging with radiopharmaceuticals labeled with three Cu radionuclides,” Nuclear Medicine and Biology, 33:217–26, 2006. [111] A. J. Reader, P. J. Julyan, H. Williams, D. L. Hastings, and J. Zweit, “EM algo- rithm system modeling by image-space techniques for PET reconstruction,” IEEE Trans. Nucl. Sci., 50(5):1392–7, October 2003. [112] P. A. Rattey and A. G. Lindgren, “Sampling the 2-D Radon transform,” IEEE Tr. Acoust. Sp. Sig. Proc., 29(4):994–1002, October 1981. [113] A. Ruangma, R. Laforest, B. Bai, and R. M. Leahy, “Characterization of USC- MAP image reconstruction on MicroPET-R4,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 3449–53, 2004. [114] M. Schwaiger, “Myocardial perfusion imaging with PET,” J. Nucl. Med., 35:693– 8, 1994. [115] W. P. Segars, “Development and application of the new dynamic nurbs-based cardiac-torso (NCAT) phantom,” PhD thesis, Univ. of North Carolina, Chpel Hill, NC., May 2001. [116] R. L. Siddon, “Fast calculation of the exact radiological path for a three- dimensional CT array,” Med. Phys., 12(2):252–255, March 1985. [117] S. Surti, J. S. Karp, L. M. Popescu, M. E. Daube-Witherspoon, and M. Werner, “Invesigation of time-of-flight benefit for fully 3-D PET,” IEEE Trans. Med. Imag., 25:529–38, 2006. [118] S. Surti, A. Kuhn, M. E. Werner, A. E. Perkins, J. Kolthammer, and J. S. Karp, “Performance of Philips Gemini TF PET/CT scanner with special consideration for its time-of-flight imaging capabilities,” J. Nuc. Med., 48:471–80, 2007. [119] L. A. Shepp and B. F. Logan, “The Fourier reconstruction of a head section,” IEEE Trans. Nucl. Sci., 21(3):21–43, June 1974. 121 [120] V . V . Selivanov, Y . Picard, J. Cadorette, S. Rodrigue, and R. Lecomte, “Detector response models for statistical iterative image reconstruction in high resolution PET,” IEEE Trans. Nucl. Sci., 47(3):1168–75, June 2000. [121] H. R. Schelbert, M. E. Phelps, S. Huang, N. S. MacDonald, H. Hansen, C. Selin, and D. E. Kuhl, “N-13 ammonia as an indicator of myocardial blood-flow,” Cir- culation, 63:1259–72, 1981. [122] D. W. Shattuck, J. Rapela, E. Asma, A. Chatzioannou, J. Qi, and R.M. Leahy, “Internet2-based 3D PET image reconstruction using a PC cluster,” Phys. Med. Biol., 47:2785–95, 2002. [123] D. Strul, R. B. Slates, M. Dahlbom, S. R. Cherry, and P. K. Marsden, “An improved analytical detector response function model for multilayer small- diameter PET scanners,” Phys. Med. Biol., 48:979–94, 2003. [124] D. L. Snyder, L. J. Thomas, and M. M. Ter-Pogossian, “A mathematical model for positron-emission tomography systems having time-of-flight measurements,” IEEE Trans. Nucl. Sci., 28(3):3575–3583, June 1981. [125] L. A. Shepp and Y . Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imag., 1(2):113–122, October 1982. [126] L. A. Shepp, Y . Vardi, J. B. Ra, S. K. Hilal, and Z. H. Cho, “Maximum likelihood PET with real data,” IEEE Trans. Nucl. Sci., 31(2):910–913, April 1984. [127] T. Tomitani, “Image reconstruction and noise evaluation in photon time-of-flight- assisted positron emission tomography,” IEEE Trans. Nucl. Sci., 28:4582–9, 1981. [128] M. M. Ter-Pogossian, D. C. Ficke, M. Yamamoto, and J. T. Hood, “Super PETT I: A positron emission tomograph utilizing photon time-of-flight information,” IEEE Trans. Med. Imag., 1(3):179–187, November 1982. [129] Y . Tai, A. Ruangma, D. Rowland, S. Siegel, D. F. Newport, and P. L. Chow, “Performance evaluation of the microPET Focus: A third-generation microPET scanner dedicated to animal imaging,” Journal of Nuclear Medicine, 46:455–63, 2005. [130] S. Vandenberghe, M. E. Daube-Witherspoon, R. M. Lewitt, and J. S. Karp, “Fast reconstruction of 3D time-of-flight PET data by axial rebinning and transverse mashing,” Phys. Med. Biol., 51:1603–21, 2006. [131] C. W. E. van Eijk, “Inorganic scintillators in medical imaging,” Phys. Med. Biol., 47:R85–R106, 2002. 122 [132] K. Vunckx, L. Zhou, S. Matej, M. Defrise, and J. Nuyts, “Fisher information- based evaluation of image quality for time-of-flight PET,” In Proc. IEEE Nuclear Science Symp. Medical Imaging Conf., pages 4129–36, 2007. [133] C. C. Watson, “An evaluation of image noise variance for time-of-flight PET,” IEEE Trans. Nucl. Sci., 54:1939–47, 2007. [134] H. N. Wagner, H. D. Burns, R. F. Dannals, D. F. Wong, B. Langstrom, T. Duelfer, J. J. Frost, H. T. Ravert, J. M. Links, S. B. Rosenbloom, S. E. Lukas, A. V . Kramer, and M. J. Kuhar, “Imaging domamine receptors in the human brain by positron tomography,” Science, 221:1264–6, 1983. [135] K. Wienhard, M. Schmand, M. E. Casey, K. Baker, J. Bao, L. Eriksson, W. F. Jones, C. Knoess, M. Lenox, M. Lercher, P. Luk, C. Michel, J. H. Reed, N. Richerzhagen, J. Treffert, S. V ollmar, J. W. Young, W. D. Heiss, and R. Nutt, “The ECAT HRRT: performance and first clinical application of the new high resolution research tomograph,” IEEE Trans. Nucl. Sci., 49(1):104–10, 2002. [136] M. Yavuz and J. A. Fessler, “Statistical image reconstruction methods for randoms-precorrected PET scans,” Med. Imag. Anal., 2(4):369–378, December 1998. [137] T. Yamaya, N. Hagiwara, T. Obi, B. Yamaguchi, N. Ohymama, K. Kitamura, T. Haseqawa, H. Haneishi, E. Yoshida, N. Inadama, and H. Murayama, “Transax- ial system models for jPET-D4 image reconstruction,” Phys. Med. Biol., 50:5339– 55, 2005. 123

Abstract (if available)

Abstract Positron emission tomography (PET) is a functional biomedical imaging technique that provides in vivo information about physiological processes within the body by reconstructing the 3D distribution of a adiolabelled tracer. It is challenging to estimate the tracer distribution accurately since photon-limited PET data have high statistical variance. Statistical iterative reconstruction methods have shown uperior image quality and better quantitative results compared to analytical reconstruction methods.However, the iterative methods involve huge computational loads, particularly for fully 3D PET, consequently limiting their routine use in practice. Clinical 3D time-of-flight (TOF) PET scanners, which become available due to recent developments of fast scintillators, involve an even higher computational cost for image reconstruction because of extra TOF information. A successful pragmatic approach that reduces computation cost for fully 3D PET is to use Fourier rebinning (FORE), reducing 3D PET data to 2D data, combined with a 2D reconstruction method such as the ordered subsets expectation-maximization (OSEM) algorithm. The combination of FORE and 2D OSEM is now routinely used in both clinical and small animal imaging systems. However, data dimensionality reduction using FORE can result in a loss in resolution or noise performance. In the first part of this dissertation, we develop fast projectors that map a 3D image into 3D PET data by using the exact inverse (Fourier) rebinning operator. The inverse rebinning operator is used for calculating 3D PET data from 2D PET data, and combined with a 2D projector to construct a fast 3D projection operator. The inverse rebinning operator is implemented cost-effectively using the fast Fourier transform (FFT), so that the computation time is reduced by an order of magnitude. We utilize this fast projector in the context of fully 3D PET maximum a posteriori (MAP) reconstruction methods.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

3D deep learning for perception and modeling

PDF

Point-based representations for 3D perception and reconstruction

PDF

3D inference and registration with application to retinal and facial image analysis

PDF

Model based PET image reconstruction and kinetic parameter estimation

PDF

3D face surface and texture synthesis from 2D landmarks of a single face sketch

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Reconstructing 3D reconstruction: a graphical taxonomy of current techniques

PDF

Machine learning methods for 2D/3D shape retrieval and classification

PDF

Feature-preserving simplification and sketch-based creation of 3D models

PDF

Face recognition and 3D face modeling from images in the wild

PDF

Labeling cost reduction techniques for deep learning: methodologies and applications

PDF

Data-driven 3D hair digitization

$Metasurfaces in 3D applications: multiscale stereolithography and inverse design of diffractive optical elements for structured light$

PDF

Metasurfaces in 3D applications: multiscale stereolithography and inverse design of diffractive optical elements for structured light

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

3D vessel mapping techniques for retina and brain as an early imaging biomarker for small vessel diseases

PDF

A data driven software platform for process automation, planning and inspection of Contour Crafting large-scale robotic 3D printing system

Asset Metadata

Creator Cho, Sang Hee (author)

Core Title Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 10/05/2009

Defense Date 07/31/2008

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Fourier rebinning,fully 3D PET,inverse rebinning,iterative reconstruction,OAI-PMH Harvest,time-of-flight PET

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Leahy, Richard M. (committee chair), Chatziioannou, Arion (committee member), Kuo, C.-C. Jay (committee member), Singh, Manbir (committee member)

Creator Email sanghcho@sipi.usc.edu,scho.2000@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1623

Unique identifier UC164630

Identifier etd-CHO-2369 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-117793 (legacy record id),usctheses-m1623 (legacy record id)

Legacy Identifier etd-CHO-2369.pdf

Dmrecord 117793

Document Type Dissertation

Rights Cho, Sang Hee

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu