Statistical lesion detection in dynamic positron emission tomography

PDF

Download

Open document

Flip pages

Summarize

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content STATISTICAL LESION DETECTION IN DYNAMIC POSITRON EMISSION TOMOGRAPHY by Zheng Li 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) May 2008 Copyright 2008 Zheng Li Dedication To my parents and my wife, whose endless encouragement made this possible. ii Acknowledgements I would like to express my deepest gratitude to my advisor, Prof. Richard M. Leahy. I have benefited greatly from your rich experience in medical imaging and signal pro- cessing. Thank you for your unwavering support during my Ph.D study, for guidance in solving research problems, for providing a friendly research environment, and for providing financial support. It has really been a pleasure working with you. I would like to thank Prof. Xiaoli Yu, who kindly guided and supported my study. I benefited a lot from her extensive knowledge in feature extraction and signal detec- tion. I also would like to thank my dissertation committee members: Prof. Shrikanth Narayanan, Prof. Remigijus Mikulevicius and Prof. Antonio Ortega for their valuable suggestions and feedback. I would like to thank Prof. Peter Conti at the USC school of medicine, for providing experiment facilities and data; as well as Prof. Jinqi Qi at UC Davis for providing assistance with image variance estimation. I would like to thank my colleague Dr. Quanzheng Li for his valuable assistance in PET reconstruction and system simulation, Dr. Dimitrios Pantazis for helping me build up my knowledge about random field theory, and Dr. Sangtae Ahn for enlightening research discussions. In the University of Southern California I have enjoyed a friendly environment and an excellent research lab. I would like to thank its past and present iii members: Dr. Belma Dogdas, Dr. Abhijit Chaudhari, Sanghee Cho, Anand Joshi, Sangeetha Somayajula, Joyita Dutta, Hua Hui and Yu-Teng Chang. Last but not the least, I would like to express my deepest thanks to my wife Nan. This work could not have been done without her enduring support and encouragement. This work was supported by the National Institute of Biomedical Imaging and Bio- engineering under Grant No. R01 EEB00363, and the BCRP Idea Award under Grant No. DAMD17-99-1-9379 iv Table of Contents Dedication ii Acknowledgements iii List of Tables vii List of Figures viii Abstract xii Chapter 1: Introduction 1 1.1 Positron Emission Tomography . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation and Organization . . . . . . . . . . . . . . . . . . . . . . . 4 Chapter 2: Dynamic FDG-PET TAC 7 2.1 Tracer Kinetics Modeling for Dynamic FDG-PET . . . . . . . . . . . . 8 2.1.1 A Homogeneous V oxel . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 A Heterogeneous V oxel . . . . . . . . . . . . . . . . . . . . . 10 2.2 Time-Frequency Properties of TACs . . . . . . . . . . . . . . . . . . . 11 2.2.1 Physiological Processf(¯;t) vs. Macro Parameter . . . . . . . 12 2.2.2 Frequency Domain Properties of TACs . . . . . . . . . . . . . 13 2.2.3 Time-Frequency Tiling of Tumor and Normal TACs . . . . . . 14 2.3 Subspace Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1 The Formation of Subspaces . . . . . . . . . . . . . . . . . . . 17 2.3.2 Subspace Identification Algorithms . . . . . . . . . . . . . . . 18 2.3.3 Subspace Identification Results . . . . . . . . . . . . . . . . . 20 Chapter 3: The Matched Subspace Method 23 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.1 Formulation of the Subspace Detection Problem . . . . . . . . 28 3.2.2 Generalized Likelihood Ratio (GLR) Test . . . . . . . . . . . . 29 3.2.3 Subspace Identification . . . . . . . . . . . . . . . . . . . . . . 35 v 3.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.1 Implementation of the Matched Subspace Detectors . . . . . . . 36 3.3.2 Noise Covariance Estimation . . . . . . . . . . . . . . . . . . . 39 3.3.3 Receiver Operation Characteristic . . . . . . . . . . . . . . . . 40 3.4 Simulations, Experiments and Results . . . . . . . . . . . . . . . . . . 42 3.4.1 Digital Phantom Simulation . . . . . . . . . . . . . . . . . . . 42 3.4.2 Physical Phantom Study . . . . . . . . . . . . . . . . . . . . . 49 3.4.3 A Clinical Example . . . . . . . . . . . . . . . . . . . . . . . . 51 3.5 Analysis and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.5.1 Detectability vs. the Number of Prompts . . . . . . . . . . . . 54 3.5.2 False Alarms and Image Artifact . . . . . . . . . . . . . . . . . 57 3.5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Chapter 4: Controlling Familywise Error Rate in Thresholding 60 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2.1 Random Field Theory . . . . . . . . . . . . . . . . . . . . . . 64 4.2.2 Controlling FWER forL 2 (y) Statistic Maps . . . . . . . . . . . 68 4.2.3 The Matched Subspace Detection and Gaussian Approximation 69 4.2.4 Normalization of the Non-stationary Gaussian Random Field . . 75 4.3 Simulations, Experiments and Results . . . . . . . . . . . . . . . . . . 76 4.3.1 Digital Phantom Simulation . . . . . . . . . . . . . . . . . . . 76 4.3.2 A Clinical Experiment . . . . . . . . . . . . . . . . . . . . . . 80 4.4 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 81 Chapter 5: Conclusion and Future Work 82 5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2.1 Lesion Detection Experiment . . . . . . . . . . . . . . . . . . 84 5.2.2 Detection Algorithm . . . . . . . . . . . . . . . . . . . . . . . 85 References 87 Appendix: Physiological process function 95 vi List of Tables 2.1 MSE for the TAC curve fitting. The left table shows the spectral analysis method. The right table shows the PCA method. . . . . . . . . . . . . . 22 3.1 Micro Parameters for True TACs Used in Digital Phantoms . . . . . . . 45 3.2 Monte Carlo Simulation DATA Acquisition Protocol . . . . . . . . . . 45 3.3 P D , P FA and AUC Values for ROC Study. P D is the probability of detection. P FA is the probability of false positive. . . . . . . . . . . . . 48 3.4 P D vs. P FA for the Different Number of Prompt Settings . . . . . . . . 58 4.1 Micro Parameters for True TACs Used in the Test Statistic Gaussianity Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 vii List of Figures 1.1 Basic principles of PET imaging. (a)»(c): Physical principles. (d)»(f): Image reconstruction. The procedure of tracer production using a cyclotron and the application of the tracer to the patient are not included in this fig- ure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 The three compartment 4-k model. . . . . . . . . . . . . . . . . . . . . 8 2.2 Clinical PET image and region of interest (ROI) TACs. Left: The last frame of a dynamic PET MAP reconstructed image. Right: TACs drawn from two lesion and two normal tissue ROIs. . . . . . . . . . . . . . . 11 2.3 The physiological process functionf(¯;t) vs. ¯. . . . . . . . . . . . . 13 2.4 A two-stage wavelets analysis filter bank. . . . . . . . . . . . . . . . . 14 2.5 The Daubechies wavelet basis functions for a 512-point, seven-level fil- ter bank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6 Wavelet tilings of tumor and normal tissue TACs. In the left column: clinical noisy and de-noised TACs. In the right column: time-frequency tilings. The upper row is for normal tissue. The lower row is for tumor tissue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.7 A Comparison of the accuracy and separability of the subspaces hHi andhSi estimated using the spectral analysis method. . . . . . . . . . . 21 2.8 A Comparison of the accuracy and separability of the subspaces hHi andhSi estimated using the PCA method. . . . . . . . . . . . . . . . . 21 3.1 Illustration for the matched subspace method. y is the observed signal, y 1 is the projection of y in the one-dimensional subspacehe 1 i which produces large mismatch. Including the basis functionse 2 ande 3 in the subspace will diminish the mismatch. . . . . . . . . . . . . . . . . . . 24 viii 3.2 Illustration of spatial pre-whitening. The5£5 dash line grid represents the 25-voxel region applying spatial pre-whitening for the center voxel. We plot two5£5 pre-whitening regions and corresponding center voxels to illustrate that the pre-whitening region moves with the center voxel. . 37 3.3 TACs used in the simulations. Left: The TACs of the blood input func- tion and the primary tumor in Phantom 1. Right: The TACs of normal tissue and metastasis in Phantom 2. All phantoms use the same true TAC for all normal tissue. The noisy TACs shown are for a 3£3 ROI drawn from a single reconstruction for the data used in the ROC study. . 43 3.4 Digital phantom images (for all dynamic images, only the last frame is shown): (1) Digital phantom 1, which contains normal tissue, one primary tumor and blood pool; (2) Digital phantom 2, which contains normal tissue and six metastases; Phantom 3 contains pure background and is not shown here. (a) The noise free digital phantom; (b) A recon- structed dynamic image; (c) A reconstructed static image. The frame duration was chosen such that this image has the maximum SNR achiev- able (in terms of the Hotelling observer) by any single frame static image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 Simulation results for Hotelling observer and Patlak method. (a) Static image Hotelling observer on 9£9 window. (b) Patlak method applied to3£3 averaged TAC. The upper row shows 2D test statistic images for the Digital Phantom 2 (Figure 3.4 (2.a)). The lower row shows profiles through these images at the level of the two sets of three metastases. . . 47 3.6 Simulation results for matched subspace detection methods. (c)L 1 (y): matched subspace detector with unknown variance white noise data model. The detector is applied to3£3 averaged TAC. (d)L 2 (y): matched sub- space detector with known noise covariance. The detector is applied to 3£3 averaged TAC. (e) L 2 (y ROI ): matched subspace detector with known noise covariance. The detector is applied to the 9 TACs within a 3£3 window. The upper row shows 2D test statistic images for the Dig- ital Phantom 2 (Figure 3.4 (2.a)). The lower row shows profiles through these images at the level of the two sets of three metastases. . . . . . . 48 3.7 Lesion detection ROC curves computed from Monte Carlo simulations. The method names shown in the legend are ranked, from high to low, based on AUC values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 ix 3.8 Noiseless TACs in the physical phantom. The decay coefficient for the background TAC, represented by the blue curve, is¯ 1 = 8:178£10 ¡4 . The decay coefficient for the lesion TAC, represented by the red curve, is¯ 2 =1:518£10 ¡4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.9 Physical phantom study. Top left: The plane containing the large tumor (Lesion 1). Top right: The plane containing the small tumor (Lesion 3). Bottom: L 1 result for Lesion 3 plane. . . . . . . . . . . . . . . . . . . . 51 3.10 Clinical breast cancer study. 1st row: A transverse plane containing the primary tumor. The ROIs used for subspace identification are shown on the left most image. 2nd row: Transverse plane containing metastatic lesion. 3rd row: Coronal plane containing the metastasis. 4th row: Sagittal plane containing the metastasis. . . . . . . . . . . . . . . . . . 53 3.11 Clinical breast cancer metastasis early detection. (A) The L 2 detection results based on the 1st scan data, the follow-up scan images and the fused images. (B) The first scan images, the follow-up scan images and the fused images. The red arrow indicates the position of the detected metastasis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.12 Lesion detection ROC curves change with the number of prompts. The method names shown in the legend are ranked, from high to low, based on AUC values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.13 Lesion detection ROC curves for the different number of prompt set- tings. The method names shown in the legend are ranked, from high to low, based on AUC values. . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 The block diagram for thresholding an inhomogeneous L 2 (y) statistic map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 The histogram and fitted Gaussian PDF of a single voxel L 2 (y) test statistic for both normal and tumor tissue. The perfect match between the histogram and the fitted Gaussian PDF implies an accurate Gaussian approximation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3 Quantile-quantile plots of the single voxel L 2 (y) test statistic versus a standard Gaussian random variable. The linearity of the quantile- quantile plot validates the Gaussian approximation for the L 2 (y) test statistic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 x 4.4 The profiles of a phantom image before and after removing the ensemble mean in each segmentation region. The upper plot shows the profile before removing the ensemble mean. The lower plot shows the profile after removing the ensemble mean. A partial volume effect causes peaks on the boundaries in the mean-value-removed image. These peaks can induce an overestimated ensemble variance. . . . . . . . . . . . . . . . 74 4.5 Histograms of the spatially sampled L 2 (y) test statistic values in three segmentation regions (Figure 4.7(b)) before (Left column), and after (Right column) removing the outliers. The fitted curves are the Gaussian PDF with the ensemble mean and variance of each region. . . . . . . . 75 4.6 Clinical images, true phantoms and reconstructed images. On the 1st row, a clinical tumor free image (left), a clinical image with one primary tumor and one metastasis (right). On the 2nd row, true phantom images. On the 3rd row, reconstructed images. . . . . . . . . . . . . . . . . . . 77 4.7 Simulation results for the thresholding algorithm. (a) TheL 2 (y) statistic map, (b) Three segmentation regions, (c) The normalizedN(0;1) GRF, (d) The binary decision map. . . . . . . . . . . . . . . . . . . . . . . . 78 4.8 Simulated FWER against preset FWER. The slope is less than one, which implies a conservative thresholding. . . . . . . . . . . . . . . . . 79 4.9 A clinical breast cancer study. The 1st column shows the transverse, sagittal and coronal views of the last dynamic frame. The 2nd»5th columns show thresholding results in the primary tumor, metastasis 1, metastasis 2 and a false alarm, respectively. The 1st row shows the original image. The 2nd row shows the normalized N(0,1) Gaussian random field. The 3rd row shows the binary decision map. The results show that the thresholding algorithm detects all three tumors but also generates one false alarm out of 31 planes. . . . . . . . . . . . . . . . . 79 xi Abstract Positron emission tomography (PET) with 18 F fluorodeoxyglucose (FDG) is improving the diagnosis, staging and treatment monitoring of a variety of cancers in comparison to alternative noninvasive imaging modalities [Con96, SLDV04]. However, detection of small lesions is limited by image resolution and low signal to noise ratio. While computer-aided detection algorithms based on static PET may assist visual detection of lesions, the temporal domain information loss in the static image could impair the detection of lesions. In addition, spatial features such as the shape and contrast of a lesion in static PET images are highly variable and difficult to characterize [HYB + 97]. Since the dynamics of tracer uptake (tracer radioactivity) are different for normal and malignant tissue, the temporal domain provides potentially important additional features for use in lesion detection. This dissertation focuses on three challenges in the detection of lesions in dynamic PET: feature extraction, detector design and thresholding procedure. A linear subspace model with additive Gaussian noise is used to describe the received time activity curve (TAC) which is a function of local tracer radioactivity versus time. Using TACs from a primary tumor region of interest (ROI) and one or multiple normal ROIs, two linear subspaces can be identified. We use two sets of subspaces instead of shape or intensity information as the features in lesion detection. A matched subspace detection algorithm based on the generalized likelihood ratio test [VO91, Sch90] was proposed to assist xii in the detection of small tumors. The algorithm differentiates tumor tissue from non- tumor background using the TACs that characterize the uptake of PET tracers. Applying a matched subspace detector with the identified subspaces to TACs on a voxel by voxel basis throughout the dynamic image produces a test statistic at each voxel which on thresholding indicates potential locations of secondary or metastatic tumors. Applying the detector to the dynamic image results in a statistic image, or a statistic map. The detector is derived for three cases: detection on a single TAC with unknown variance white noise, detection on a single TAC with known noise covariance, and detection on multiple TACs within a small ROI with known noise covariance. To make the binary decision on the presence or absence of a lesion at each voxel, the continuous-value statistic map is processed by a constant familywise error rate (FWER) thresholding algorithm. The statistic map is modeled as an inhomogeneous Gaussian random field. The thresholding procedure includes three steps. First, the PET image is segmented into several homogeneous regions. Then the statistic map is normalized to a zero mean unit variance Gaussian random field. Finally a maximum statistic threshold- ing with a fixed FWER is applied. To evaluate the proposed matched subspace detector, a receiver operation character- istic (ROC) study for dynamic PET tumor detection was designed. The detector uses a dynamic sequence of frame-by-frame 2-dimensional (2-D) MAP reconstructions as input. We compared the performance of subspace detectors with that of a Hotelling observer [BGG + 92] applied to a single frame image and parametric Patlak [PBF83] method. Our Monte Carlo simulations show that the matched subspace detector is supe- rior to the 2-D static Hotelling observer and the Patlak method for the lesion detection task. The application of detection approaches to clinical PET data was also demon- strated. To evaluate the proposed thresholding algorithm, we used a digital phantom gen- erated from clinical dynamic images. Simulations show that with a preset FWER=5%, xiii the proposed thresholding method achieves a 97.1% true positive rate, with 1.8% FWER. We also demonstrate here the application of the proposed approach to clinical PET data from a breast cancer patient with metastatic disease. xiv Chapter 1 Introduction 1.1 Positron Emission Tomography Positron emission tomography (PET) is a medical imaging modality for producing 3- dimensional (3-D) images of the spatial distribution of biochemical tracers within bio- logical systems [LQ00]. It enables physicians to non-invasively and quantitatively deter- mine metabolic activity by tracing a radioactive isotope [CP96]. PET is routinely used in both clinical and small animal imaging on applications ranging from lesion detection to gene expression. The basic principles and image reconstruction algorithms have been described in many books and review papers [LQ00, CP96, KS87, CJS93]. Physical Principles In a PET study, a radiolabeled biologically active compound (tracer) is injected into or inhaled by a subject. The positron emitting isotope, which is attached to the tracer, then circulates through the blood stream. Positrons are emitted from the nucleus of unstable radioisotopes attached to the tracer, because they have an excessive number of protons and a positive charge. A positron emitted from a decaying nucleus travels a short distance before colliding with the electron of a nearby atom (see Figure 1.1(a) 1 ). The distance that a positron travels ranges from 0.3mm to 2.7mm for clinically used isotopes and tracers, depending on its energy. The annihilated energy is carried by two 511-keV photons that travel in approximately opposite directions (180 ± § 5 ± ). The slight deviation from the opposite directions is due to the residual kinetic energy 1 Crump Institute for Molecular Imaging, UCLA. 1 of the positron. The physical basis for positron emission tomography lies in the fact that the annihilation photons are emitted simultaneously in opposite directions. This pair of photons, after traveling in opposite directions along a straight line through the human body, can be recorded by radiation detectors. In most scanners, the detectors consist of a combination of scintillators and photomultiplier tubes (PMTs). Typically 64 scintillation detectors will be coupled to four PMTs as shown in Figure 1.1(c). Scintillators convert the high energy 511KeV photons into a large number of low energy photons which are collected and amplified by the PMTs. Arranging the detec- tor blocks in a circular fashion produces a ring scanner. Additional rings of blocks are added to increase the axial filed of view (FOV) of the scanner. If both members of the detector pair receive photon signals within a very narrow time interval, a coincidence event will be detected, which allows us to approximately determine the line of response (LOR) containing the positron emitter (see Figure 1.1(b)). After collecting enough data, which usually means that millions of coincidence events are detected, the two or three dimensional source distribution can be reconstructed using certain algorithms. Image Reconstruction Consider the simple 2D arrangement as illustrated in Figure 1.1(b) and (c). In the absence of an attenuating object, the total number of coincidences between any detector pair is approximately a line integral through the 2D source dis- tribution. The set of line integrals of a 2D function forms a Radon transform [SL74], and therefore, the coincidence data is readily sorted into a Radon transform or sino- gram format as illustrated in Figure 1.1(d). Once the sinogram data has been collected, the source distribution can be reconstructed using a conventional filtered backprojec- tion (FBP) algorithm which inverts the Radon transform of a 2D function [SL74]. For the 3D detector arrangement, the inverse Randon transform cannot be applied directly. With the rebinning technique, this problem can be solved by using 2D reconstruction 2 algorithms with little loss in performance. These reconstruction algorithms fall into a category called the analytical approach. Although the analytical approaches are fast, the accuracy of the reconstruction is lim- ited by the implicit approximations in the line integral model on which the reconstruc- tion formulae are based. In contrast, statistical methods adopt arbitrarily accurate mod- els for mapping between the source volume and the sinogram. Statistical approaches also allow the explicit modeling of statistical noise associated with photon detection [LQ00]. Statistical approaches typically achieve much better image quality than analyt- ical methods, but they also greatly increase the computation load. For example, a single 3D scan from a typical PET system could produce 10 7 » 10 8 sinogram elements with 10 6 image elements to be estimated. The well-studied statistical algorithms include Maximum Likelihood (ML) [SV82] and the Maximum a Posterior (MAP) [QLC + 98] algorithms. The Ordered Subset Expectation Maximum (OSEM) algorithm [HL94] is widely used to speed up convergence for ML estimators. Dynamic PET Dynamic PET aims at imaging biological functions over a certain period of time, and providing an activity estimate for each time point during that period. Traditionally, dynamic PET image reconstruction is performed by dividing the data set into multiple frames and reconstructing the individual frames independently. Conven- tional static image reconstruction algorithms, either analytical or statistical approaches, can be employed to reconstruct each frame. The new approach of reconstructing contin- uous dynamic PET image (4D reconstruction: 3D spatial + 1D time) based on listmode data [LAAL07] is currently of interest. Dynamic PET studies are essential to the investigation of basic chemical and func- tional processes in physiology, anatomy, molecular biology and pharmacology. This technique has been used to image blood volume and flow, local cerebral metabolic rates 3 of glucose (LCMRGlc) [HPH + 80, Sul93] and other biological processes. To character- ize the tracer kinetic behaviors, compartmental models [HPH + 80, HCH + 83] are usually employed. 1.2 Motivation and Organization Positron emission tomography (PET) with 2-[fluorine-18]-fluoro-2-deoxy-D-glucose ( 18 F-FDG, commonly abbreviated to FDG) is improving the diagnosis, staging and treat- ment monitoring in a variety of cancers in comparison to alternative noninvasive imag- ing modalities [Con96, SLDV04]. Although conventional static PET imaging is highly sensitive in tumor detection, further improvement is needed, because even a small per- centage of the false negatives result in high risk and cost. Conventional visual inspection and diagnosis is potentially inaccurate due to limited spatial resolution and low lesion- to-background contrast, particularly for small tumors. This motivated research into the development of algorithms for both static and dynamic PET to assist in lesion detec- tion, which can help assess cancer treatment efficacy or detect tumors in early stages. Computer aided detection or diagnosis techniques have been widely studied in mam- mography, CT and MRI [GKA01, Doi07]. But there is considerably less application to PET images due to the relative inaccessibility of PET equipment [DCGJ99]. In static PET, lesion detection algorithms [HYB + 97, DCGJ99, YZS + 04] were proposed to assist in visual inspection. While these algorithms work to some extent, the temporal domain information loss during static image reconstruction could reduce the detectability of lesions. Since the dynamics of tracer uptake are different for normal and malignant tis- sue, the temporal domain provides potentially important additional features which may improve lesion detection performance. The research in this thesis is aimed at developing 4 new lesion detection methods using dynamic FDG-PET. The research is organized into the following sections: Chapter 2 presents the dynamic PET time activity curve (TAC) modeling technique and the concept of using a linear subspace for the TAC. Using the compartment model and wavelet decomposition, we analyze the time-frequency properties of the TAC, which provide evidence for the separability of normal and tumor TAC’s linear subspaces. The formulation of the linear subspace and subspace identification techniques are introduced. In chapter 3, we describe a matched subspace detection algorithm in dynamic PET. The hypothesis testing model is presented, followed by the derivation of the generalized likelihood ratio test (GLRT). We demonstrate our algorithm with digital phantom, phys- ical phantom and clinical data studies. A receiver operation characteristic (ROC) study is conducted to evaluate and compare the performance of the proposed detector to the existing lesion detection methods. Chapter 4 introduces the use of random fields to threshold the test statistical map. We first validate the Gaussian approximation for the test statistic produced by the matched subspace method. With the given image segmentation, the statistical map is normalized to a homogenous zero mean unit variance Gaussian random field. By using the random field theory, the number of resolution elements and the threshold value for a given fam- ilywise error rate are estimated. The algorithm is demonstrated using digital phantom and clinical data. We conclude this work and discuss the possible directions for future research in Chapter 5. 5 (a) Annihilation (b) Coincidence (e) Sinogram (f) Reconstructed Image Reconstruction Algorithms γ γ f(x,y) θ 2 (d) Line Integral (c) Detector θ 3 θ 4 θ 1 Figure 1.1: Basic principles of PET imaging. (a)»(c): Physical principles. (d)»(f): Image reconstruction. The procedure of tracer production using a cyclotron and the application of the tracer to the patient are not included in this figure. 6 Chapter 2 Dynamic FDG-PET TAC Dynamic PET imaging refers to the measurement of tracer uptake over time. An image of kinetic tracer activity distribution is a good starting point for obtaining more use- ful biological information. The kinetic behavior of the tracer concentration in a cer- tain region can be measured by the time activity curve (TAC), which represents the counts/second/voxel (or counts/second/ml) in the given region as a function of time. We will focus on FDG-PET dynamic studies because FDG is the most widely used tracer in clinical as well as small animal PET scans. Using tracer kinetic modeling techniques, biological information such as the local influx constant [PBF83] or distribution volume ratio [Log00], can be estimated from TAC data. The TAC also plays an important role in tumor diagnosis. It has been observed that a tumor’s TAC is different from that of normal tissue [ZW96, ZPL + 01]. Checking the TAC can facilitate tumor diagnosis, par- ticularly of a small size or low contrast tumor where diagnosis based on a static image is difficult. In this chapter, two tracer kinetic modeling methods are introduced. By using these kinetic models, the frequency domain differences between tumor and nor- mal TACs are discussed. The goal is to demonstrate the separability and provide a method to differentiate between tumor and normal TACs based on the linear subspace idea. 7 C p (t) C e (t) C m (t) K 1 k 2 k 3 k 4 Plasma Tissue Figure 2.1: The three compartment 4-k model. 2.1 Tracer Kinetics Modeling for Dynamic FDG-PET 2.1.1 A Homogeneous Voxel Under the homogeneous tissue assumption, the FDG tracer distribution in tissue can be modeled by the three compartments 4-k model (Figure 2.1). The left compartment represents the blood pool region. C p (t) represents the FDG concentration in arterial plasma, which is considered to be a deterministic function. There are several techniques available to obtain C p (t): (i) normalized population-based input function, (ii) invasive arterial line sampling, (iii) ROIs of aorta or large vessels. The middle compartment contains the FDG in the tissue, with the concentration C e (t). The right compartment contains the phosphorylate FDG (FDG-6-P) in the tissue, with the concentrationC m (t). The middle and right compartments together represent the tissue region. The parameters K 1 »k 4 1 are the rate constants defined as: K 1 , the transport rate of FDG in plasma to tissue. k 2 , the transport rate of FDG in tissue to plasma. k 3 , the phosphorylation rate of FDG to FDG-6-P. 1 All tracer concentration functions (C p (t),C e (t) andC m (t)) and rate constant values (K 1 »k 4 ) are decay-corrected quantities. 8 k 4 , the dephosphorylation rate of FDG-6-P to FDG. According to the 3-compartmental 4-k model [HPH + 80], the molecular kinetics of FDG and physiological status in homogeneous tissues are characterized by the following differential equations: dC e (t) dt =K 1 C p (t)¡(k 2 +k 3 )C e (t)+k 4 C m (t) C m (t) dt =k 3 C e (t)¡k 4 C m (t) (2.1) It can be shown [HPH + 80] that the uptake of the tracer for a homogeneous tissue C tis is: C tis (t) =C e (t)+C m (t)=[M 1 e ¡¯ 1 t +M 2 e ¡¯ 2 t ]C p (t) (2.2) where ¯ 1 = 1 2 h (k 2 +k 3 +k 4 )¡ p (k 2 +k 3 +k 4 ) 2 ¡4k 2 k 4 i ¯ 2 = 1 2 h (k 2 +k 3 +k 4 )+ p (k 2 +k 3 +k 4 ) 2 ¡4k 2 k 4 i M 1 = K 1 ¯ 2 ¡¯ 1 (k 3 +k 4 ¡¯ 1 ) M 2 = K 1 ¯ 2 ¡¯ 1 (¯ 2 ¡k 3 ¡k 4 ) (2.3) The total activity measured by the PET scanner (assuming there is no noise), C(t), including the activity in the blood, can be represented by: C(t)=V p C p (t)+C tis (t)=V b C p (t)+[M 1 e ¡¯ 1 t +M 2 e ¡¯ 2 t ]C p (t); (2.4) where V p denotes vascular space within the tissue, which means that the FDG concen- tration in blood is assumed to beV p C p (t) at timet [PMS86]. 9 2.1.2 A Heterogeneous Voxel Due to the limited spatial resolution of the scanner, the tissue concentration in a sin- gle voxel measured by the scanner can be better modeled by the mass-weighted aver- age concentration of radioactivity in many homogeneous tissues [LSM + 93],[SMS91]. Given that a mixed (heterogeneous) voxel contains ¹ J homogeneous subregions, each of which has a weighting coefficient ofw j . The total radioactivity in a heterogeneous voxel measured by a PET scanner, ¹ C(t), can be represented as: ¹ C(t) = V p C p (t)+ ¹ C tis (t) = V p C p (t)+ ¹ J X j=1 w j [M 1j e ¡¯ 1j t +M 2j e ¡¯ 2j t ]C p (t) (2.5) Notice that C p (t) ¼ C p (t)¯e ¡¯t when ¯ is a large number 2 . For clinical data, the above approximation is very accurate when ¯ is close to or larger than 1 sec. ¡1 . With this approximation, the first term in (2.5) can also be expressed as a weighted C p (t)e ¡¯t term. Combining constant terms and rearranging the subscriptj, (2.5) can be rewritten as: ¹ C(t)= J X j=1 ® j e ¡¯ j t C p (t) (2.6) whereJ =2 ¹ J+1,® j are determined byM 1j ,M 2j andw j . Equation (2.6) is the general format used to model the TAC measured by a PET scanner, while (2.4) can be regarded as a special case withJ =3. 2 Assume C p (t) is non-negative, C p (0) = 0 and ¯ is a large number. 8t; t > 0, C p (¿) ¯e ¡¯¿ j ¿=t = R t 0 C p (¿)¯e ¡¯(t¡¿) d¿ = R t 0 C p (¿)de ¡¯(t¡¿) = C p (t)¡ R t 0 e ¡¯(t¡¿) dC p (¿). Because ¯ is a large number, e ¡¯(t¡¿) is close to zero outside the region [t ¡ ¢t;t], where ¢t is a small number. R t 0 e ¡¯(t¡¿) dC p (¿) ¼ R t t¡¢t e ¡¯(t¡¿) dC p (¿), because e ¡¯(t¡¿) ¼ 1 for t 2 [t¡ ¢t;t], R t t¡¢t e ¡¯(t¡¿) dC p (¿)¼ C p (t)¡C p (t¡¢¿)¼ 0. Therefore, the convolution result is approximately equal toC p (t). 10 0 500 1000 1500 2000 2500 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Time (sec.) Activity BG1 BG2 LES1 LES2 20 40 60 80 100 120 20 40 60 80 100 120 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 LES2 LES1 BG2 BG1 Figure 2.2: Clinical PET image and region of interest (ROI) TACs. Left: The last frame of a dynamic PET MAP reconstructed image. Right: TACs drawn from two lesion and two normal tissue ROIs. 2.2 Time-Frequency Properties of TACs It has been observed that FDG accumulation in tumors tends to rise with time, while accumulation in most normal tissues decreases, with the possible exception of bone marrow in fasting patients [WZH + 93, JWA + 95, ZW96]. As a result, the TAC of malig- nant tissue increases with time; whereas a normal tissue TAC diminishes with time. This pattern has been shown to be true in at least 80% of untreated primary cancers [WZH + 93, ZW96]. Figure 2.2 illustrates tumor and non-tumor TACs drawn from a lung cancer patient image. The image shown is the last frame (28 th ) of the dynamic sequence. Four 3£3 regions of interest (ROI) are chosen. LES1 is the primary lung tumor, LES2 is the metastasis, BK1 is non-tumor tissue in the mirror symmetrical posi- tion of LES2, and BK2 is non-tumor tissue adjacent to the metastasis. This figure shows that both tumor TACs increase with time, whereas non-tumor TACs diminish. In this section, we first analyze temporal differences (increasing vs. diminishing) in the frequency domain in order to characterize the behaviors of individual physiological process functions. Then, by decomposing the TACs into physiological process func- tions, or alternatively wavelet basis functions, we show that tumor and normal TACs 11 occupy different time-frequency bands, implying that tumor and normal TACs can be approximated by different basis function sets. These time-frequency domain properties are useful in the subspace identification described in section 2.3.2. 2.2.1 Physiological Processf(¯;t) vs. Macro Parameter The physiological process function [Sul93, NQAL02, CJ93]f(¯;t) is defined as: f(¯;t)=C p (t)e ¡¯t ;t¸0 (2.7) where C p (t) is the blood input function, and ¯, ¯ > 0, is the decay coefficient. Some constraints can be added to noise free C p (t): C p (t) is a non-negative function which first increases with time and then decreases. In mathematical terms this means: (i) C p (t) must be non-negative. i.e., C p (t) ¸ 0;t ¸ 0, and (ii) there exists a T 0 for each individualC p (t) such thatdC p (t)=dt¸0 fort·T 0 anddC p (t)=dt< 0 fort>T 0 . The first constraint is validated by the physical meaning ofC p (t). The second constraint has been verified in [FW93]. We prove that for a given time T;T > T 0 , (i) there exists a threshold value ¯ ¤ for each individual blood input function C p (t) such that @f(¯ ¤ ;t)=@tj t=T = 0, (ii) For ¯, 0 · ¯ < ¯ ¤ , @f(¯;t)=@tj T 0 <t·T > 0, and (iii) For ¯ > ¯ ¤ , @f(¯;t)=@tj t=T < 0 (see Appendix). Put simply, we can say that at time T , for ¯ · ¯ ¤ , f(¯;t) increases with time; whereas for¯ > ¯ ¤ , f(¯;t) decreases with time. Figure 2.3 illustrates how the function f(¯;t) changes when ¯ varies from 0 to 1. All curves, including C p (t), have been normalized to unit energy in order to fit them all in one plot. T 0 = 130 sec., T =3240 sec.,¯ ¤ =4:5£10 ¡4 (min ¡1 ). The figure also shows thatf(¯;t) is very close to the blood input function when¯ =1. 12 0 500 1000 1500 2000 2500 3000 3500 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.0e+000 1.8e−004 4.5e−004 1.2e−003 3.0e−003 7.9e−003 3.0e−002 Time (sec.) Activity f(β * , t) Blood Input Function t=T 0 t=T β<β * β>β * C p (t) β = 1.0e+000 Figure 2.3: The physiological process functionf(¯;t) vs. ¯. 2.2.2 Frequency Domain Properties of TACs With equation (2.6) and (2.7), the TAC ¹ C(t) measured by a PET scanner can be written as: ¹ C(t)= J X j=1 ® j e ¡¯ j t C p (t)= J X j=1 ® j f(¯ j ;t); t¸0 (2.8) From a signal processing point of view, the physiological function f(¯;t) is an output of the arterial input function C p (t) passing through a linear time invariant system with the impulse response h(¯;t) = e ¡¯t ; t ¸ 0. Applying the Fourier transform to both sides of equation (2.8), ¹ C(!)= J X j=1 ® j C p (!)H(¯ j ;!) (2.9) where C p (!) is the Fourier transform of C p (t); H(¯ j ;!) is the Fourier transform of h(¯ j ;t): H(¯;!)=1=(¯+j!); jH(¯;!)j=1= p ¯ 2 +! 2 (2.10) The smaller the value of ¯, the lower the filter’s cut-frequency. Malignant tissue TACs increase with time, whereas normal tissue TACs decrease with time [ZPL + 01]. 13 2 2 y[n] h H [n] h L [n] Stage 1 Channel 1, v 1 [n] 2 2 h H [n] h L [n] Stage 2 Channel 2, v 2 [n] Channel 3, v 3 [n] More details Figure 2.4: A two-stage wavelets analysis filter bank. When decomposing a TAC into the linear combination of f(¯;t) functions, the dom- inant f(¯;t) terms of a tumor’s TAC should be increasing functions and accordingly be associated with¯ < ¯ ¤ ; while the dominantf(¯;t) terms of a normal tissue’s TAC should be decreasing functions and accordingly associated with¯ > ¯ ¤ . Therefore, the TAC of malignant tissue has less high frequency components than normal tissue TAC. 2.2.3 Time-Frequency Tiling of Tumor and Normal TACs To verify the conclusion from the last section, the time-frequency tiling of TAC is plot- ted to show energy distribution in different frequency band and time slots. Wavelet decomposition is employed to decompose the TAC signal. Figure 2.6 illustrates that the time-frequency tilings calculated from clinical TACs are matched with properties derived in the last section. The discrete-time wavelet decomposition and synthesis can be expressed in vector format as: Decomposition :v =T a y Synthesis :y =T s v (2.11) 14 50 100 150 200 250 300 350 400 450 500 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Level=1 Level=2 Level=3 Level=4 Level=5 Level=6 Level=7 N Wavelet Base Function Value Figure 2.5: The Daubechies wavelet basis functions for a 512-point, seven-level filter bank. where y is the input signal, v is the vector containing wavelet coefficients, T a is the analysis matrix andT s is the synthesis matrix whose columnsÃ i ;i=1:::k, are discrete- time wavelet basis functions,k is the total number of basis functions: T s =[Ã 1 ;Ã 2 ;:::;Ã k ] (2.12) The discrete-time wavelet decomposition can be implemented using a subband lin- ear filter bank [VK95]. Suppose h L [n] represents the low-pass subband filter, h H [n] represents the high-pass subband filter, y[n] represents the input signal, an v i [n] rep- resents the output wavelet coefficient for channel i. Figure 2.4 illustrates a two-stage cascading structure of wavelet decomposition. The design ofh H [n] andh L [n] depends on which specific basis function set is chosen. In this section, Daubechies [VK95] basis function set and the MATLAB r wavelet toolbox function wfilters.m are used to com- pute wavelet coefficients. Figure 2.5 shows the Daubechies basis functions for seven levels (seven frequency bands), from high frequency (left) to low frequency (right). 15 0 500 1000 1500 2000 2500 0 0.005 0.01 0.015 0.02 0.025 0.03 Time (sec.) (a) TAC Value Original Normal TAC Denoised Normal TAC Time (sec.) (b) Frequency Band Number 500 1000 1500 2000 2500 1 2 3 4 5 6 7 8 High frequency 0 500 1000 1500 2000 2500 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Time (sec.) (c) TAC Value Original Tumor TAC Denoised Tumor TAC Time (sec.) (d) Frequency Band Number 500 1000 1500 2000 2500 1 2 3 4 5 6 7 8 High frequency Figure 2.6: Wavelet tilings of tumor and normal tissue TACs. In the left column: clinical noisy and de-noised TACs. In the right column: time-frequency tilings. The upper row is for normal tissue. The lower row is for tumor tissue. By plotting the amplitude of wavelet coefficients in a 2D format, we get the so called time-frequency tiling which illustrates the energy distribution of the signal in the time- frequency domain. Figure 2.6 illustrates the time-frequency tiling of normal and tumor TACs from the clinical data. The X-axis represents time and the Y-axis represents the frequency band number, which is the channel number in Figure 2.4. The de-noised TACs from the model-based least mean square error (LMSE) curve fitting are also included in the same figure to show the increasing or decreasing trend of the curves. Figure 2.6 shows that (i) the tumor TAC has fewer high frequency components than that of normal tissue, which is matched by the frequency domain properties mentioned in section 2.2.2, and (ii) the energy of the tumor TAC increases in the later time slots, while the energy of the normal TAC decreases, which is a direct inference from [ZPL + 01]. 16 2.3 Subspace Processing According to the analysis in the previous section, the TACs from different tissues can be approximated using different sets of basis function. In this section, the algorithms used to identify basis functions (linear subspace) will be discussed. The time-frequency distribution patterns of tumor and normal TACs are used in one subspace identification process. 2.3.1 The Formation of Subspaces Let N be the number of frames in a dynamic PET image sequence. t i , i = 1:::N, is the middle time point for each frame. The measured noisy TAC is represented by the N£1 column vectory =[y(1);y(2);:::;y(N)] T . Depending on the tissue type (tumor or normal), the TACy can be modeled as: Normal tissue : y =SÁ+n; S 4 =[s 1 ;s 2 ;:::;s q ] (2.13) Tumor tissue : y =Hµ+n; H 4 =[h 1 ;h 2 ;:::;h p ] (2.14) wheres 1 ;:::;s q are independentN£1 column vectors, which span aq-dimension linear subspace denoted byhSi, or in matrix formatS2R N£q . Á = [Á(1); Á(2);¢¢¢ ;Á(q)] T is the coordinate of the true but unknown normal TAC in subspace hSi. Similarly, h 1 ;:::;h p are independentN£1 column vectors, which span ap-dimension linear sub- space denoted byhHi, or in matrix formatH 2 R N£p . µ = [µ(1); µ(2);¢¢¢ ;µ(p)] T is the coordinate of the true but unknown tumor TAC in subspacehHi. It is important to note that p and q may not be the same. In this paper, the noisen is assumed to be additive Gaussian noise. 17 2.3.2 Subspace Identification Algorithms Given a P -pixel tumor or normal tissue region of interest (ROI), the subspace identifi- cation problem is to identifyH orS respectively. Using equation (2.14), for aP -pixel tumor tissue ROI, the data model can be written as, : Y =H£+N Y =[y 1 ;:::;y P ]; H=[h 1 ;:::;h P ] £=[µ 1 ;:::;µ P ]; N=[n 1 ;:::;n P ] (2.15) where Y is known, H, £ and N are unknown. The goal is to estimate H. We can divide the subspace identification algorithms into two categories: with time-frequency distribution prior knowledge derived in the previous section, or without. Subspace identification without prior knowledge 4-K Nonlinear Least Square Error (LSE) Fitting Modeling the averaged TAC with (2.2), the parametersK 1 »k 4 andV p can be estimated using the non-linear least square curve fitting technique. The details of this method are discussed in several publications [LH74, HPH + 80]. This method is widely studied and used, and commercial software is available. There are several drawbacks. Firstly, the blood input function is required. However, arterial blood sampling is not a routine clinical procedure, whereas the image-based ROI method cannot be performed if the imaging field of view (FOV) does not include the heart or any large artery. Secondly, the algorithm is very sensitive to noise and initial values, which may introduce a large estimation error. Principle Component Analysis We compute the singular value decomposition ofY, and rank order the singular values 18 as(¸ 2 1 ;e 1 ),(¸ 2 2 ;e 2 ),:::,(¸ 2 N ;e N ), where¸ 2 1 ¸ ¸ 2 2 ¸ :::¸ ¸ 2 N ¸ 0. We assume that the signal to noise ratio is high in the ROIs used in subspace identification. Therefore, the order of principal components is determined by the signal rather than noise, and the left most singular vectors correspond to the signal basis functions. We can then construct the tumor signal subspace matrix,H, with columns equal to the firstp left singular vectors computed from the tumor region. To determine an appropriate value for p, we use a minimum number of bases such that the accumulated energy exceeds a preset threshold value´ PCA : p=min p P p i=1 ¸ 2 i P N i=1 ¸ 2 i ¸´ PCA (2.16) We then select one or more ROIs for normal tissue and follow a similar procedure to identify an appropriate value forq and the normal tissue’s signal subspace matrixS. Subspace identification with prior knowledge Spectral Analysis Method In this method, the continuous format TACy(t) is modeled as [CJ93, MMC + 96]: y(t)= N k X k=1 ® k C p (t)exp(¡¯ k t); ®¸0; ¸·¯ k ·1 (2.17) whereN k is the maximum number of basis functions allowed in the model, and¸ is the decay constant for the radioisotope. The¯ k values are fixed and are chosen to cover the spectrum of expected kinetic behavior, from the slowest possible clearance ¸ (1:51£ 10 ¡4 sec. ¡1 ) to the fastest anticipated dynamic (e.g. 1 sec. ¡1 ). The problem here is to determine the values of ® k that best fit the measured data given a set of ¯ k values 19 predefined at the interval [¸;1]. With the given blood input function C p (t), the above equation can be rewritten in matrix format as: y =A®; A jk = 1 t j+1 ¡t j Z t j+1 t j C p (t)exp(¡¯ k t)dt (2.18) whereA2R N£N k andt 1 ;:::;t N+1 are the starting and ending frame times. The equa- tion can be solved for the coefficient vector ® using a non-negative least squares algo- rithm to minimize the cost function(A®¡y) T (A®¡y) [MMC + 96, LH74]. Only the convolution terms with significant coefficients are considered as basis functions, while the convolution terms with negligible coefficients are ignored. We modified the ¯ k range selection in our spectral analysis implementation as fol- lows. The¯ k range for a tumor is set as[¸;¯ ¤ ], whereas that for normal is set as[¯ ¤ ;1]. 15» 20¯ k values are spaced logarithmically within each range. The blood input func- tion C p (t) is used as a common basis function (ie. one column in matrix A) for both tissue types. 2.3.3 Subspace Identification Results Subspace identification procedures were tested on a clinical data set which was acquired using a CTI ECAT-953 scanner with a dynamic data sampling protocol (15sec£8, 30sec£4, 1min£11, 5min£5; total 40min, 28 frames). 128£128 images were recon- structed using the MAP algorithm [QLH + 98]. The patient had one primary lung tumor and one metastasis, which were confirmed by other clinical methods. Figure 2.2 shows the clinical image and TACs sampled from the selected ROIs. hHi was estimated from tumor region5£5 LES1 (Figure 2.2), andhSi was estimated from normal tissue region 5£ BG2 (Figure 2.2). The results show that the identified subspaces are accurate and the two subspaces are separable. 20 0 1000 2000 0 0.02 0.04 0.06 Time(sec) TAC intensity Primary Tumor LES1 TAC LES1 Rep. by <H> Rep. by <S> 0 1000 2000 0 0.01 0.02 0.03 Time(sec) TAC intensity Background BG TAC BG Rep. by <H> Rep. by <S> 0 1000 2000 0 0.01 0.02 0.03 Time(sec) TAC intensity Metastasis LES2 TAC LES2 Rep. by <H> Rep. by <S> Figure 2.7: A Comparison of the accuracy and separability of the subspaceshHi and hSi estimated using the spectral analysis method. 0 1000 2000 0 0.01 0.02 0.03 0.04 0.05 Time(sec) TAC intensity Primary Tumor LES1 TAC LES1 Rep. by <H> Rep. by <S> 0 1000 2000 0 0.01 0.02 0.03 0.04 Time(sec) TAC intensity Background BG TAC BG Rep. by <H> Rep. by <S> 0 1000 2000 0 0.01 0.02 0.03 0.04 Time(sec) TAC intensity Metastasis LES2 TAC LES2 Rep. by <H> Rep. by <S> Figure 2.8: A Comparison of the accuracy and separability of the subspaceshHi and hSi estimated using the PCA method. Spectral Analysis From the blood input function,¯ ¤ was estimated numerically as¯ ¤ = 1:37£10 ¡3 . The estimated basis functions forhHi wereC p (t)e ¡1:5£10 ¡4 t ,C p (t)e ¡1:05£10 ¡3 t and C p (t). The estimated basis functions forhSi wereC p (t)e ¡3£10 ¡3 t ,C p (t)e ¡1:44£10 ¡2 t and C p (t). Figure 2.7 shows the fitted curves using the identified subspaces. Table.2.1 shows the mean square errors of the curve fitting. Principle Component Analysis 21 The threshold value ´ was set at 5%. The resulting p and q values were both 3. Figure 2.8 shows the fitted curves using the identified subspaces. Table 2.1 shows the curve fitting errors. The MSEs in the PCA were on a magnitude order of one smaller than those using the spectral analysis method. However, this does not mean that the PCA technique is superior to the spectrum analysis method. This is because the noise pattern (Figure 2.8, the zigzag shape of the TACs) of LES1 TAC was fitted by the PCA basis functions, but was not fitted by the spectral analysis method. If we check the fitted LES2 TAC, PCA also produces a large MSE, because the noise pattern in the region LES2 is different from that in the region LES1. SA hHi hSi LES1 TAC 8:7£10 ¡4 1:5£10 ¡3 LES2 TAC 5:1£10 ¡4 6:5£10 ¡4 BG2 TAC 5:4£10 ¡4 1:3£10 ¡4 PCA hHi hSi LES1 TAC 2:9£10 ¡4 2:3£10 ¡3 LES2 TAC 4:3£10 ¡4 1:4£10 ¡3 BG2 TAC 1:3£10 ¡3 1:9£10 ¡4 Table 2.1: MSE for the TAC curve fitting. The left table shows the spectral analysis method. The right table shows the PCA method. Both the spectral analysis and PCA results show that the identified tumor subspace hHi can fit the tumor TAC reasonably accurate but cannot fit the normal TAC closely, whereas the identified normal subspacehSi can fit the normal TAC accurately, but can- not fit the tumor TAC closely. These results indicate that we can assume that the two identified subspaces are reasonably accurate and separable. 22 Chapter 3 The Matched Subspace Method The matched subspace detector is applied to a wide range of problems in radar, sonar and data communication, wherein the signal is constrained to lie in a multidimensional linear subspace. The matched filter detector is “matched” to a signal which is assumed to lie in a one-dimensional subspace. Thus it is “mismatched” to signals which lie outside this subspace. The obvious solution to this problem of mismatch is to match to signals that can lie in higher dimensional subspaces, see Figure 3.1 for an illustration of this concept. Herey is the observed signal,y 1 is the projection ofy in the one-dimensional subspace he 1 i which produces a large mismatch. When extending the one-dimensional subspace he 1 i to a two-dimensional subspacehe 1 ;e 2 i, the mismatch ofy 1;2 will be diminished. A further extension tohe 1 ;e 2 ;e 3 i will minimize the mismatch in this case. Detectors doing this are called the matched subspace detectors. The matched subspace detectors are suited to problems where the signal to be detected is comprised of a linear combination of modes whose weight are unknown. In light of the dynamic models described in the previous chapter, it is natural to apply the matched subspace detector to TACs for the lesion detection task. 3.1 Introduction Positron emission tomography (PET) with 2-[fluorine-18]-fluoro-2-deoxy-D-glucose (FDG) is improving the diagnosis, staging and treatment monitoring of a variety of can- cers relative to the alternative noninvasive imaging modalities [Con96, SC91]. Although 23 e 3 e 2 e 1 y y 1,2 y 1 Figure 3.1: Illustration for the matched subspace method. y is the observed signal, y 1 is the projection of y in the one-dimensional subspace he 1 i which produces large mismatch. Including the basis functions e 2 and e 3 in the subspace will diminish the mismatch. the conventional static PET imaging provides high sensitivity in tumor detection, further improvement is needed, because even a small percentage of false negatives result in high risk and cost. Conventional visual inspection and diagnosis is potentially inaccurate due to the limited spatial resolution and low lesion-to-background contrast, particularly for small tumors. This motivated the research into the development of algorithms for both static and dynamic PET to assist in lesion detection, which can help to assess cancer treatment efficacy or to detect tumors in early stages. Computer aided detection or diag- nosis has been widely studied for mammography, CT and MRI [GKA01, Doi07]. But there is considerably less application to PET images due to the relative inaccessibility of PET equipment [DCGJ99]. By using signal detection and image processing techniques, several algorithms for static PET lesion detection [HYB + 97, DCGJ99, YZS + 04] have been proposed to assist in visual inspection of tumors. While these algorithms work to some extent, the temporal domain information loss during static image reconstruction can reduce the lesion detectability. Since the dynamics of tracer uptake are different for normal and malignant tissue, the temporal domain provides potentially important additional features which may improve lesion detection performance. 24 Dynamic PET images are routinely used to estimate kinetic parameters for inves- tigation of delivery, transport and biochemical transformation of molecules within the body. Traditionally, kinetic parameter estimation is performed using non-linear regres- sion [LH74]. Based on a compartmental model [HPH + 80, HCH + 83], the quantifica- tion of individual model parameters, e.g. K 1 ;k 2 ;k 3 ;k 4 andV p for a three compartment model, can be estimated by fitting the model to the tissue time activity curve (TAC) over a selected region of interest (ROI). Differentiating primary colorectal tumors from nor- mal colon tissue based on quantitative parameters has been studied in [SKP + 07]. The result shows that combining all kinetic parameters for the classification analysis using support vector machine achieves high accuracy. However, due to the high noise level of single voxel kinetics, parameter images are usually of poor quality in small lesions. Computing parameters for each voxel throughout the whole image is still a challenging problem, especially when the blood input function is unknown. Patlak graphic analysis [PBF83, PB85] is a popular method to estimate the influx constant Ki = K 1 k 3 k 2 +k 3 using a simplified compartmental model. The Ki image is a single parameter approach that can assist in lesion detection since tumors often have higher influx rates compared to normal tissues [KV AN + 05]. However, use of Patlak analysis may be limited in small lesion cases, where the Ki parameter estimate is corrupted by partial volume effects from surrounding tissue. An alternative approach to the analysis of dynamic PET is to use a data driven method to extract clinically interesting information. Principle components analysis (PCA) has been widely used for this purpose [Bar80, PBBL94, TSDS + 03]. PCA is concerned with explaining the variance-covariance structure of a set of variables through a linear combinations of these variables [PBBL94]. In PCA, the dynamic image sequence is decomposed into the summation of orthogonal basis functions. Projecting the dynamic image sequence onto these basis functions produces PCA images. Visual 25 inspection of these images can assist in lesion detection, however simply using indi- vidual components for detection is limited by the fact that the dynamics of normal and malignant tissues are not orthogonal. In most circumstances, the projected energy of the tumor TAC spreads over several basis functions. It is difficult by visual inspection to combine these distributed energy to detect tumors. Research efforts have been made in dynamic PET based image segmentation. Sev- eral papers [WFMF02, BGYW03] use clustering techniques to automatically segment tissues in dynamic PET and potentially replace manual ROI delineation. Although it is not their main purpose, [WFMF02] gives examples of segmenting relatively large tumor regions. However, the homogenous assumption used in the paper may not hold between different tumors, especially between the primary tumor and metastases. Lack of an accu- rate noise model may also reduce it’s performance for low SNR cases.Consequently, this algorithm may not be an ideal tool to detect small metastasis in high noise. With the assumption of bilateral symmetry existing in normal brains but absent in pathological tissues, a segmentation based asymmetry feature tumor detection approach is proposed in [CFC03]. However, the bilateral symmetry may not be an appropriate assumption in chest. In addition, this algorithm also lacks an accurate noise model. In this chapter, we describe a computer aided detection method that applies statistical hypothesis testing to the time activity curves (TACs) at each voxel to distinguish tumors from background. We aim at detecting metastatic lesions, which often have low tumor- to-background contrast, small tumor size and are embedded in high statistical noise. In contrast to the most of algorithms mentioned above, this approach is formulated as a detection problem and thus produces a single statistic image which on thresholding indicates the lesion locations. A linear subspace model and a generalized likelihood ratio (GLR) test [Pro00, Sch90, Fuk90] are combined to derive a matched subspace 26 detector (MSD) [VO91, SF94, Sch90]. While developed for use in digital communica- tion systems, this approach has also been applied to activation detection in functional MRI [AKKK99, SdD05]. A lesion detection method developed for dynamic PET using matched subspace technique was presented in [HY99], where it was hypothesized that, for a given patient, tumor TACs may vary among different lesions (e.g. primary vs. metastatic lesions) but occupy approximately the same linear subspace. Similarly, TACs for normal tissue are assumed to occupy a second linear subspace which does not coin- cide with the tumor subspace. Using TACs from a primary tumor region of interest (ROI) and one or more normal ROIs, each identified by a human observer, these two linear subspaces can be identified. Applying the MSD to these subspaces on a voxel by voxel basis throughout the dynamic image produces a test statistic at each voxel which on thresholding indicates potential locations of metastatic lesions. The idea of presenting TACs of different tissues using different basis function sets is used in PCA. However, the purpose of PCA is to extract basis functions, the purpose here is to detect tumor signals given known basis functions. One difference between our method and the classical GLRT lies in how the likelihood ratio is defined. In the classical GLRT [AKKK99, SdD04, SdD05], the likelihood ratio (i) is defined as a ratio for a constrained parameter set £ 0 , £ 0 ½ £ over an unconstrained parameter set £ [CB01]; while in the alternative definition we use here, the likelihood ratio (ii) is defined as a ratio for one parameter set£ 1 over another parameter set£ 0 [Sch90, Pro00, Fuk90]. Consequently, definition (i) leads to F or Â 2 distributed test statistics, while definition (ii) does not. In the earlier MSD approach [HY99], the noise was assumed to be white and Gaus- sian with unknown variance. However, noise in reconstructed PET images will be time varying and spatially correlated. In this chapter we modify the data model to allow a time varying and spatially correlated noise model. The algorithm is also extended 27 to detection using multiple TACs within a small region of interest. A number of tech- niques have been described to estimate reconstructed noise covariance from a single data set [Hue84, Car93, QL00], depending on which reconstruction algorithm is employed. Here we use a frame-by-frame MAP reconstruction [QLH + 98] method and estimate the covariance using the plug-in method described in [QL00, LAQ + 04]. The normal and tumor subspaces are estimated from the training ROIs using data driven principal components analysis (PCA) [TSDS + 03]. To evaluate the sensitivity and specificity char- acteristics of this detector in comparison to others, we use a Monte Carlo simulation of a digital phantom and the receiver operating characteristic (ROC). We also present exam- ples of application of the proposed detection approach to clinical dynamic PET data for a breast cancer patient with metastatic disease. 3.2 Theory 3.2.1 Formulation of the Subspace Detection Problem Let N be the number of frames in a sequence of reconstructed dynamic PET images. The center time point for each frame ist 1 ;t 2 ;:::;t N . The measured TAC at a given voxel j is represented by anN£1 column vectory j =[y j (1);y j (2);:::;y j (N)] T . The lesion detection problem is formulated by deciding between two possible hypotheses. The null hypothesis H 0 is that y j consists of a normal tissue time activity curve x 0j and noise n 0j » N(0;§ 0j ). The lesion present hypothesisH 1 is thaty j consists of a tumor time activity curvex 1j and noisen 1j »N(0;§ 1j ), i.e., H 0 :y j =x 0j +n 0j H 1 :y j =x 1j +n 1j (3.1) 28 The noise covariance § 0j and § 1j may be known or unknown matrices depending on the selected noise model as discussed later. The amplitude and even the shape of the signal x 0j and x 1j may varies from voxel to voxel, but we assume that they obey the linear subspace model: x 0j =SÁ j x 1j =Hµ j (3.2) where S 2 R N£q , Á j 2 R q£1 and H 2 R N£p , µ j 2 R p£1 . The normal tissue’s TACx 0j lies in a q-dimensional subspace ofR N which is spanned by the columns of S and denoted byhSi. Á j = [Á j (1); Á j (2);¢¢¢ ;Á j (q)] T is the coordinate of the true but unknown normal TACx 0j with respect to the subspacehSi. Similarly, the true but unknown tumor’s TACx 1j lies in the p-dimensional subspacehHi with coordinate µ j . Methods for subspace identification forH andS are discussed in section 3.2.3. In the detection problem,H andS are assumed known. Substituting (3.2) into (3.1), we have the data model for our hypothesis test: H 0 :y j =SÁ j +n 0j H 1 :y j =Hµ j +n 1j (3.3) 3.2.2 Generalized Likelihood Ratio (GLR) Test The GLR test for the hypothesis test in (3.3) are derived below for three cases: Case 1, L 1 (y): Detection for a single voxel TAC, with the assumption thatn 0j andn 1j are white Gaussian noise with a single unknown variance parameter for each TAC. Case 2, L 2 (y): Detection for a single voxel TAC, with the assumption that the noise is Gaus- sian with known temporal covariance matrices § 0j and § 1j . Note that for frame by frame reconstruction, these covariances will be diagonal since the noise is independent 29 between frames, however the diagonal elements can vary from frame to frame. Case 3, L 2 (y ROI ): Detection for multiple TACs within a region of interest (ROI) with known spatial and temporal covariance. Case 1,L 1 (y): Detection for a single TAC with unknown white noise variance. Given the data model in (3.1) and assuming that the covariance matrices§ 0j =¾ 2 0j I and § 1j =¾ 2 1j I, where¾ 2 0j and¾ 2 1j are unknown, the GLR is given by: ^ l 1 (y j )= l( ^ µ j ;H;^ ¾ 2 1j I;y j ) l( ^ Á j ;S;^ ¾ 2 0j I;y j ) (3.4) wherel(®;M;§;y) is the likelihood function defined as: l(®;M;§;y)=(2¼) ¡N=2 j§j ¡1=2 ¢ exp ¡1=2(y¡M®) T § ¡1 (y¡M®) (3.5) Substituting§=¾ 2 I in (3.5) leads to the simplification: l(®;M;¾ 2 I;y)=(2¼¾ 2 ) ¡N=2 exp ¡1=2¾ ¡2 ky¡M®k 2 (3.6) Let ^ Á j , ^ µ j , ^ ¾ 2 0j and ^ ¾ 2 1j be the maximum likelihood estimates (MLEs) ofÁ j ,µ j ,¾ 2 0j and ¾ 2 1j , respectively. It can be shown that: ^ Á j =(S T S) ¡1 S T y j ; ^ ¾ 2 0j = 1 N ky j ¡S ^ Á 2 j k 2 ^ µ j =(H T H) ¡1 H T y j ; ^ ¾ 2 1j = 1 N ky j ¡H ^ µ 2 j k 2 (3.7) 30 Substituting (3.7) into (3.4): ^ l 1 (y j )= l( ^ µ j ;H;^ ¾ 2 1j I;y j ) l( ^ Á j ;S;^ ¾ 2 0j I;y j ) =( ^ ¾ 2 1j ^ ¾ 2 0j ) ¡N=2 (3.8) It is convenient to replace the GLR by the N=2 th root [Sch90] for the equivalent test statistic: L 1 (y j )= ^ l 1 (y j ) 2=N = ^ ¾ 2 0j ^ ¾ 2 1j = y T P ? S y y T P ? H y (3.9) whereP ? H = I¡H(H T H) ¡1 H T is the projection onto the orthogonal complement of subspacehHi andP ? S = I¡S(S T S) ¡1 S T the projection onto the orthogonal comple- ment ofhSi. Case 2,L 2 (y): Detection from a single TAC with known noise covariance In this case, S, H, § 0j and § 1j are assumed known. Methods for estimating these parameters from the data are discussed below. The GLR for this case is: ^ l 2 (y j )= l( ^ µ j ;H;§ 1j ;y j ) l( ^ Á j ;S;§ 0j ;y j ) (3.10) where ^ Á j =argmax Á j l(Á j ;S;§ 0j ;y j ) ^ µ j =argmax µ j l(µ j ;H;§ 1j ;y j ) (3.11) 31 Pre-whiteningH 0 andH 1 in (3.3) using§ ¡1=2 0j and§ ¡1=2 1j respectively: ~ y 0j =§ ¡1=2 0j y j ; ~ S 0j =§ ¡1=2 0j S ~ y 1j =§ ¡1=2 1j y j ; ~ H 1j =§ ¡1=2 1j H (3.12) we can rewrite (3.3) as: H 0 : ~ y 0j = ~ S 0j Á j +n H 1 : ~ y 1j = ~ H 1j µ j +n (3.13) wheren»N(0;I) is additive white Gaussian noise. The GLR is then: ^ l 2 (y j )= l( ^ µ j ;H;§ 1j ;y j ) l( ^ Á j ;S;§ 0j ;y j ) = l( ^ µ j ; ~ H 1j ;I;~ y 1j ) l( ^ Á j ; ~ S 0j ;I;~ y 0j ) (3.14) Following the derivation in [SF94], the maximum likelihood estimates of ^ Á j and ^ µ j are: ^ Á j =( ~ S T 0j ~ S 0j ) ¡1 ~ S 0j ~ y 0j ^ µ j =( ~ H T 1j ~ H 1j ) ¡1 ~ H 1j ~ y 1j (3.15) Substituting (3.15) into (3.14), and replacing the GLR with its log [Sch90], we have L 2 (y j )=2ln ^ l 2 (y j ) =ln( j§ 0j j j§ 1j j )+ ~ y T 0j P ? ~ S0j ~ y 0j ¡ ~ y T 1j P ? ~ H1j ~ y 1j (3.16) where P ? ~ S0j =I¡ ~ S 0j ( ~ S T 0j ~ S 0j ) ¡1 ~ S T 0j P ? ~ H1j =I¡ ~ H 1j ( ~ H T 1j ~ H 1j ) ¡1 ~ H T 1j (3.17) 32 Finally inserting (3.17) and (3.12) into (3.16) produces the test statistic for Case 2: L 2 (y j )=ln( j§ 0j j j§ 1j j )+y T j § ¡1 1j [H(H T § ¡1 1j H) ¡1 H T ]§ ¡1 1j y j ¡y T j § ¡1 0j [S(S T § ¡1 0j S) ¡1 S T ]§ ¡1 0j y j +y T j (§ ¡1 0j ¡§ ¡1 1j )y j (3.18) As described below, we estimate the noise covariance directly from the data for the iteratively reconstructed MAP images [QLH + 98] using the method described in [QL00]. The estimated covariance will be the same when computing the likelihood for the lesion and no-lesion cases (§ 0j =§ 1j =§ j ; 8j), in which case (3.18) simplifies to L 2 (y j )=y T j § ¡1 j [H(H T § ¡1 j H) ¡1 H T ¡S(S T § ¡1 j S) ¡1 S T ]§ ¡1 j y j (3.19) Case 3, L 2 (y ROI ): Detection for multiple TACs within an ROI with known spatial and temporal covariance. When computing the test statistic over a region of interest we must take the spatial as well as temporal correlation into account. The GLR in this case keeps the identical form of (3.18) but with vector y now the concatenation of the TACs for each voxel in the ROI. The subspace operatorsH andS are modified to account for the additional TACs through the Kronecker product of their single TAC counterparts with an identity matrix of dimension equal to the number of voxels in the ROI. Similarly the noise covariance matrices are modified to model both spatial and temporal covariance. Assume the ROI to apply the detection contains M voxels. Let m;m = 1:::M be the voxel index. Concatenating the TACs y m ; m = 1:::M, within the ROI to form a long column vectory ROI 2R N¢M£1 : y ROI =[y T 1 ; y T 2 ;:::; y T M ] T (3.20) 33 For hypothesis H 0 , the concatenated coefficient vector Á ROI 2 R N¢q£1 and the noise vectorn 0ROI 2R N¢M£1 are denoted by Á ROI =[Á T 1 ;Á T 2 ;:::;Á T M ] T n 0ROI =[n T 01 ;n T 02 ;:::;n T 0M ] T (3.21) whereÁ m andn 0m are the coefficient vector and the noise vector at voxelm respectively. Defining a block diagonal matrixS2R N¢M£q¢M : S=I M£M S=diag(S; S;:::;S | {z } ) with totalM sub-blocks (3.22) Therefore, the data model under hypothesisH 0 is rewritten as: H 0 :y ROI =SÁ ROI +n 0ROI (3.23) The noise termn 0ROI »N(0;§ 0ROI ), where the covariance matrix§ 0ROI 2R N¢M£N¢M can be partitioned intoM 2 sub-blocks§ 0lm 2R N£N ; l;m =1:::M, each represents the covariance matric betweenn 0l andn 0m . In this thesis, the MAP [QLH + 98] algorithm is used for the dynamic image reconstruction. The noise covariances are estimated using the method provided in [QL00]. We can also defineH,µ ROI ,n 1ROI and rewrite the data model for hypothesisH 1 : H 1 :y ROI =Hµ ROI +n 1ROI (3.24) 34 Inserting equation (3.18) with the newly formed data models forH 0 andH 1 , the detector on an ROI is obtained: L 2 (y ROI )=y T ROI § ¡1 1ROI [H(H T § ¡1 1ROI H) ¡1 H T ]§ ¡1 1ROI y ROI ¡y T ROI § ¡1 0ROI [S(S T § ¡1 0ROI S) ¡1 S T ]§ ¡1 0ROI y ROI +ln( j§ 0ROI j j§ 1ROI j ) (3.25) 3.2.3 Subspace Identification In the derivation of the detector, we do not impose any constraints on the form of the TAC subspaces, with the exception that H and S must be of full column rank. This allows flexibility in how the subspaces are selected. We assume that the subspaces will be identified through training on regions selected by a human observer for the presence of primary tumors (H) and background areas where we are confident there are no lesions (S). Using these training regions we can estimate the subspaces by locally fitting parameters to a compartmental model [HPH + 80], using a spectral analy- sis [CJ93, TML + 94] approach in which the order of the underling kinetic model need not be specified, or in an entirely data driven form by applying principle components analysis (PCA) [PBBL94, TSDS + 03]. The basic principles of these algorithms can be found in Chapter 2. In the results presented here we use PCA, which avoids the need to specify a dynamic model or have knowledge of the tracer input function. Consider an background ROI with P voxels. We arrange the TACs for these voxels as the columns of a matrix Y = [y 1 ;y 2 ;:::;y P ]. Now compute the singular value decomposition of Y and rank order the singular value - singular vector pairs as(¸ 2 1 ;e 1 ); (¸ 2 2 ;e 2 ); :::; (¸ 2 N ;e N ), where ¸ 2 1 ¸ ¸ 2 2 ¸ ::: ¸ ¸ 2 N ¸ 0. Here we assume that the signal to noise ratio is high in the ROIs used for subspaces identifications. Therefore, the order of principal components over these ROIs is determined by signal rather than noise. Therefore, the left most 35 singular vectors are corresponding to the signal basis functions. We then construct the background signal subspace matrix,S, with columns equal to the firstq left singular vec- tors computed respectively from background regions. To determine an appropriate value of q, we use the minimal number of bases such that the accumulated energy exceeds a preset threshold value´ PCA : q =min q P q i=1 ¸ 2 i P N i=1 ¸ 2 i ¸´ PCA (3.26) By selecting another ROI for the primary tumor and following a similar procedure, we can identify the tumor signal subspace matrixH. 3.3 Method 3.3.1 Implementation of the Matched Subspace Detectors The detectors were implemented and applied to simulated and clinical PET data as described in Section 3.4. We implemented the single voxel detectors (L 1 (y) andL 2 (y)) using TACs averaged over a 3£ 3 ROI centered at each voxel, and then scanned the detector over the entire image to produce a 2D spatial image of the GLR test statis- tic. Thresholding of the statistic can be used as the basis for computer assisted lesion detection. The implementation ofL 2 (y ROI ) does not follow equation (3.25). The main concern is that the computation ofL 2 (y ROI ) involves an inversion of the matrix§ ROI (assuming § ROI =§ 0ROI =§ 1ROI ), which is computationally intensive. For example, on a small3£3 ROI with 30 frames, § ROI is a (3¢3¢30)£(3¢3¢30) matrix. Because § ROI is not a block diagonal matrix and is spatially variant, the time consuming matrix inversion has to be calculated for each location. To solve this problem, we maintain separate TACs 36 The 5x5 region to apply pre- whitening for the center voxel The center voxel When the center voxel moves, the prewhitening region moves as well. Figure 3.2: Illustration of spatial pre-whitening. The 5£ 5 dash line grid represents the 25-voxel region applying spatial pre-whitening for the center voxel. We plot two 5£5 pre-whitening regions and corresponding center voxels to illustrate that the pre- whitening region moves with the center voxel. for each voxel and compute the GLR test statistic over the specified ROI. A frame-by- frame spatial pre-whitening operation precedes calculation of the likelihood ratio for L 2 (y ROI ). Because dynamic frames are independent, the pre-whitening on each frame can be done separately. The idea here is to rearrange the data so that the covariance matrix is a block diagonal matrix. This procedure reduces the computational load by an order of magnitude. To prewhiten the image on one frame, we calculate the covariance over a 5£ 5 voxel region and use its inverse to prewhiten the voxel at the center of the region as illustrated in Figure 3.2. We then shift the5£5 region across the image, pre-whitening each voxel in turn. This pre-whitening operation modifiesy,S,H to~ y, ~ S, ~ H respectively. On one voxel m, the hypothesis test is changed to: H 0 : ~ y m = ~ S m Á m +n H 1 : ~ y m = ~ H m µ m +n (3.27) 37 wheren 2 R N£1 is the unit variance white Gaussian noise. With the assumption that the noise covariances are the same forH 0 andH 1 , the GLR test on voxelm is: L 2 (y m )= ~ y T m ~ H m ( ~ H T m ~ H m ) ¡1 ~ H T m ~ y m ¡~ y T m ~ S m ( ~ S T m ~ S m ) ¡1 ~ S T m ~ y m (3.28) The terms ~ y m 2 R N£1 , ~ H m 2 R N£p and ~ S m 2 R N£q are computed as follows. Assume the matrix K m;n 2 R D£D is the spatial covariance matrix for the D-voxel pre-whitening region centered at voxel m on then th frame (in Figure 3.2, theD-voxel pre-whitening region is illustrated by the 5x5 dashed grid,D = 25). The vectory m;n 2 R D£1 is the spatial image in the pre-whitening region. The pre-whitening matrix is K ¡1=2 m;n . In practice,K ¡1=2 m;n is replaced by 1 2 [K ¡1=2 m;n +(K ¡1=2 m;n ) T ] to guarantee symmetry. Then we draw the pre-whitening matrix’s center row k m;n , which corresponds to the center voxel. ~ y m and ~ S m are computed by: ~ y m =[~ y m (1);~ y m (2);:::;~ y m (N)] T ~ y m (n)=k m;n y m;n (3.29) ~ S m = 2 6 6 6 4 ~ s (1) m . . . ~ s (N) m 3 7 7 7 5 ; ~ s (n) m 2 R 1£D ; ~ s (n) m =k m;n 2 6 6 6 4 s (n) . . . s (n) 3 7 7 7 5 (3.30) wheres (n) is the n th row of the matrixS. ~ H m can be computed similarly. Notice that the white noise n in (3.28) is independent with noise on other voxels. Therefore, the 38 GLR test L 2 (y ROI ) on the M-voxel ROI is the summation of the GLR test on each individually pre-whitened voxel: L 2 (y ROI )= M X m=1 ~ y m ~ H m ( ~ H T m ~ H m ) ¡1 ~ H m ~ y m ¡~ y m ~ S m ( ~ S T m ~ S m ) ¡1 ~ S m ~ y m (3.31) 3.3.2 Noise Covariance Estimation Images for each frame were reconstructed using the preconditioned conjugate gradient method for MAP estimation described in our earlier paper [QLH + 98]. For L 1 (y), the variance is assumed unknown so that no covariance estimates are needed. To compute the variances and covariances of the images for use in the detectorsL 2 (y) andL 2 (y ROI ) we used the plug-in approximation techniques described in [QL00] and [LAQ + 04]. These methods use the observed data to compute an estimate of the Fisher informa- tion matrix for the observed data, which is used in turn to estimate the covariance of the reconstructed MAP images. Monte-Carlo simulations [QL00], [LAQ + 04] indicate that the approximations are reasonably accurate. We assume that the lesions we are looking for are sufficiently small that their presence will have minimal impact on the covariance of the image [QH01]. Thus, under hypothesesH 0 andH 1 , we assume§ j =§ 0j =§ 1j , for all voxel indices j. Since we use an average over a 3£ 3 ROI in L 2 (y) we first compute the spatial covariance of the reconstructed images for each frame. We then use this to compute the variance of the 3£3 voxel average. Repeating this procedure for each frame results in an estimate of the time varying variance for each frame in the TAC at each voxel. ForL 2 (y ROI ) we compute the spatial covariance over the5£5 ROI for each frame and prewhiten as described above. 39 3.3.3 Receiver Operation Characteristic The receiver operating characteristic(ROC) curve is a standard tool to evaluate and com- pare different detection methods. By varying a threshold applied to the test statistic, two performance measures are estimated: the sensitivity or true positive fraction (TPF), and specificity or false positive fraction (FPF). The ROC curve is a plot of TPF vs. FPF as a detection threshold varies. When comparing two detection methods, the one whose ROC curve gives higher sensitivity at matched specificity for all points on the curve is the better detector. A simple metric to compare methods is the area under the curve (AUC), where the method with larger AUC is superior. The goal of the ROC study in this chapter is to determine, first, whether the addi- tion of dynamic information in the form of a voxel-wise TAC improves metastatic lesion detection performance in PET; and second, whether matched subspace detection is supe- rior to the graphical Patlak method, which also uses dynamic data, in the lesion detec- tion task. Consequently we present a study in which we compare false vs. true positive performance in an ROC study using simulated PET data for three classes of detectors: (i) the three matched subspace detectors described in Section 3.2; (ii) a 2D Hotelling observer applied to single frame reconstructions; (iii) detection based on theKi param- eter computed using Patlak analysis. We would expect the dynamic methods to outperform static detectors, since they have more information available, however it is important in this preliminary evaluation to compare against a method that reflects current clinical practice, in the sense that a single frame is used for detection. To optimize performance of the static detector in our comparisons we used the Hotelling observer [BGG + 92, BYRM93] with signal and background known exactly, since this is the ideal observer under our assumed Gaussian noise model. To select the optimal frame period for the static observer, we evaluated performance of the detector as a function of the start and end point of the frame and 40 chose the duration that maximized the area under the ROC curve for the study described below. To compute the filter template for the Hotelling observer, a signal-known-exactly, background-known-exactly (SKE/BKE) model was used. Letg st 0 andg st 1 denote the reconstructed static image under hypothesesH 0 andH 1 respectively, andg st denote the reconstructed static image. The test statistic for the Hotelling observer is [BGG + 92]: ¸ HO (g st )=[E(g st 1 )¡E(g st 0 )] T § ¡1 st g st (3.32) where§ st represents the spatial covariance matrix for the imageg st . It has been demon- strated that reconstructions of noiseless data provide reasonable approximations of the mean of reconstructed PET images [QL00, FR96]. We therefore computed the tem- plate [E(g st 1 )¡E(g st 0 )] from MAP reconstructions of noiseless data with and with- out lesions. The reconstructed spatial covariance matrix was estimated in the manner described above in section 3.3.1 [QL00]. In the implementation used to generate the results presented below, the observer template was computed over9£9 windows cen- tered at each of the tumor locations. By sliding this window throughout the entire image, we are also able to produce a test statistic image of the Hotelling observer for qualita- tive comparison with the outputs of the matched subspace detectors. To generate the ROC curve we consider the observer output at the true lesion locations only. A positive detection occurs if¸ HO >´, where´ is the decision threshold. Patlak analysis [PBF83, PB85] is based on a three compartment model in which a single uptake constant Ki 4 = K 1 k 3 k 2 +k 3 is computed from the observed TAC using linear regression. Larger values of Ki are often associated with tumors when compared to normal tissue [KV AN + 05] and hence a single parametric image can be used for tumor detection. Here, we model estimation of the blood input function required by the Patlak 41 method by sampling a blood pool ROI in the reconstructed image sequence. As with the subspace detectors, we reduce noise in the voxelwise estimates ofKi by averaging the TACs over a 3£3 ROI centered at each voxel before applying the Patlak method. Sliding the ROI window throughout the entire image produces a parametric image of ^ K i . In assessing performance for the ROC study, we threshold the value ^ K i only at each true lesion location to generate the ROC curve. A tumor is detected if ^ Ki > ´ at the lesion location, where´ is the threshold. The implementation of the matched subspace detection for three cases are described in section 3.3.1. The test statistics are computed in each3£3 ROI. Similar to generation of ROC curves for Patlak method, we consider the outputs at the true lesion locations only. Note that in contrast to the static Hotelling observer, the dynamic test statistic does not assume spatial distribution of the tumor is known, i.e. the detector uses only the temporal characteristics of the data for detection. While compared to Patlak method, the dynamic test statistic does not use blood input function. 3.4 Simulations, Experiments and Results 3.4.1 Digital Phantom Simulation With the aim of future in vivo evaluation using xenograft tumor models in mice, we simulated a 2D scanner with the same geometry as the Siemens Focus 220 small animal scanner. The reconstructed image size was128£128 and the sinogram dimension was 288£ 252. A total of 27 frames were simulated for a 60-minute dynamic scan. The last 20 min (40min - 60min), which maximized the area under the ROC curves for the Hotelling observer, was used for the static reconstruction. The scan protocols for both dynamic and static scans are shown in Tab.3.2. Each frame was reconstructed using 42 0 500 1000 1500 2000 2500 3000 3500 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Time (Sec.) Activity, CPS/pixel Noiseless Blood Reconstructed Blood Noiseless Primary Tumor Reconstructed Primary Tumor 0 500 1000 1500 2000 2500 3000 3500 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Time (Sec.) Activity, CPS/pixel Noiseless Background Reconstructed Background Noiseless Metastasis Reconstructed Metastasis Figure 3.3: TACs used in the simulations. Left: The TACs of the blood input function and the primary tumor in Phantom 1. Right: The TACs of normal tissue and metastasis in Phantom 2. All phantoms use the same true TAC for all normal tissue. The noisy TACs shown are for a3£3 ROI drawn from a single reconstruction for the data used in the ROC study. MAP [QLC + 98] initialized with two-iterations of OSEM with 9 subsets followed by 20 iterations of the MAP PCG algorithm. Three digital phantoms were generated: Phantom 1, Figure 3.4.(1), which contained blood pool, normal tissue (background), and the primary tumor, was used to extract the blood input function for the Patlak method and estimate ^ S and ^ H for matched subspace detection. Phantom 2, Figure 3.4.(2), consisting of normal tissue and six small lesions, representing the case of metastatic disease in the patient. Placing multiple metastases in one phantom saves time in both sinogram simulation and image reconstruction. Since the six metastases are spaced reasonably far apart, the noise in each, and hence prob- ability of detection, is effectively independent. The hypothesis test was performed on each metastasis. Phantom 3 (not shown) was composed of background only, i.e. the same as Figure 3.4.(2) without the six lesions. To simulate real physiological processes, TACs assigned in the digital phantom were obtained from a clinical dynamic FDG-PET breast cancer study. ROIs were drawn in the clinical images to define background, pri- mary tumor and metastatic regions. The micro-parameters in each type of tissue were estimated by fitting to a three compartment model [HPH + 80] using a Newton-Raphson 43 Primary Tumor Background Heart (Blood) (1.a) Background 6 Metastases (2.a) (1.b) (2.b) (1.c) (2.c) Figure 3.4: Digital phantom images (for all dynamic images, only the last frame is shown): (1) Digital phantom 1, which contains normal tissue, one primary tumor and blood pool; (2) Digital phantom 2, which contains normal tissue and six metastases; Phantom 3 contains pure background and is not shown here. (a) The noise free digital phantom; (b) A reconstructed dynamic image; (c) A reconstructed static image. The frame duration was chosen such that this image has the maximum SNR achievable (in terms of the Hotelling observer) by any single frame static image. 44 Table 3.1: Micro Parameters for True TACs Used in Digital Phantoms Normal Tissue Primary Tumor Metastasis K 1 , ml/g/min 4.5x10 ¡2 3.6x10 ¡2 2.7x10 ¡2 k 2 , /min 3.1x10 ¡1 2.1x10 ¡1 7.2x10 ¡2 k 3 , /min 3.4x10 ¡3 1.2x10 ¡1 2.4x10 ¡2 k 4 , /min 2.0x10 ¡4 1.8x10 ¡3 6.0x10 ¡4 V p , ml/g 1.1x10 ¡1 1.7x10 ¡2 1.6x10 ¡2 procedure [LH74]. The arterial blood activity curve was approximated using an ROI in the left ventricle [GSH + 89]. The resulting micro-parameters are listed in Tab. 4.1. Table 3.2: Monte Carlo Simulation DATA Acquisition Protocol Scan Type Scan Interval Frame Duration Frame Number Dynamic scan 0 - 5 min 30 sec/frame 10 Dynamic scan 5 - 8 min 1 min/frame 3 Dynamic scan 8 - 20 min 2 min/frame 6 Dynamic scan 20 - 60 min 5 min/frame 8 Static scan 40 - 60 min 20 min 1 Substituting the arterial blood curve and estimated micro-parameters back into the three compartment model yielded the fitted TACs in the primary tumor, metastases, and normal tissues. The fitted and reconstructed TACs are shown in Figure 3.3. The fitted primary and metastasis TACs were used in the digital phantoms while the TACs of normal tissue in the phantom were scaled by 2 from the clinical results in order to lower the lesion to background contrast, so that the metastasis became more difficult to discern in the reconstructed images. The total counts in the 2D phantom 3 (pure background) was 0:37£10 6 for the last 5-min dynamic frame and 1:57£10 6 for the static 20-min frame. 45 The principle component analysis method was employed to estimate rank 3 sub- spaces ^ H and ^ S from the MAP reconstructions of simulated Phantom 1, see Figure 3.4(1.b). The lesion subspace ^ H was estimated from a 7£7 ROI centered in the pri- mary tumor. The normal subspace ^ S was estimated from a 7£7 ROI centered in the normal tissue region. The PCA method does not need a blood input function, therefore the proposed detection procedure can potentially be applied in clinical cases where the left ventricle is not present in the acquired images. The blood input function needed by the Patlak method was extracted by averaging the TACs in a7£7 ROI centered in the blood pool region as described previously. We compared detection results for five methods: (a) the Hotelling observer for static images; (b) Patlak analysis in3£3 ROIs in which the estimatedKi parameter is thresh- olded for lesion detection; (c)L 1 (y): matched subspace detector with an unknown vari- ance white noise data model; (d) L 2 (y): the matched subspace detector with known noise covariance; and (e)L 2 (y ROI ): the matched subspace detector applied to multiple TACs with known noise covariance. Examples of the 2D maps of the statistics for each of these detectors are shown Figure 3.4.1 and Figure 3.6. In the ROC study we use the statistic value only at the center of the six true lesion locations to determine the presence or absence of a lesion. To perform the quantitative ROC study, 200 pseudo-random dynamic Poisson data sets were generated for Phantoms 2 and 3. In each reconstructed image, separate hypoth- esis tests were conducted on the six possible metastasis locations. There were a total of 6£200 = 1200 lesion present cases and an equal number of lesion absent cases. ROC curves for each detector are plotted in Figure 3.7. We also tabulate the areas under the ROC curves and probability of detection as a function of the probability of false positive in Tab. 3.3. These results show that all of the dynamic detectors outperform the static observer. TheL 1 (y) detector performs similarly to, although somewhat better than, the 46 (a) Hotelling Row#65 Row#80 (b) Patlak 0 20 40 60 80 100 120 140 0 50 100 150 200 250 300 350 Row#65 Row#80 0 20 40 60 80 100 120 140 −1 −0.5 0 0.5 1 1.5 2 x 10 −4 Row#65 Row#80 Figure 3.5: Simulation results for Hotelling observer and Patlak method. (a) Static image Hotelling observer on9£9 window. (b) Patlak method applied to3£3 averaged TAC. The upper row shows 2D test statistic images for the Digital Phantom 2 (Figure 3.4 (2.a)). The lower row shows profiles through these images at the level of the two sets of three metastases. Patlak analysis. However, when introducing more accurate models for the uncertainty in the TACs with methods L 2 (y) and L 2 (y ROI ), we see substantial gains in performance. If the total count number is increased, the performances of all detectors will improve accordingly, however the ranking in terms of performance remains unchanged. It should be noted that all of the dynamic detectors use only information extracted from the data and reconstructed images while the static Hotelling observer uses information about the true background and lesion activity and shape. The simulation results show, perhaps unsurprisingly, that the dynamic detectors out- perform a static detector, even if the static detector assumes the shape and activity of the tumor is known while the dynamic detectors make no prior assumption about the 47 (c)L 1 (y) Row#65 Row#80 (d)L 2 (y) (e)L 2 (y ROI ) 0 20 40 60 80 100 120 140 0 0.5 1 1.5 2 2.5 Row#65 Row#80 0 20 40 60 80 100 120 140 −16 −14 −12 −10 −8 −6 −4 −2 0 2 4 Row#65 Row#80 0 20 40 60 80 100 120 140 −120 −100 −80 −60 −40 −20 0 20 40 Row#65 Row#80 Figure 3.6: Simulation results for matched subspace detection methods. (c) L 1 (y): matched subspace detector with unknown variance white noise data model. The detector is applied to 3£3 averaged TAC. (d) L 2 (y): matched subspace detector with known noise covariance. The detector is applied to3£3 averaged TAC. (e)L 2 (y ROI ): matched subspace detector with known noise covariance. The detector is applied to the 9 TACs within a 3£3 window. The upper row shows 2D test statistic images for the Digital Phantom 2 (Figure 3.4 (2.a)). The lower row shows profiles through these images at the level of the two sets of three metastases. Table 3.3: P D ,P FA and AUC Values for ROC Study. P D is the probability of detection. P FA is the probability of false positive. Methods P D P D P D AUC (P FA =5%) (P FA =1%) (P FA =0:1%) Hotelling 74.7% 57.3% 41.3% 1¡4:67£10 ¡2 Patlak 83.3% 66.5% 43.7% 1¡3:02£10 ¡2 MSDL 1 (y) 86.8% 74.0% 55.5% 1¡2:43£10 ¡2 MSDL 2 (y) 98.5% 95.0% 88.1% 1¡2:92£10 ¡3 MSDL 2 (y ROI ) 99.5% 98.5% 95.3% 1¡8:21£10 ¡4 48 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P FA P D L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer Figure 3.7: Lesion detection ROC curves computed from Monte Carlo simulations. The method names shown in the legend are ranked, from high to low, based on AUC values. tumor shape. The results also show the matched subspace detectors outperform the Pat- lak method, especially for the known noise covariance cases. Furthermore, the signal subspaces for the TACs used in the detector are estimated directly from the data without using the blood input function. Consequently the signal subspace detector is an entirely data driven procedure that can be applied independently to each new dynamic study for detection of metastatic lesions provided primary tumor and background ROIs are identified by a human observer for determination of the subspacesH andS. 3.4.2 Physical Phantom Study Detection errors have two major sources. One is the effect of noise on the computed statistic. Another is a systematic error resulting from inaccurate subspace identifica- tion. Here, we want to test the matched subspace detection under the circumstance that the subspace identification is accurate. Moreover, the detector is applied to physical phantom data collected by a real scanner. 49 0 1000 2000 3000 4000 5000 6000 7000 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Time, (sec) Activity e −β 1, Background e −β 2, Lesions Figure 3.8: Noiseless TACs in the physical phantom. The decay coefficient for the background TAC, represented by the blue curve, is ¯ 1 = 8:178£ 10 ¡4 . The decay coefficient for the lesion TAC, represented by the red curve, is¯ 2 =1:518£10 ¡4 . The physical phantom (see Figure 3.9) was built using two tracer solutions. 11 C- acetate with the decay coefficient ¯ 1 = 8:178£ 10 ¡4 sec ¡1 was used for the kidney shaped normal tissue region. FDG with the decay coefficient¯ 2 = 1:518£10 ¡4 sec ¡1 was used for three spherical lesions: large size Lesion 1, medium size Lesion 2 and small size Lesion 3. Because the true TACs in this physical phantom are pure expo- nential functions with known decay coefficients (see Figure 3.8), the one-dimensional subspacesH andS used in the detector are error free. The physical phantom was scanned using a CTI ECAT-953 scanner with dynamic data acquisition protocol (15sec£8, 30sec£4, 1min£21, 5min£3; total 40min, 36 frames). The images were reconstructed frame-by-frame using MAP algorithm. The subspace identification procedure was skipped because the normal and tumor TAC’s subspaces are known. The reconstructed images and L 1 detection result are shown in Figure 3.9. The results show that given accurate subspaces, the matched subspace detec- tor is effective in detecting lesions, even under low SNR conditions. This result shows 50 20 40 60 80 100 120 20 40 60 80 100 120 Background Lesion 1 20 40 60 80 100 120 20 40 60 80 100 120 Lesion 3 Figure 3.9: Physical phantom study. Top left: The plane containing the large tumor (Lesion 1). Top right: The plane containing the small tumor (Lesion 3). Bottom: L 1 result for Lesion 3 plane. our detector working well with experimental data in addition to the simulated data pre- sented in the earlier ROC study. 3.4.3 A Clinical Example Performance Demonstration We conclude with an example of this detector applied to a breast cancer patient with metastatic disease. The data were acquired using a CTI ECAT-953 scanner with a dynamic data sampling protocol (15sec£8, 30sec£4, 1min£21, 5min£3; total 40min, 36 frames). The reconstructed image size was 128£ 128 with 31 planes. A total of 51 36 frames were collected over 40 minutes. A MAP reconstruction of the last frame (a 5min scan starting at 35min) is shown in the 1st column of Figure 3.10. The transverse plane in the first row shows one primary tumor and that in the second row shows a metastatic lesion. Coronal and sagittal sections through the metastasis are show in the subsequent rows. ROIs were drawn in this plane to estimate the tumor and normal sub- spaces. Because different non-tumor (normal) tissues have different TACs, two ROIs were allocated in the image (see Figure 3.10, the upper left image) for normal subspace identification. Again, in order to avoid using a blood input function, PCA was used for subspace identification. The blood input function needed by the Patlak method was drawn from an ROI in the left ventricle (not shown). The outputs formed by scanning the detectors across the image volume are shown in columns 2»5 of Figure 3.10 for the Patlak method and the three matched subspace detectors. These results indicate that the dynamic detectors can enhance contrast in metastatic lesions relative to other areas, and hence illustrate its potential for assisting in the detec- tion of metastases. The results also show the matched subspace detectors appear to outperform the Patlak method, especially for the known noise covariance cases, in the sense that the statistic output is maximum, in the three orthogonal planes shown, at the confirmed location of the metastatic lesion. BothL 1 andL 2 detectors suppress normal tissue background, including the heart, to enhance the signal to noise ratio in both the primary and metastatic lesions. In contrast, while the Patlak method does enhance con- trast in the lesion, other areas of normal activity also produce relatively large values of Ki. The result for the L 1 detector shows false alarms near the body surface. These are probably due to a combination of motion and attenuation correction related artifacts and the limitations of the white noise assumption. In fact, the TAC curves in these surface voxels showed increases with time, which mimics behavior in the tumor. Therefore, the 52 0 0.3 0.1 0.2 Metastasis Metastasis Metastasis 0 4.5e-3 1.5e-3 3.0e-3 0.8 2.9 1.5 2.2 -0.3 2.4 0.6 1.5 Primary Tumor BG ROI 1 BG ROI 2 -3 23 6 14 0 0.15 0.05 0.1 0 2.1e-3 0.7e-3 1.4e-3 0.8 1.8 1.1 1.5 -0.3 1.8 0.4 1.2 -3 18 4 11 0 0.27 0.09 0.18 0 0.10 0.03 0.06 0 4.5e-3 1.5e-3 3.0e-3 0 1.8e-3 0.6e-3 1.2e-3 -3 18 4 11 -3 18 4 11 -0.3 1.8 0.4 1.2 -0.3 1.8 0.4 1.2 0.8 1.8 1.1 1.5 0.8 1.8 1.1 1.5 Transverse Transverse Coronal Sagittal Artifact Artifact Artifact Image Patlak L 1 (y) L 2 (y) L 2 (y ROI ) Primary Tumor Metastasis Metastasis Metastasis Figure 3.10: Clinical breast cancer study. 1st row: A transverse plane containing the primary tumor. The ROIs used for subspace identification are shown on the left most image. 2nd row: Transverse plane containing metastatic lesion. 3rd row: Coronal plane containing the metastasis. 4th row: Sagittal plane containing the metastasis. matched subspace detectors tended to produce high values at these voxels. These effects were reduced using the L 2 detectors, since in these cases the covariance is estimated from the data, and the higher uncertainty at the tissue boundary resulted in a lowered likelihood ratio. The L 2 detectors do still show some regions in which the likelihood ratio is relatively high, although with a smaller value than at the confirmed metastatic lesion location. These could result in false positives, depending on the threshold that is used for detection. Metastasis Early Detection In addition to the tumors shown in Fig.3.10, the proposedL 2 detector also detected one 53 suspicious metastasis which was not diagnosed by the doctor. However, based on a single scan, we can not tell whether the detected signal is a false alarm or a true positive. The patient got a follow-up scan several months later without any therapy, which showed more metastases. Using the follow-up scan data, the doctor confirmed that the metastasis detected by the L 2 detector based on the first scan was a true positive. The first scan image and the detection result were co-registered to the follow-up scan image using RView 1 and are shown in Fig.3.11. This result demonstrates that the proposed L 2 detector can detect metastasis earlier than the regular diagnosis method based on PET and CT (The CT image used in the diagnosis is not shown here). 3.5 Analysis and Discussion 3.5.1 Detectability vs. the Number of Prompts The performance of the proposed matched subspace detection depends on the image SNR, while the image SNR changes with the number of prompts. Understanding how the number of prompts affects the performance of the proposed detection algorithm is important. In this subsection, we use the Monte-Carlo simulation to give an numerical answer to this question. The digital phantoms, simulation settings and ROC study are the same as those in 3.4.1. The true phantom images were scaled such that the number of prompts at the last frame was set to 0.2 million, 0.4 million, 0.6 million and 1 million respectively. ROC curves for each prompt setting are plotted in Figure 3.12. We also tabulate the AUCs and probability of detection (P D ) as a function of the probability of false positive (P FA ) in Tab.3.4. These results show that, for all detection method the detectability will increase 1 RView is a volumetric display and rigid co-registration tool. It is available online at http://www.colin- studholme.net/software/rview/rvmanual/rvman.htm 54 Transverse Sagittal Coronal A Transverse Sagittal Coronal Detection Follow-up Scan Image Fused Image 1st Scan Image Follow-up Scan Image Fused Image Detected Metastasis B Figure 3.11: Clinical breast cancer metastasis early detection. (A) The L 2 detection results based on the 1st scan data, the follow-up scan images and the fused images. (B) The first scan images, the follow-up scan images and the fused images. The red arrow indicates the position of the detected metastasis. 55 The Number of Prompts=0.2M 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P FA P D L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer The Number of Prompts=0.4M 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P FA P D L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer The Number of Prompts=0.6M 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P FA P D L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer The Number of Prompts=1.0M 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P FA P D L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer Figure 3.12: Lesion detection ROC curves change with the number of prompts. The method names shown in the legend are ranked, from high to low, based on AUC values. 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Prompt Number at the Last Frame (Million) P D , (for P FA =1%) L 2 (y ROI ), known Σ L 2 (y), known Σ L 1 (y), white noise with unknow σ 2 Patlak Hotelling Observer Figure 3.13: Lesion detection ROC curves for the different number of prompt settings. The method names shown in the legend are ranked, from high to low, based on AUC values. 56 as the number of prompt increases. The matched subspace detection algorithms, espe- ciallyL 2 detectors, have substantial advantage over other methods for all four cases. 3.5.2 False Alarms and Image Artifact The result for the L 1 detector shows false alarms near the body surface. These are probably due to a combination of motion and attenuation correction related artifacts and the limitations of the white noise assumption. In fact, the TAC curves in these surface voxels showed increases with time, which mimics behavior in the tumor. Therefore, the matched subspace detectors tended to produce high values at these voxels. This issue became severer forL 1 detector due to two reasons. First,L 1 test statistic is of the form of a ratio, which is invariant to amplitude changes. Signals laying in the subspace with large or small amplitude will be equally detected. Second, L 1 detector does not use prior knowledge of noise statistics and could be more vulnerable to noise interference. These effects were reduced using the L 2 detectors, since in these cases the covariance is estimated from the data, and the higher uncertainty at the tissue boundary resulted in a lowered likelihood ratio. The L 2 detectors do still show some regions in which the likelihood ratio is relatively high, although with a smaller value than at the confirmed metastatic lesion location. When the anatomic images are available, like those provided by PET-CT systems, the false alarms appearing near the body surface or between arms can be screened out. 3.5.3 Conclusions We have described generalized likelihood ratio tests for lesion detection using matched subspace detectors applied to dynamic PET data. A computer observer ROC study was used to evaluate and compare the matched subspace detection algorithms with a static Hotelling computer observer and a Patlak-based method for lesion detection. The area 57 Table 3.4: P D vs. P FA for the Different Number of Prompt Settings The Number of Prompts = 0.2 Million Methods P D P D P D AUC (P FA =5%) (P FA =1%) (P FA =0:1%) Hotelling 15.1% 31.9% 52.8% 1¡1:17£10 ¡ 1 Patlak 14.4% 37.1% 60.4% 1¡9:32£10 ¡ 2 L 1 (y) 35.0% 55.3% 74.0% 1¡5:49£10 ¡ 2 L 2 (y) 55.5% 73.4% 88.4% 1¡1:96£10 ¡ 2 L 2 (y ROI ) 75.0% 87.4% 95.3% 1¡7:98£10 ¡ 3 The Number of Prompts = 0.4 Million Methods P D P D P D AUC (P FA =5%) (P FA =1%) (P FA =0:1%) Hotelling 27.9% 45.0% 67.0% 1¡6:69£10 ¡ 2 Patlak 36.4% 60.7% 83.2% 1¡3:45£10 ¡ 2 L 1 (y) 48.4% 69.8% 88.9% 1¡2:57£10 ¡ 2 L 2 (y) 79.8% 95.9% 98.7% 1¡2:88£10 ¡ 3 L 2 (y ROI ) 96.2% 98.7% 99.7% 1¡6:01£10 ¡ 4 The Number of Prompts = 0.6 Million Methods P D P D P D AUC (P FA =5%) (P FA =1%) (P FA =0:1%) Hotelling 51.1% 68.2% 87.8% 1¡2:44£10 ¡ 2 Patlak 68.8% 83.8% 96.1% 1¡8:72£10 ¡ 3 L 1 (y) 86.4% 94.1% 98.3% 1¡4:14£10 ¡ 3 L 2 (y) 98.8% 99.6% 99.9% 1¡2:35£10 ¡ 4 L 2 (y ROI ) 99.8% 100.0% 100.0% 1¡2:54£10 ¡ 5 The Number of Prompts = 1 Million Methods P D P D P D AUC (P FA =5%) (P FA =1%) (P FA =0:1%) Hotelling 83.6% 92.2% 98.0% 1¡4:03£10 ¡ 3 Patlak 91.4% 96.2% 99.1% 1¡1:71£10 ¡ 3 L 1 (y) 96.3% 99.2% 99.8% 1¡6:49£10 ¡ 4 L 2 (y) 100.0% 100.0% 100.0% 1¡2:08£10 ¡ 6 L 2 (y ROI ) 100.0% 100.0% 100.0% 1¡7:98 58 under the ROC curve indicates the following descending order of detection performance: matched subspace detector over a small ROI with known covariance, the matched sub- space detector with known covariance, the matched subspace detector with unknown but constant variance, the Patlak method, and the Hotelling static observer. An example with application to a clinical breast cancer study indicates that this approach has poten- tial to assist in metastatic lesion detection in a clinical setting. While it is unsurprising that the dynamic observer outperforms a static observer, this result does provide support for the argument that lesion detection may be enhanced by using dynamic information. The simulation and clinical example also indicate that the proposed matched subspace detection will outperform the parametric Patlak algorithm for the lesion detection task, without needing access to the blood input function. 59 Chapter 4 Controlling Familywise Error Rate in Thresholding Detection of small lesions in fluorodeoxyglucose (FDG) positron emission tomogra- phy (PET) is limited by image resolution and low signal to noise ratio. In chapter 3, we described a matched subspace detection method that uses the time activity curve (TAC) to distinguish tumors from background in dynamic FDG PET. Applying this algo- rithm on a voxel by voxel basis throughout the dynamic image produces a test statistic image or “map”, which on thresholding indicates the potential locations of secondary or metastatic tumors. In this chapter, a thresholding method with a given familywise error rate (FWER) for the matched subspace detection statistic map is proposed. The statistic map is modeled as an inhomogeneous Gaussian random field, which we have validated using simulations. The thresholding procedure includes three steps. Firstly, the PET image is segmented into several homogeneous regions. Then the statistic map is normalized to a zero mean unit variance Gaussian random field. Finally a maximum statistic thresholding with a fixed FWER is applied. Because the lesions, if they exist, are assumed to be small, the distribution and smoothness of the statistic map remains the same, whether a lesion is present or absent. Therefore, the FWER thresholds for the maximum statistic are estimated from the statistic map to be tested, as if they are estimated from the statistic map in a case where there is no lesion. We demonstrated and evaluated the proposed thresholding method using digital phantoms generated from 60 clinical dynamic images. We also show an example of the application of the proposed approach to clinical PET data from a breast cancer patient with metastatic disease. 4.1 Introduction In chapter 3, we described a computer aided detection method that applies statistical hypothesis testing to the dynamic PET time activity curves (TACs) at each voxel to distinguish tumors from background. The result of this method is a continuous val- ued statistic map, which on thresholding indicates possible locations of secondary or metastatic tumors. An observer receiver operating characteristic (ROC) study demon- strated the superior performance of the matched subspace detector over static PET Hotelling observer [BYRM93] and dynamic PET Patlak [PB85] methods in lesion detection tasks. In this chapter, we aim to detect significant tumor signals in a statis- tic map with a preset false alarm rate. In cases where we have no prior information about distributions under the H 0 and H 1 hypotheses, statistical maps must be analyzed everywhere voxel by voxel for sta- tistically significant activations. In the medical imaging community, this problem has been studied to locate the significance of the apparent signal observed in noisy images, and various developed methods have been used for different applications. In PET, [WEMN92, WMN + 96] familywise error rate (FWER) control methods have been devel- oped to detect the activation regions in cerebral blood flow (CBF) difference images under baseline and stimulation conditions. In functional magnetic resonance imaging (fMRI), [FWF + 94, FAH + 96] studied the algorithms to obtain delineation and detect activation areas involved in cognitive processes. In magnetoencephalography (MEG), [PNBL05] studied the objective assessment of current density maps (CDMs) to identify regions of activation. A survey of these techniques can be found in [Lan96, NH03]. 61 The analysis of statistic images involves testing hypotheses at each of thousands of volume elements, or voxels, for statistically significant experimental effects. This raises the possibility of a large number of false positives as a result of multiple hypothesis test- ing. For example, the naive thresholding of 10,000 independent voxels at an ® = 5% threshold would be expected to generate 500 false positives. To effectively control the number of false positives over all tests, we must therefore consider the multiple hypoth- esis testing problem. Many false positive measures have been proposed in this context, including the familywise error rate (FWER), expected false discovery rate (FDR), per- comparison error rate and per-family error rate [WMN + 96, NH03]. In this study, we choose the familywise error rate (FWER), i.e. the probability of one or more false pos- itives when testing multiple hypotheses (one per voxel throughout the image), as the false positive measure. The simplest approach to controlling FWER is the Bonferroni correction method [NH03, HT87]. This method produces conservative thresholds, unless the tests are inde- pendent which is rarely true in practice. Other methods that consider the spatial depen- dence of the data make inferences based on the global maximum distribution. FWER is directly related to the maximum statistic. One or more voxels, T i , will exceed the thresholdt under the null hypothesisH 0 only if the maximum exceeds the threshold: FWER=P( [ i fT i >ugjH 0 ) =P(T max >ujH 0 ) =1¡F maxjH 0 (t) =1¡(1¡® 0 )=® 0 (4.1) 62 where F maxjH 0 is the cumulative density function of the maximum statistic under the null hypothesis. Therefore, we can control FWER if we choose the thresholdt to be in the(1¡®)100 th percentile of the maximum distribution. Random field (RF) theory methods approximate the upper tail of the maximum dis- tribution F max using the expected value of the Euler characteristic (EC) of the thresh- olded image [WMN + 96, Adl81]. This technique is typically used in fMRI and MEG studies. RF theory relies on several assumptions: that the image has the same parametric distribution at each spatial location, that the point spread function has two derivatives at the origin, that there is sufficient smoothness of the random field, and that there is a sufficiently high threshold for the asymptotic results to be accurate. Most of these assumptions are true in our case, except for image homogeneity. We therefore employ a pre-processing procedure to normalize an inhomogeneous statistic map to a homoge- nous one. In this chapter, we first introduce a Gaussian approximation for the single voxel test statistic. Though the random field theory is not limited to Gaussian distributed statistics, we found that the Gaussian distribution best fits our test statistic. The entire statistic map is modeled as an inhomogeneous Gaussian field. This statistic map is then normalized to a homogenous field with pre-estimated segmentation information. Finally, maximum statistic thresholding with a fixed familywise error rate is applied [WMN + 96]. Digi- tal phantom simulations and a clinical example are used to demonstrate the proposed thresholding algorithm. 63 4.2 Theory 4.2.1 Random Field Theory In this subsection, we review basic concepts and equations for the random field theory. Theoretical results are based on continuous random fields, and under sufficient smooth- ness assumptions, can be applied to discrete statistical maps, such as our lesion detection maps. We present the unified approach, as introduced in [WMN + 96], to assess the prob- ability that the maximum value in the entire statistic map exceeds a given threshold u under the null hypothesis. Mathematical details can be found in [NH03, Adl81]. Consider a continuous Gaussian random fieldZ(s) defined ons2 ½R D , where D is the dimension of the process. LetZ(s) have zero mean and unit variance, as would be required by the null hypothesis. For a threshold value u applied to the field, the suprathreshold regions are known as the ‘excursion set’,A u =\fs:Z(s)>ug. The Euler characteristic (EC) Â(A u )´ Â u is a topological measure of the excursion set. It counts the number of connected suprathreshold regions or ‘clusters’, minus the number of ’holes’ plus the number of ‘hollows’ [NH03]. A closed-form approximation to com- pute the expectation of EC was given in [Adl81], and then generalized in [WMN + 96] with the introduction of resolution elements (RESELs). Under reasonable approximations, the expectation of the EC can be used to control the FWER in discrete statistical mapsT i , wherei is the image pixel index. We require the discrete random fields to be sufficiently smooth, so that they adequately sample their continuous counterparts. For high thresholds, the ‘holes’ and ‘hollows’ disappear and the EC has only ‘clusters’. For even higher thresholds, the EC is either 0 or 1, indicating the presence of a cluster, or equivalently the presence of voxels that exceed a threshold. 64 Therefore, the chance of making at least one false positive under the null hypothesis is given by the expected value of the EC: FWER=P( [ i fT i ¸ugjH 0 ) =P(T max ¸ujH 0 ) ¼E(Â u ) ¼ D X d=0 R d ¢½ d (u) (4.2) where d is the dimensionality index, R d is the d-dimensional RESEL count, a unitless quantity that depends only on topological features of the random field, and½ d (u) is the EC density that depends on the type of statistical field (such asz,t,X 2 , and Hotelling’s T 2 ) [WMN + 96]. The summation of the lower dimensional terms (d < D) compensate for the case the excursion set touches the boundary. For a Gaussian random field we have: ½ 0 (u)= Z 1 u 1 (2¼) 1=2 e ¡t 2 =2 dt ½ 1 (u)= (4ln2) 1=2 2¼ e ¡u 2 =2 ½ 2 (u)= (4ln2) (2¼) 3=2 ue ¡u 2 =2 ½ 3 (u)= (4ln2) 3=2 (2¼) 2 (u 2 ¡1)e ¡u 2 =2 (4.3) The RESEL countR d is given by: R d =¸()j¤j 1=2 (4ln2) ¡d=2 (4.4) 65 where ¸() is the Lesbegue measure of the search region, i.e. the volume V in three dimensions or area S in two dimensions. ¤ is the variance-covariance matrix of the gradient of the random field: ¤, varfrTg. For example, in 3-dimensional space it is defined as: ¤= 2 6 6 6 4 var( @ 2 T @ 2 x ) cov( @T @x ; @T @y ) cov( @T @x ; @T @z ) cov( @T @y ; @T @x ) var( @ 2 T @ 2 y ) cov( @T @y ; @T @z ) cov( @T @z ; @T @x ) cov( @T @z ; @T @y ) var( @ 2 T @ 2 z ) 3 7 7 7 5 (4.5) The square root of its determinant,j¤j 1=2 , is a measure of roughness of the random field. Because the termj¤j 1=2 lacks interpretability, Worsley [WMN + 96] proposed a repa- rameterization based on the convolution of a white noise field with a Gaussian kernel. With proper selection of the FWHM of the Gaussian kernel, the resulting random field can match the smoothness of the original data. Therefore, the FWHM is a representation of the smoothness of the original random field. If the Gaussian kernel has variance matrix§ k , the roughness of the convolved white noise field is ¤ = § ¡1 k =2 [Hol94]. In the 3-dimensional space, if § k has¾ 2 x , ¾ 2 y , ¾ 2 z on the diagonal and zeros elsewhere, then: j¤j 1=2 =j§ k j ¡1=2 ¢2 ¡3=2 =(¾ x ¾ y ¾ z ) ¡1 ¢2 ¡3=2 (4.6) A one-dimensional Gaussian kernel with standard deviation ¾ has FWHM equal to ¾ p 8ln2. By allowing different variance in each dimension, the roughness of the 3- dimensional random field becomes: j¤j 1=2 = (4ln2) 3=2 FWHM x FWHM y FWHM z (4.7) 66 Statistic Map L 2 (y) Normalize to N(0,1) RF Image Segmentation Max Statistic Thresholding for GRF Figure 4.1: The block diagram for thresholding an inhomogeneousL 2 (y) statistic map. A RESEL in three dimensions is defined as a spatial element with extent FWHM x £FWHM y £FWHM z . To compute the RESEL count of the random field with dimensiond=3 and volume¸()=V , we combine equations (4.7) and (4.4): R 3 = V FWHM x FWHM y FWHM z (4.8) To match the smoothness of the original data, [Jen00, KPF + 99] propose algorithms to compute FWHM x , FWHM y and FWHM z . To compute the FWER of a 3-dimensional Gaussian random field at a threshold u, we substitute equations (4.3) and (4.8) into equation (4.2) with D = 3. For simplicity, we ignore alld<3 terms, because they are typically very small for regular random field shapes. FWER 3D = V ¢(u 2 ¡1)e ¡u 2 =2 (4ln2) 3=2 (2¼) 2 ¢ FWHM x ¢ FWHM y ¢ FWHM z (4.9) Similarly, the FWER value for a 2-dimensional field with areaS can be derived as: FWER 2D = S¢ue ¡u 2 =2 (4ln2) 3=2 (2¼) 3=2 FWHM x ¢ FWHM y (4.10) Given a preset FWER, typically5%, we solve equations (4.9) and (4.10) to compute the desirable thresholdu for a Gaussian statistical map. 67 Table 4.1: Micro Parameters for True TACs Used in the Test Statistic Gaussianity Vali- dation Normal Tissue Primary Tumor Metastasis K 1 , ml/g/min 4.5x10 ¡2 3.6x10 ¡2 2.7x10 ¡2 k 2 , /min 3.1x10 ¡1 2.1x10 ¡1 7.2x10 ¡2 k 3 , /min 3.4x10 ¡3 1.2x10 ¡1 2.4x10 ¡2 k 4 , /min 2.0x10 ¡4 1.8x10 ¡3 6.0x10 ¡4 V p , ml/g 1.1x10 ¡1 1.7x10 ¡2 1.6x10 ¡2 4.2.2 Controlling FWER forL 2 (y) Statistic Maps Random field theory relies on several assumptions, including the parametric distribution (usually a Gaussian distribution) at each spatial location, and the homogeneity of the field. However, neither assumption holds for anL 2 (y) statistic map. To the best of our knowledge, the single voxel test statisticL 2 (y) does not have a closed-form probability density function. Furthermore, due to PET image inhomogeneity caused by physiolog- ical differences in different tissue types, the L 2 (y) statistic map is not a homogeneous random field. We address the first issue by introducing a Gaussian approximation, and the second one by introducing a piecewise homogeneity and a normalization step. These two issues are discussed in detail in subsections 4.2.3 and 4.2.4, respectively. The complete thresholding process can be divided into three steps (see Figure 4.1): image segmentation, statistic map normalization and fixed FWER maximum statistic thresholding. In the segmentation step, the image is segmented into a set of nominally homogeneous regions, such that the statistic map in each region is approximately a homogeneous Gaussian random field. The image-based segmentation can be done auto- matically, either by co-registering PET and CT data and using the anatomical informa- tion from the CT image for segmentation [BMGM + 97, HHR01, BGS + 00], or by apply- ing TAC-based segmentation in a dynamic PET image [WFMF02, BGYW03]. For the 68 sake of simplicity, we use TAC-based manual segmentation in the sample data shown in this chapter. The statistic map is spatially downsampled within each segmented region to compute the ensemble mean and variance. Then, the statistic map is normalized with the estimated mean and variance for each region to generate an N(0;1) Gaussian ran- dom field. The detailed normalization procedure is provided in section 4.2.4. Finally, this normalized statistic map is thresholded for a preset FWER using the random field theory described in section 4.2.1. 4.2.3 The Matched Subspace Detection and Gaussian Approxima- tion The random field theory introduced in section 4.2.1 is for Gaussian statistics only. Unfortunately, the test statistic generated by the matched subspace detection is not Gaus- sian, and may not have a known parametric form. In a classical matched subspace detector, the likelihood ratio is defined for the constrained parameter set £ 0 , £ 0 ½ £ over an unconstrained parameter set£ [CB01], and consequently results in anF orÂ 2 distributed statistic depending on the noise model. In the proposed matched subspace detector with known noise covariance, the likelihood ratio is defined for one parame- ter set £ 1 over another parameter set £ 0 [Sch90, Pro00, Fuk90], and results in a non- parametric distribution. In this subsection, the matched subspace detection algorithm is briefly reviewed and its test statistic properties are studied. We provide ideal conditions under which the test statistic can be approximated by a Gaussian random variable. Our simulations verify that the Gaussian approximation is accurate in practice. LetN be the number of frames in a dynamic PET image sequence, andy2R N£1 be the measured TAC at one voxel. A linear subspace model with additive Gaussian noise is used to describe the received TAC. The noise component iny is a zero mean Gaussian random vector with the covariance matrix§2R N£N , which is estimated using a plug 69 estimator applied to the sinogram data [QL00]. The subspace matrices, H 2 R N£p for tumor tissue andS 2 R N£q for normal tissue, are identified from a primary tumor region of interest (ROI) and one or more normal ROIs, respectively. A description of the matched subspace detection theory and the subspace identification algorithm can be found in [LLYL06]. The test statistic produced by this detector under the known noise variance case is: L 2 (y)=y T § ¡1 [H(H T § ¡1 H) ¡1 H T ¡S(S T § ¡1 S) ¡1 S T ]§ ¡1 y (4.11) Define ~ y =§ ¡1=2 y ~ H=§ ¡1=2 H; ~ S=§ ¡1=2 S (4.12) P ~ H = ~ H( ~ H T ~ H) ¡1 ~ H T ; P ~ S = ~ S( ~ S T ~ S) ¡1 ~ S T We can rewriteL 2 (y) as: L 2 (y)= ~ y T P ~ H ~ y¡ ~ y T P ~ S ~ y (4.13) whereP ~ H andP ~ S are projection matrices. The eigen-decomposition ofP ~ H produces: P ~ H =U ~ H ¤ ~ H U T ~ H ; ¤ ~ H =diagf1;1;::1 | {z } ;0;:::;0g with totalr H 1s (4.14) 70 The total number of 1s in¤ ~ H equals the rank of the matrixP ~ H , or equivalently, the rank of the matrixH,r H . Define: z ~ H =U T ~ H ~ y =U T ~ H § ¡1=2 y; m ~ H =U T ~ H § ¡1=2 E(y) (4.15) Given that y is a Gaussian random vector, and z ~ H is a unit variance white Gaussian random vector: z ~ H »N(m ~ H ;I), the first term in (4.13) equals: ~ y T P ~ H ~ y =z T ~ H ¤ H z ~ H = r H X k=1 (z ~ H (k)) 2 (4.16) The summation term P r H k=1 (z ~ H (k)) 2 obeys ar H -degree of freedom non-centralÂ 2 dis- tribution, with the mean and variance [Pro00]: mean:r H +m T ~ H m ~ H variance:2(r H +2m T ~ H m ~ H ) (4.17) Similarly, we have: P ~ S =U ~ S ¤ ~ S U T ~ S ; ¤ ~ S =diagf1;1;::1 | {z } ;0;:::;0g with totalr S 1s z ~ S =U T ~ S § ¡1=2 y; m ~ S =U T ~ S § ¡1=2 E(y) (4.18) ~ y T P ~ S ~ y = r S X l=1 (z ~ S (l)) 2 71 The summation P r S l=1 (z ~ S (l)) 2 obeys another non-central Â 2 distribution with an r S - degree of freedom, where r S is the rank of the matrixS. Substituting equation (4.16) and (4.19) into equation (4.13) results in: L 2 (y)= r H X k=1 z 2 ~ H(k) ¡ r S X l=1 z 2 ~ S(l) (4.19) When r H and r S are large numbers, by applying the Lindeberg central limit theorem [CB01], the two summation terms are approximately Gaussian random variables. Fur- thermore, ifP ~ H ? P ~ S , the two summation terms are independent. ThenL 2 (y) can be approximated as a Gaussian random variable with mean r H ¡r S +m T ~ H m ~ H ¡m T ~ S m ~ S and variance2(r H +r S )+4(m T ~ H m ~ H +m T ~ S m ~ S ). Thatr H andr S are large numbers and P ~ H ?P ~ S are the ideal conditions under which theL 2 (y) statistic can be approximated using a Gaussian random variable. The orthogonality and the matrices rank numbers (r H andr S ) are dependent on the ROIs from which they are identified. Usually, these two conditions do not strictly hold in practice. For example, in the clinical cancer studies we conducted in [LLYL06], the rank numbers were between 3 and 5, and orthogonality did not hold. To solve the problem, we relax these assumptions and assume thatL 2 (y) still approximately obeys a Gaussian distribution. We then investigate this Gaussian approximation using a Monte Carlo simulation. Two digital phantoms (not shown here), one with a primary tumor and one with a metastasis, were generated by inserting a primary tumor and a metastasis region into a flat background, respectively. To simulate real physiological processes, TACs assigned in the digital phantom were obtained from a clinical dynamic FDG-PET breast cancer study. ROIs were drawn in the clinical images to define background, primary tumor and metastatic regions. The micro-parameters in each type of tissue were estimated by fitting 72 −25 −20 −15 −10 −5 0 5 10 0 0.05 0.1 0.15 0.2 0.25 Test Statistic Value PDF H 0 : Histogram H 0 : Fitted Gaussian H 1 : Histogram H 1 : Fitted Gaussian Figure 4.2: The histogram and fitted Gaussian PDF of a single voxelL 2 (y) test statistic for both normal and tumor tissue. The perfect match between the histogram and the fitted Gaussian PDF implies an accurate Gaussian approximation. −4 −3 −2 −1 0 1 2 3 4 −15 −10 −5 0 5 10 Standard Normal Quantiles Quantiles of Input Sample Statistic for Normal Tissue Statistic for Tumor . Figure 4.3: Quantile-quantile plots of the single voxelL 2 (y) test statistic versus a stan- dard Gaussian random variable. The linearity of the quantile-quantile plot validates the Gaussian approximation for theL 2 (y) test statistic. to a three compartment model [HPH + 80] using the Newton-Raphson procedure [LH74]. The arterial blood activity curve was approximated using an ROI in the left ventricle [GSH + 89]. The resulting micro-parameters are listed in Tab. 4.1. To simulate the real detection procedure, subspace matricesH andS were identified from the ROIs drawn from the reconstructed primary tumor phantom. TheL 2 (y) test statistic for the normal tissue and the tumor were computed by applying the L 2 (y) detector to the TAC at one background voxel, and one metastasis voxel in the reconstructed image, respectively. 1200 simulated test statistic data were generated for each case. The histograms for the L 2 (y)test statistic for both background and tumor are shown in Figure 4.2. The 73 0 20 40 60 80 100 120 140 160 0 50 100 150 200 250 X−axis Activity True Noiseless Reconstruction Reconstruction 0 20 40 60 80 100 120 140 160 −60 −40 −20 0 20 40 60 80 X−axis Activity Noiseless Reconstruction Reconstruction Figure 4.4: The profiles of a phantom image before and after removing the ensemble mean in each segmentation region. The upper plot shows the profile before removing the ensemble mean. The lower plot shows the profile after removing the ensemble mean. A partial volume effect causes peaks on the boundaries in the mean-value-removed image. These peaks can induce an overestimated ensemble variance. Gaussian PDF function with the ensemble mean and the ensemble variance estimated from the data samples was overlapped on the histogram for each case. The perfect match between the histogram and the Gaussian PDF validates the Gaussian approximation. To further verify the Gaussian approximation, a quantile-quantile plot of 1200 simulated data against standard Gaussian distribution is shown in Figure 4.3. The linearity of plot shows that the Gaussian approximation is accurate. Therefore, we can model the single voxelL 2 (y) test statistic as a Gaussian random variable, and theL 2 (y) statistic map as a Gaussian random field. 74 −40 −20 0 20 0 0.1 0.2 0.3 0.4 0.5 −10 −5 0 5 0 0.05 0.1 0.15 0.2 −400 −200 0 200 0 0.2 0.4 0.6 0.8 −80 −60 −40 −20 0 0 0.05 0.1 0.15 0.2 −400 −300 −200 −100 0 0 0.02 0.04 0.06 0.08 0.1 0.12 −400 −200 0 200 0 0.05 0.1 0.15 0.2 Region I Region II Region III Percentage Value Value Region I Region II Region III Figure 4.5: Histograms of the spatially sampledL 2 (y) test statistic values in three seg- mentation regions (Figure 4.7(b)) before (Left column), and after (Right column) remov- ing the outliers. The fitted curves are the Gaussian PDF with the ensemble mean and variance of each region. 4.2.4 Normalization of the Non-stationary Gaussian Random Field We address two issues in the statistic map normalization step: partial volume effect and outliers in the statistical distribution. Partial volume effects are artifacts that occur where multiple tissue types contribute to a single voxel, resulting in a blurring across bound- aries. For any segmentation region, even if the segmentation itself is correct, removing the ensemble mean of the region will result in a peak near the boundary, which affects the smoothness of the mean-value-removed image. These peaks also cause an overes- timated ensemble variance. Figure 4.4 illustrates the image profiles and demonstrates how the partial volume effect affects smoothness and variance. For this reason, a three voxel-width boundary is excluded from the data samples, to estimate the variance in the 75 statistic map. However, the excluded boundary areas are not left blank in the normalized field, because leaving blank boundaries will change the smoothness of the random field. Instead, we replaced the boundary values with neighbor voxels values in the normalized random field. The second issue is that outliers in the statistical distribution affect variance estima- tion. Although the image is segmented into nominally homogeneous regions, the homo- geneity within each region is only held to some extent. The heterogenous voxels inside the region produce outliers in the test statistical distribution, which can cause an over- estimated variance. We use an interquartile range based criterion to detect and remove outliers from the data samples [Job99]. The interquartile range (IQR) is a measure of statistical dispersion, being equal to the difference between the third and first quartiles. The data sample value beyond 5£ IQR from the median value of all data samples is considered an outlier, and is excluded from the data samples used in variance estima- tion. Figure 4.5 shows histograms of the L 2 (y) test statistic in the three segmentation regions, before and after removing outliers. 4.3 Simulations, Experiments and Results 4.3.1 Digital Phantom Simulation A digital phantom Monte Carlo simulation was conducted to verify the proposed thresh- olding procedure. We simulated a 2D scanner with the same geometry as the Siemens Focus 220 small animal scanner. The reconstructed image size was 128£ 128, and the sinogram dimension was 288£ 252. Two digital phantoms, one with a primary tumor and one with three metastases, were generated. To simulate clinical background images, we used a tumor-free clinical dynamic PET image (Figure 4.6, 1st row, left) as a true normal background template. This image, which was scanned using a Siemens 76 Clinic Image Background Clinic Image Tumors Primary Tumor Metastasis True Phantom Primary Tumor True Phantom Metastases Reconstruction Primary Tumor Reconstruction Metastases Figure 4.6: Clinical images, true phantoms and reconstructed images. On the 1st row, a clinical tumor free image (left), a clinical image with one primary tumor and one metastasis (right). On the 2nd row, true phantom images. On the 3rd row, reconstructed images. ECAT-953 system, has a total of 28 frames for a 40-minute scan duration (15sec£8, 30sec£4, 1min£11, 5min£5; total 40 min, 28 frames). To simulate clinical tumor activity, we used a dynamic PET image from a lung cancer patient with metastatic dis- ease as the tumor template (Figure 4.6, 1st row, right). By inserting the primary tumor drawn from the tumor template into the background template, a primary tumor phantom was generated (Figure 4.6, 2nd row, left). Due to the partial volume effect, TACs in the small size metastasis are corrupted by the background TACs. Therefore, we did not directly insert the clinic metastasis region to generate the metastasis phantom. Instead, this phantom image was generated by inserting three small flat regions with averaged 77 a b c d I II III Figure 4.7: Simulation results for the thresholding algorithm. (a) The L 2 (y) statistic map, (b) Three segmentation regions, (c) The normalizedN(0;1) GRF, (d) The binary decision map. clinical metastatic TAC into the background template (Figure 4.6, 2nd row, right). All three phantoms were scaled by the same factor, such that the total count at the last frame of the background phantom was set to one million. Each frame was forward projected and reconstructed using the MAP algorithm [QLH + 98], initialized with two iterations of OSEM with 9 subsets followed by 20 iterations of the MAP PCG algorithm. The regularization parameter was set to 0.1. The clinical images, true phantom images and reconstructed images are shown in Figure 4.6. To approximate a real detection proce- dure, the subspace matrices H and S were identified from the reconstructed primary tumor image [LLYL06]. The proposed three-step thresholding algorithm with a preset FWER=5% was applied to theL 2 (y) statistical map computed using a matched subspace detector. TheL 2 (y) statistic map, the segmentation result, the normalized statistic map and the binary decision map are all shown in Figure 4.7. 78 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 Preset FWER FWER Simulation Result Figure 4.8: Simulated FWER against preset FWER. The slope is less than one, which implies a conservative thresholding. Primary Tumor Metastasis 1 Metastasis 2 False Alarm Image N(0,1) GRF Binary Decision Primary Tumor Metastasis 1 Metastasis 2 False alarm Transverse Coronal Sagittal Original Images Figure 4.9: A clinical breast cancer study. The 1st column shows the transverse, sagittal and coronal views of the last dynamic frame. The 2nd»5th columns show thresholding results in the primary tumor, metastasis 1, metastasis 2 and a false alarm, respectively. The 1st row shows the original image. The 2nd row shows the normalized N(0,1) Gaus- sian random field. The 3rd row shows the binary decision map. The results show that the thresholding algorithm detects all three tumors but also generates one false alarm out of 31 planes. 79 To compute the simulated FWER with the preset FWER=5%, one thousand pseudo- random dynamic Poisson data sets were generated for the phantoms. In each recon- structed image, a matched subspace detector followed by the proposed thresholding algorithm were applied. A total of one thousand binary decision maps were generated. The true positive rate wasP D = 97:1%, with FWER=1.8%. When changing the preset FWER from 1% to 30%, the simulated FWER values are plotted in Figure 4.8. The ideal curve should be a straight line with a slope that equals one. In fact, the slope is less than one, which implies a conservative thresholding. Because of partial volume effects and normal tissue inhomogeneity, we expect to overestimate the roughness of the statistic map, which may explain the conservative threshold. 4.3.2 A Clinical Experiment We conclude with an example of this thresholding algorithm applied to the case of a breast cancer patient with metastatic disease. This data was obtained using a Siemens ECAT-953 scanner with a dynamic data sampling protocol (15sec£8, 30sec£4, 1min£21, 5min£3; total 40min, 36 frames). The reconstructed image size was 128£128 with 31 planes. The transverse, sagittal and coronal views of the last frame are shown in the first column of Figure 4.9. Following the procedure described in [LLYL06], we applied the matched subspace detector to generate anL 2 (y) statistical map. The 3D image is manually segmented into three nominally homogenous regions based on TACs: the chest, lung and heart. The thresholding was applied plane-by-plane to theL 2 (y) sta- tistical map using the algorithm described in section 4.2.2. The thresholding algorithm detected all three tumors, including one primary tumor and two metastases, but also gen- erated one false alarm. Figure 4.9 shows the original images, the normalized Gaussian random field and the binary decision maps. To facilitate comparison, only the center plane for each tumor is shown in the figure. 80 4.4 Discussion and Conclusion Using random field theory, we provide a method to threshold the dynamic PET matched subspace detection statistical map. This map is modeled as an inhomogeneous Gaus- sian random field. A three-step thresholding algorithm is presented to achieve a fixed FWER in lesion detection. The digital phantom simulations demonstrate that with a preset FWER=5%, the proposed thresholding algorithm achieves a 97.1% positive rate and 1.8% FWER. Random field theory relies on several assumptions, including that the statistical map has sufficient smoothness and that the voxels have the same parametric distribution at each spatial location. Because of partial volume effects and normal tissue inhomogeneity, we overestimate the roughness of the image. This results in a conserva- tive threshold and a 1.8% error rate which is less than nominal. A conservative threshold has less detection power, but does not give inflated false positives, which is a desirable feature in the detection of lesions. 81 Chapter 5 Conclusion and Future Work 5.1 Conclusion This dissertation has focused on three issues of computer-aided lesion detection in dynamic PET: feature extraction, detector design and thresholding procedure. The features of the time activity curve (TAC) were studied for the lesion detection task. Starting from the kinetic model and the observed temporal domain differences, the frequency domain differences between the tumor and normal TACs were studied. Three feature extraction algorithms were proposed to identify features from the training images with pre-defined primary tumor and normal tissue regions of interest (ROIs). Clinical data sets confirmed that the identified features, i.e. two linear subspaces, were accurate and separable. We theorized that a matched subspace detection algorithm could be used to assist in the detection of small tumors in dynamic PET. A linear subspace model and a gener- alized likelihood ratio test were used to derive the detection algorithm. This algorithm uses the linear subspaces identified from the feature extraction process as a prior knowl- edge and differentiates tumors from background tissue based on the TAC. Three noise models, namely unknown variance white noise, known variance colored noise and spa- tial correlated noise with known covariances, were considered. Applying the detection algorithm on a voxel by voxel basis throughout the dynamic image produces a test sta- tistical image or ’map’ which on thresholding indicates potential locations of secondary or metastatic tumors. According to our simulation and clinical case study, even without 82 using the blood input function, the proposed algorithm is superior to both the Patlak algorithm and the static PET Hotelling observer in the lesion detection task. This result also provides support for the argument that lesion detection may be enhanced by using dynamic information. To make the binary decision on the presence or absence of a lesion at each voxel based on the statistical map, a constant familywise error rate (FWER) thresholding algorithm was proposed. The statistical map was modeled as an inhomogeneous Gaus- sian random field which was validated using a simulation. The thresholding proce- dure includes three steps. First, the PET image is segmented into several homogeneous regions. Second, the statistical map is normalized to a zero mean unit variance Gaussian random field. Finally, a maximum statistic thresholding with a fixed FWER is applied. The proposed thresholding method was demonstrated using both digital phantom and clinical images. Simulations show that with a preset FWER of 5%, the proposed thresh- olding method achieves a 97.1% true positive rate with a 1.8% FWER. The conservative thresholding result is caused by the inhomogeneity in normal tissue. We also show an example of the application of our proposed approach to clinical PET data from a breast cancer patient with metastatic disease. All publications related to this work are listed below: ² Zheng Li, Quanzheng Li, Xiaoli Yu, Peter S. Conti and Richard M. Leahy, “Lesion Detection in Dynamic FDG-PET Using Matched Subspace Method”, submitted to IEEE Trans. Med. Imag., under revision. ² Zheng Li, Quanzheng Li, Dimitrios Pantazis and Richard M. Leahy, “Controlling Familywise Error Rate for Matched Subspace Detection in Dynamic PET”, to be sub- mitted to IEEE Trans. Med. Imag. 83 ² Zheng Li, Quanzheng Li, Xiaoli Yu, Peter S. Conti and Richard M. Leahy, “Perfor- mance of Matched Subspace Detectors for Dynamic FDG PET”, IEEE Nuclear Science Symposium, Medical Imaging Conference (NSS-MIC), 2006. ² Zheng Li, Quanzheng Li, Xiaoli Yu, Peter S. Conti and Richard M. Leahy, “Matched Subspace Detection for Dynamic PET: An ROC Phantom Study for MAP Reconstruc- tion”, IEEE International Symposium on Biomedical Imaging, Symposium Proceedings pp. 1332-1335, Apr., 2006 ² Zheng Li and X. Yu, “Computer Observer Design for Dynamic PET Image Using Time Domain Filtering”, IEEE NSS-MIC, Conference Record pp. 1700-1703, Oct., 2005. ² Zheng Li and X. Yu, “Exploring Physiological Differences in Time-Frequency Domain to Improve Tumor Detectability for Dynamic PET”, IEEE NSS-MIC, Oct., 2004. ² Zheng Li and Xiaoli Yu, “Exploring Frequency Differences of Physiological Processes to Enhance Dynamic FDG-PET without Blood Function”, IEEE NSS-MIC, V ol. 4, pp. 2844-2848, Oct., 2003. ² X. Yu, Z. Li, H. Jadvar and P. Conti. “Assessment of Noninvasive Blood Time Activity Extraction for Dynamic Oncology”, J. Nucl. Med., 43:102P, May 2002 5.2 Future Work 5.2.1 Lesion Detection Experiment We have demonstrated that the proposed algorithm improves the detection of small metastatic lesions in small animal studies and in one clinical sample. We plan to con- tinue our small animal study and extend our clinical studies. In a small animal longitudi- nal study, the selected human cancer cells will be injected into mice. Dynamic PET, CT 84 and optical scans will be carried out using microPET, microCAT and Xenogen biolumi- nescence systems to monitor the tumor growth. The detection results generated by the matched subspace detection algorithm on earlier PET scan data will be confirmed using later PET, CT and optical scans. The aim of this study is to determine how early the tumor can be detected using the proposed lesion detection algorithm and as a precursor to studies in humans. In the clinical part of our study, we plan to collect breast cancer patient dynamic FDG-PET data using a high resolution PET-CT scanner . In comparison to small animal studies, it is not possible to monitor tumor growth over a long periods in clinical exper- iments. We plan to apply the proposed detection algorithm on dynamic data to detect small tumors and verify the detection results from the CT-PET combined image by an experienced doctor. The goal of this study is to use the detection algorithm in more clinical trials to test its efficacy. 5.2.2 Detection Algorithm The maximum likelihood estimator in the current matched subspace detector does not have positive constraints on the basis function coefficients. Although the unconstrained estimator is correct for PCA basis functions, in the cases where physiological process functions are used as basis functions, positive constraints are needed. In this case, we need to either find an analytical solution, or use the iterative solution for the maximum likelihood estimator. The statistical properties of the detector output and the threshold- ing algorithm would need to be changed accordingly. In practice, using a single detector to make a decision does not usually guarantee a high detection rate and low false positive rate. One solution is to combine the outputs from multiple detectors into one decision. Some research has studied the support vector machine (SVM) for lesion detection in dynamic PET. Preliminary results show that the 85 SVM has potential in detecting lesions, but also has some drawbacks. In comparison to the matched subspace detector, the SVM is more sensitive to amplitude information. As a result, the SVM is less prone to generate artifacts near the body boundary, but is prone to miss metastases within a high intensity background. Combining the outputs from the matched subspace detector and the SVM into one decision may address the disadvantage of both algorithms. 86 References [Adl81] R.J. Adler. The Geometry of Random Fields. Wiley, 1981. [AKKK99] B. A. Ardekani, J. Kershaw, K. Kashikura, and I. Kanno. Activation detection in functional MRI using subspace modeling and maximum like- lihood estimation. IEEE Trans. Med. Imag., 18(2):101–114, 1999. [Bar80] D. Barber. The use of principal components in the quantitative analysis of gamma camera dynamic studies. Phys. Med. Biol., 25:283–292, 1980. [BGG + 92] H. H. Barrett, T. Gooley, K. Girodias, J. Rolland, T. White, and J. Yao. Linear discriminants and image quality. Image and Vision Computing, 10(6):451–460, 1992. [BGS + 00] M. S. Brown, J. G. Goldin, R. D. Suh, M. F. McNitt-Gray, L. E. Greaser, A. Sapra, K.-T. Li, J. W. Sayre, K. Martin, and D. R. Aberle. Knowledge- based segmentation of thoracic computed tomography images for assess- ment of split lung function. Med. Phys., 27(3):592–598, 2000. [BGYW03] J. G. Brankov, N. P. Galatsanos, Y . Yang, and M. N. Wernick. Segmenta- tion of dynamic PET or fMRI images based on a similarity metric. IEEE Trans. Nucl. Sci., 50(5):1410–1414, 2003. [BMGM + 97] M. S. Brown, M. F. McNitt-Gray, N. J. Mankovich, J. G. Goldin, J. Hiller, L. S. Wilson, and D. R. Aberie. Method for segmenting chest CT image data using an anatomicalmodel: preliminary results. IEEE Trans. Med. Imag., 16:828–839, 1997. [BYRM93] H. H. Barrett, J. Yao, J. P. Rolland, and K. J. Myers. Model observers for assessment of image quality. Proc. Natl. Acad. Sci. USA, 90:9758–9765, 1993. [Car93] R. E. Carson. An approximation formula for the variance of PET region- of-interest values. IEEE Trans. Med. Imag., 12(2):240–250, 1993. 87 [CB01] George Casella and Roger L. Berger. Statistical Inference. Duxbury Press, 2nd edition, 2001. [CFC03] Z. Chen, D. Feng, and W. Cai. Pathological lesion detection in 3D dynamic PET images using asymmetry. In Proceedings of the 12th ICIAP. IEEE Computer Society, 2003. [CJ93] V . J. Cunningham and T. Jones. Spectral analysis of dynamic PET stud- ies. J. Cereb. Blood Flow Metab., 13:15–23, 1993. [CJS93] Z. H. Cho, J. Jones, and M. Singh. Foundations of Medical Imaging. John Wiley & Sons Inc., 1993. [Con96] P. S. Conti. PET and [ 18 F]-FDG in oncology: A clinical update. Nuclear Med. and Bio., 23:717–735, 1996. [CP96] S. R. Cherry and M. E. Phelps. Brain Mapping: The Methods, chapter 8. Academic Press, 1996. [DCGJ99] C. Dykstra, A. Celler, K. Greer, and R. Jaszczak. The use of image morphing to improve the detection of tumors in emission imaging. IEEE Trans. Nucl. Sci., 46:673–679, 1999. [Doi07] Kunio Doi. Computer-aided diagnosis in medical imaging: Historical review, current status and future potential. Comp. Med. Imag. and Graph- ics, 31:198–211, 2007. [FAH + 96] K. J. Friston, J. Ashburner, A. Holmes, J. B. Poline, R. Turner, R. Wise, C. Price, R. Dolan, P. Fletcher, J. Coull, and C. Frith. Statistical paramet- ric mapping and functional neuroimaging. Technical report, Functional Imaging Laboratory, May 1996. [FR96] J. A. Fessler and W. L. Rogers. Spatial resolution properties of penalized- likelihood image reconstruction: Space-invariant tomographs. IEEE Trans. Med. Imag., 5(9):1346–1358, Sep 1996. [Fuk90] Keinosuke Fukunaga. Introduction to Statistical Pattern Recognition. Academic Press, 2nd edition, 1990. [FW93] D. Feng and X. Wang. A computer simulation study on the effects of input function measurement noise in tracer kinetic modeling with positron emission tomography. Computer Biology Medicine, 23(1):57– 68, 1993. 88 [FWF + 94] K. J. Friston, K. J. Worsley, R. S. J. Frackowiak, J. C. Mazziotta, and A. C. Evans. Assessing the significance of focal activations using their spatial extent. Human Brain Mapping, 1:210–220, 1994. [GKA01] M. L. Giger, N. Karssemeijer, and S. G. Armato. Guest Editorial computer-aided diagnosis in medical imaging. IEEE Trans. Med. Imag., 20:1205–1208, 2001. [GSH + 89] S. S. Gambhir, M. Schwaiger, S. C. Huang, J. Krivokapich, H. R. Schel- bert, C. A. Nienaber, and M. E. Phelps. Simple noninvasive quantification method for measuring myocardial glucose utilization in humans employ- ing positron emission tomography and fluorine-18 deoxyglucose. J. Nucl. Med., 30(3):359–366, 1989. [HCH + 83] S. C. Huang, R. E. Carson, E. J. Hoffman, J. Carson, N. MacDonald, J. R. Barrio, and M. E. Phelps. Quantitative measurement of local cerebral blood flow in humans by positron emission tomography and 15o-water. J. Cereb. Blood Flow Metab., 3:141–153, 1983. [HHR01] S. Hu, E. A. Hoffman, and J. M. Reinhardt. Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images. IEEE Trans. Med. Imag., 20(6):490–498, 2001. [HL94] H. M. Hudson and R. S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag., 13:100–108, 1994. [Hol94] A. P. Holmes. Statistical issues in functional brain mapping. PhD thesis, University of Glasgow, Glasgow, 1994. [HPH + 80] S. C. Huang, M. E. Phelps, E. J. Hoffman, K. Sideris, C. J. Selin, and D. E. Kuhl. Noninvasive determination of local cerebral metabolic rate of glucose in man. Am. J. Physiol. Endocrinol. Metab., 238:E69–E82, 1980. [HT87] Y . Hochberg and A.C. Tamhane. Multiple Comparison Procedures. Wiley, 1987. [Hue84] R. H. Huesman. A new fast algorithm for the evaluation of regions of interest and statistical uncertainty in computed tomography. Phys. Med. Biol., 29(5):543–552, 1984. [HY99] C. C. Huang and X. Yu. A new method of computer-aided feature iden- tification for lesion detection in PET-FDG dynamic study. In Nuclear Science Symposium, Conference Record, volume 3, pages 1277–1281, 1999. 89 [HYB + 97] C. C. Huang, X. Yu, J.R. Bading, , and P.S. Conti. Computer-aided lesion detection with statistical model-based features in PET images. IEEE Trans. Nucl. Sci., 44(6):2509–2521, 1997. [Jen00] Mark Jenkinson. Estimation of smoothness from the residual field. Tech- nical report, Oxford Centre for Functional Magnetic Resonance Imaging of the Brain (FMRIB), 2000. [Job99] J. D. Jobson. Applied Multivariate Data Analysis: Regression and Exper- imental Design, volume 1. Springer, 1999. [JWA + 95] T. Jansson, J. E. Westlin, H. Ahlstrom, A. Lilja, B. Langstrom, and J. Bergh. Positron emission tomography studies in patients with locally advanced and/or metastatic breast cancer: a method for early therapy evaluation? Journal of Clinical Oncology, 13:1470–1477, 1995. [KPF + 99] S. J. Kiebel, J.-B. Poline, K. J. Friston, A. P. Holmes, and K. J. Worsley. Robust smoothness estimation in statistical parametric maps using stan- dardized residuals from the general linear model. NeuroImage, 10:756– 766, 1999. [KS87] A. C. Kak and M. Slaney. Principles of Computerized Tomography Imag- ing. IEEE Press, 1987. [KV AN + 05] L. M. Kenny, D. M. Vigushin, A. Al-Nahhas, S. Osman, S. K. Luthra, S. Shousha, R. C. Coombes, and E. O. Aboagye. Quantification of cellu- lar proliferation in tumor and normal tissues of patients with breast cancer by [ 1 8f]fluorothymidine-positron emission tomography imaging: evalu- ation of analytical methods. Cancer Research, 65:10104–10112, Nov. 2005. [LAAL07] Q. Li, E. Asma, S. Ahn, and R. M. Leahy. A fast fully 4D incremental gradient reconstruction algorithm for list mode PET data. IEEE Trans. Med. Imag., 26(1):58–67, 2007. [Lan96] N. Lange. Tutorial in biostatistics, statistical approaches to human brain mapping by fMRI. Statistics in Med., 15:389–428, 1996. [LAQ + 04] Q. Li, E. Asma, J. Qi, J. R. Bading, and R. M. Leahy. Accurate estima- tion of the Fisher information matrix for the PET image reconstruction problem. IEEE Trans. Med. Imag., 23(9):1057–1064, 2004. [LH74] C. L. Lawson and R.J. Hanson. Solving Least Squares Problems, chap- ter 23, page 161. Prentice-Hall, 1974. 90 [LLYL06] Zheng Li, Quanzheng Li, Xiaoli Yu, and Richard Leahy. Matched sub- space detection for dynamic PET: an ROC phantom study for map recon- struction. In IEEE ISBI, Symposium Proceedings, pages 1332–1335, 2006. [Log00] Jean Logan. Graphical analysis of PET data applied to reversible and irreversible tracers. Nucl. Med. and Bio., 27:661–670, 2000. [LQ00] R. M. Leahy and J. Qi. Statistical approaches in quantitative positron emission tomography. Statistics and Computing, 10:147–165, 2000. [LSM + 93] G. Lucignani, K. C. Schmidt, R. M. Moresco, G. Striano, F. Colombo, L. Sokoloff, and F. Fazio. Measurement of regional cerebral glucose uti- lization with fluorine-18-FDG and PET in heterogeneous tissues: Theo- retical considerations and practical procedure. J. Nucl. Med., 34(3):360– 369, 1993. [MMC + 96] S. R. Meikle, J. C. Matthews, V . J. Cunningham, D. L. Bailey, L. Livier- atos, T. Jones, and P. Price. Spectral analysis of PET projection data. In IEEE Nuclear Science Symposium, Conference Record., volume 3, pages 1888–1892, 1996. [NH03] T.E. Nichols and S. Hayasaka. Controlling the familywise error rate in functional neuroimaging: A comparative review. Statistical Methods in Medical Research, 12(5):419–446, 2003. [NQAL02] T. E. Nichols, J. Qi, E. Asma, and R. M. Leahy. Spatiotemporal recon- struction of list-mode PET data. IEEE Trans. Med. Imag, 21(4):396–404, 2002. [PB85] C. Patlak and R. G. Blasberg. Graphical evaluation of blood-to-brain transfer constants from multiple-time uptake data. generalizations. J. Cereb. Blood Flow Metab., 5:584–590, 1985. [PBBL94] F. Pedersen, M. Bergstrome, E. Bengtsson, and B. Langstrom. Principal component analysis of dynamic positron emission tomography images. European J. of Nucl. Med., 21:1285–1292, 1994. [PBF83] C. Patlak, R. G. Blasberg, and J. D. Fenstermacher. Graphical evaluation of blood-to-brain transfer constants from multiple-time uptake data. J. Cereb. Blood Flow Metab., 3:1–7, 1983. [PMS86] M. E. Phelps, J. C. Mazziotta, and H. R. Schelbert. Positron Emis- sion Tomography and Autoradiology, Principles and Applications for the Brain and Heart. Raven Press, 1986. 91 [PNBL05] D. Pantazis, T. Nichols, S. Baillet, and R. M. Leahy. A comparison of random field theory and permutation methods for the statistical analysis of MEG data. Neuroimage, 25(2):383–394, 2005. [Pro00] John G. Proakis. Digital Communications. McGraw Hill, 4th edition, 2000. [QH01] J. Qi and R. H. Huesman. Theoretical study of lesion detectability of MAP reconstruction using computer observers. IEEE Trans. Med. Imag., 20(8):815–822, Aug. 2001. [QL00] J. Qi and R. M. Leahy. Resolution and noise properties of MAP recon- structions in fully 3D PET. IEEE Trans. Med. Imag., 19(5):493–506, May 2000. [QLC + 98] 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–1013, 1998. [QLH + 98] 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, 1998. [SC91] L. G. Strauss and P. S. Conti. The applications of PET in clinical oncol- ogy. J. Nucl. Med., 32(4):623–648, 1991. [Sch90] L. L. Scharf. Statistical Signal Processing - Detection, Estimation, and Time Series Analysis. Addison Wesley, 1990. [SdD04] J. Sijbers and A.J. den Dekker. Generalized likelihood ratio tests for com- plex fMRI data. In Medical Imaging 2004, Proceeding of SPIE, volume 5369, pages 652–663, 2004. [SdD05] J. Sijbers and A. J. den Dekker. Generalized likelihood ratio tests for complex fMRI data: a simulation study. IEEE Trans. Med. Imag., 24(5):604–611, 2005. [SF94] L. L. Scharf and B. Friedlander. Matched subspace detectors. IEEE Trans. Signal Processing, 42(8):2146–2157, Aug. 1994. [SKP + 07] L. G. Strauss, S. Klippel, L. Pan, K. Schonleben, U. Haberkorn, and A. Dimitrakopoulou-Strauss. Assessment of quantitative FDG PET data in primary colorectal tumours: which parameters are important with respect to tumour detection. J. Nucl. Med., 34:868–877, 2007. 92 [SL74] L. Shepp and B. Logan. The Fourier reconstruction of a head section. IEEE Trans. Nucl. Sci., 21:21–33, 1974. [SLDV04] L. Schrevens, N. Lorent, C. Dooms, and J. Vansteenkiste. The role of PET scan in diagnosis, staging, and management of non-small cell lung cancer. The Oncologist, 9:633–643, 2004. [SMS91] K. Schmidt, G. Mies, and L. Sokoloff. Model of kinetic behavior of deoxyglucose in heterogeneous tissues in brain: a reinterpretation of the significance of parameters fitted to homogeneous tissue models. J. Cereb. Blood Flow Metab., 11:10–24, 1991. [Sul93] F. O. Sullivan. Imaging radiotracer model parameters in PET: a mixture analysis approach. IEEE Trans. Med. Imag., 12:399–412, 1993. [SV82] L. Shepp and Y . Vardi. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imag., 1(2):113–122, 1982. [TML + 94] F. Turkheimer, R. M. Moresco, G. Lucignani, L. Sokoloff, F. Fazio, and K. Schmidt. The use of spectral analysis to determine regional cerebral glucose utilization with positron emission tomography and [ 18 F]fluorodeoxyglucose: theory, implementation, and optimization pro- cedures. J. Cereb. Blood Flow Metab., 14:406–422, 1994. [TSDS + 03] T. Thireou, L. G. Straussb, A. Dimitrakopoulou-Straussb, G. Kontax- akisc, S. Pavlopoulosa, and A. Santos. Performance evaluation of prin- cipal component analysis in dynamic FDG-PET studies of recurrent col- orectal cancer. Computerized Medical Imaging and Graphics, 27:43–51, 2003. [VK95] Martin Vetterli and Jelena Kovacevic. Wavelets and subband coding. Prentice Hall, 1995. [VO91] M Viberg and B Ottersten. Sensor array processing based on subspace fitting. IEEE Trans. Signal Processing, 39(5):1110–1121, May 1991. [WEMN92] K. J. Worsley, A. C. Evans, S. Marrett, and P. Neelin. A three- dimensional statistical analysis for CBF activation studies in human brain. Journal of Cerebral Blood Flow and Metabolism, 12:900–918, 1992. [WFMF02] K. P. Wong, D. Feng, S. R. Meikle, and M. J. Fulham. Segmentation of dynamic PET images using cluster analysis. IEEE Trans. Nucl. Sci., 49(1):200–207, 2002. 93 [WMN + 96] K.J. Worsley, S. Marrett, P. Neelin, A.C. Vandal, K.J. Friston, and A.C. Evans. A unified statistical approach for determining significant signals in images of cerebral activation. Brain Mapping, 4(58-73), 1996. [WZH + 93] R. L. Wahl, K. Zasadny, M. Helvie, G. D. Hutchins, B. Weber, and R. Cody. Metabolic monitoring of breast cancer chemohormonotherapy using positron emission tomography: initial evaluation. Journal of Clin- ical Oncology, 11(2101-2111), 1993. [YZS + 04] H. Ying, F. Zhou, A. F. Shields, O. Muzik, D. Wu, and E. I. Heath. A novel computerized approach to enhancing lung tumor detection in whole-body PET images. In Proc. of 26th Annual International Conf. IEEE EMBS, Sep. 2004. [ZPL + 01] H. Zhuang, M. Pourdehnad, E. S. Lambright, A. J. Yamamoto, M. Lanuti, P. Li, P. D. Mozley, M. D. Rossman, S. M. Albelda, and A. Alavi. Dual time point 18F-FDG PET imaging for differentiating malignant from inflammatory. J. Nucl. Med., 42:1412–1417, 2001. [ZW96] K. R. Zasadny and R. L. Wahl. Enhanced FDG-PET tumor imaging with correlation-coefficient filtered influx-constant images. J. Nucl. Med., 37(2), 1996. * 94 Appendix A Physiological process function In this appendix, we derive the conditions under which the physiological process func- tion has a positive slope. The conclusions are used for the model based feature extraction in section 2.3.2. We assume that all functions in this proof are continuous and differentiable. Assume that the blood input functionC p (t) satisfies the following conditions: (i)C p (t) 8 > > > > > < > > > > > : =0; fort<0 ¸0; fort¸0 ; (ii) dC p (t) dt 8 > > > > > < > > > > > : >0; fort<T 0 =0; fort=T 0 <0; fort>T 0 (A.1) The physiological process function is defined as: f(¯;t),C p (t)e ¡¯t = Z t 0 C p (¿)e ¡¯(t¡¿) d¿ = Z t 0 C p (t¡¿)e ¡¯¿ d¿; for ¯¸0; t¸0 (A.2) where ¯ ¸ 0 is the decay coefficient, and t represents the time. We will prove four statements: (I)8T >T 0 ,9¯ ¤ s.t. @f(¯ ¤ ;t) @t j t=T =0 (II)80<t<T , @f(¯ ¤ ;t) @t >0 (III)8¯ >¯ ¤ , @f(¯;t) @t j t=T <0 (IV)80·¯ <¯ ¤ , @f(¯;t) @t j 0<t·T >0 95 Letg(¯;t)= @f(¯;t) @t . Using (A.2),g(¯;t) can be written as: g(¯;t)= @f(¯;t) @t =C p (t)¡¯f(¯;t) (A.3) or g(¯;t)= @(C p (t)¡¯f(¯;t)) @t =C p (0)e ¡¯t + Z t 0 C 0 p (¿)e ¡¯(t¡¿) d¿ =C p (0)e ¡¯t +e ¡¯t Z t 0 e ¯¿ C 0 p (¿)d¿ (A.4) (I)8T >T 0 ,9¯ ¤ s.t. @f(¯ ¤ ;t) @t j t=T =0 Proof: Fixingt = T , g(¯;T) is a continuous function of¯. From (A.3), when¯ = 0, g(0;T) = C p (T) > 0. If we can show thatg(¯;t)j t=T < 0 for some¯ > 0, then there must exist a¯ ¤ such thatg(¯ ¤ ;T)=0. LetT + 0 be a time point s.t. T 0 <T + 0 <T . Using equation (A.4): g(¯;T)=e ¡¯T · C p (0)+ Z T 0 e ¯¿ C 0 p (¿)d¿ ¸ =e ¡¯T " C p (0)+ Z T 0 0 e ¯¿ C 0 p (¿)d¿ + Z T + 0 T 0 e ¯¿ C 0 p (¿)d¿ + Z T T + 0 e ¯¿ C 0 p (¿)d¿ # ·e ¡¯T " C p (0)+ Z T 0 0 e ¯¿ C 0 p (¿)d¿ + Z T T + 0 e ¯¿ C 0 p (¿)d¿ # (A.5) There exists a positive number 0 < M < 1 s.t. C 0 p (t) · ¡M; 8t 2 [T + 0 ;T]. There exist another positive number0<N <1 s.t. C 0 p (t)·N;8t2[0;T 0 ]. Follow (A.5): 96 g(¯;T)·e ¡¯T " C p (0)+ Z T 0 0 e ¯¿ Nd¿¡ Z T T + 0 e ¯¿ Md¿ # =e ¡¯T " C p (0)+N e ¯T 0 ¡1 ¯ ¡M e ¯T ¡e ¯T + 0 ¯ # < e ¡¯T ¯ h C p (0)¯e ¯T + 0 +Ne ¯T + 0 +Me ¯T + 0 ¡Me ¯T i = Me ¡¯(T¡T + 0 ) ¯ · C p (0)¯+N +M M ¡e ¯(T¡T + 0 ) ¸ (A.6) Notice that T ¡T + 0 > 0, the term Cp(0)¯+N+M M increases linearly with ¯, whereas the term e ¯(T¡T + 0 ) increases exponentially with ¯. So, there exist a ¯ g<0 , ¯ g<0 > 0, such that g(¯ g<0 ;T) < 0. Since g(0;T) > 0, g(¯ g<0 ;T) < 0 and g(¯;T) is a continuous function with respect to¯, there must exist a¯ ¤ 2[0;¯ g<0 ] satisfyingg(¯ ¤ ;T)=0. (II)80<t<T , @f(¯ ¤ ;t) @t >0 Proof: When fix ¯ ¤ , g(¯ ¤ ;t) is a function of t. By using (A.4), it’s easy to show that g(¯ ¤ ;t)>0 fort2(0;T 0 ]. Considering the case fort2(T 0 ;T), from (A.3), taking the partial derivative ofg(¯;t) with respect tot yields: @g(¯;t) @t =C 0 p (t)¡¯ @f(¯;t) @t =C 0 p (t)¡¯g(¯;t) (A.7) We have shown that g(¯ ¤ ;t)j t=T = 0. If g(¯ ¤ ;t) is not always positive in (T 0 ;T), suppose there exist a t 0 2 (T 0 ;T), s.t. g(¯ ¤ ;t 0 ) · 0. Then there must also exist a t 00 2 (t 0 ;T) s.t. @g(¯ ¤ ;t) @t j t=t 00 ¸ 0,C 0 p (t 00 ) < 0 andg(¯ ¤ ;t 00 ) = 0, which contradicts with (A.7). Therefore g(¯ ¤ ;t) > 0 for t 2 (0;T 0 ). Combine the two cases above, we have 80<t<T , @f(¯ ¤ ;t) @t >0 97 (III)8¯ >¯ ¤ , @f(¯;t) @t j t=T <0 Proof: We first show that @ 2 f(¯;t) @¯@t j ¯=¯ ¤ <0; t2(0;T]. Continue from (A.4), @ 2 f(¯;t) @¯@t =¡tC p (0)e ¡¯t ¡ Z t 0 (t¡¿)C 0 p (¿)e ¡¯(t¡¿) d¿ (A.8) This equation implies that @ 2 f(¯;t) @¯@t < 0, 8¯ ¸ 0 and 8t 2 (0;T 0 ]. To show @ 2 f(¯;t) @¯@t j ¯=¯ ¤ < 0 for t 2 (T 0 ;T), the proof by contradiction is used. Proceed from (A.8): @ 3 f(¯;t) @t@¯@t =¡C p (0)e ¡¯t +¯tC p (0)e ¡¯t ¡ Z t 0 C 0 p (¿)e ¡¯(t¡¿) d¿ +¯ Z t 0 (t¡¿)C 0 p (¿)e ¡¯(t¡¿) d¿ =¡C p (0)e ¡¯t ¡ Z t 0 C 0 p (¿)e ¡¯(t¡¿) d¿ +¯ · tC p (0)e ¡¯t + Z t 0 (t¡¿)C 0 p (¿)e ¡¯(t¡¿) d¿ ¸ =¡ @f(¯;t) @t ¡¯ @ 2 f(¯;t) @¯@t (A.9) Suppose there exists a t 0 2 (T 0 ;T), such that @ 2 f(¯;t) @¯@t j ¯=¯ ¤ ;t=t 0 = 0. Without loss of generality, assume thatt 0 , if it exists, is the minimum one among all possible solutions. So @ 2 f(¯;t) @¯@t <0 for¯ =¯ ¤ ; t2(T 0 ;t 0 ). Using equation (A.9) @ 3 f(¯;t) @t@¯@t j ¯=¯ ¤ ;t=t 0 =¡ @f(¯;t) @t j ¯=¯ ¤ ;t=t 0 <0 (A.10) where the ’<’ sign comes from the conclusion of statement (II). Remind that @ 2 f(¯;t) @¯@t <0 fort2 (0;T 0 ]. Ast increases fromT 0 toT , due to our assumption, t 0 is the first point that the function @ 2 f(¯;t) @¯@t equals zero. This implies that @ 3 f(¯;t) @¯@ 2 t ¸ 0j ¯=¯ ¤ ;t=t 0 , which 98 contradicts with (A.10). The original assumption must be wrong, and there does not exist at 0 s.t. @ 2 f(¯;t) @¯@t j ¯=¯ ¤ ;t=t 0 =0 fort2(0;T). Therefore, @ 2 f(¯;t) @¯@t j ¯=¯ ¤ <0; t2(0;T) (A.11) To show @ 2 f(¯;t) @¯@t j ¯=¯ ¤ ; t=T <0, we use equation (A.8), @ 2 f(¯;t) @¯@t j ¯=¯ ¤ ;t=T =¡TC p (0)e ¡¯ ¤ T ¡ Z T 0 (T ¡¿)C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ =¡T · C p (0)e ¡¯ ¤ T + Z T 0 C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ ¸ + Z T 0 ¿C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ =¡Tg(¯ ¤ ;T)+ Z T 0 ¿C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ (A.12) Becauseg(¯ ¤ ;T)=0, we have: @ 2 f(¯;t) @¯@t j ¯=¯ ¤ ;t=T = Z T 0 ¿C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ = Z T 0 0 ¿C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ + Z T T 0 ¿C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ < Z T 0 0 T 0 C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ + Z T T 0 T 0 C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ =T 0 Z T 0 C 0 p (¿)e ¡¯ ¤ (T¡¿) d¿ =T 0 £ g(¯ ¤ ;T)¡C p (0)e ¡¯ ¤ T ¤ =¡T 0 C p (0)e ¡¯ ¤ T ·0 (A.13) I.e. @ 2 f(¯;t) @t@¯ j ¯=¯ ¤ ; t=T <0. Combine this with (A.11), we get: @ 2 f(¯;t) @t@¯ j ¯=¯ ¤ <0; t2(0;T] (A.14) 99 Next, we will show statement (III) by using (A.14). Suppose that (III) does not hold, then @f(¯;t) @t j t=T crosses zero at least one time for ¯ > ¯ ¤ . Let ¯ ¤+ be the small- est one among these points. Notice that, same as ¯ ¤ , ¯ ¤+ satisfies (II), (A.14) and @ 2 f(¯;t) @t@¯ j ¯=¯ ¤+ ; t=T < 0, which contradicts with the assumption that ¯ ¤+ is the smallest value that makes @f(¯;t) @t j t=T cross zero. Therefore, 8¯ >¯ ¤ ; @f(¯;t) @t j t=T <0 (A.15) (IV) 80·¯ <¯ ¤ , @f(¯;t) @t j 0<t·T >0 Proof: Using the same technique in (III), we can show that80·¯ <¯ ¤ , @f(¯;t) @t j t=T > 0. Thus, we only need to show: 80·¯ <¯ ¤ ; @f(¯;t) @t j 0<t<T >0 (A.16) Use the proof by contradiction. Suppose (A.16) does not hold, since @f(¯;t) @t is a contin- uous function, there must exist a ¯ ¤¡ ;¯ ¤¡ 2 [0;¯ ¤ ) and a T ¡ ;T ¡ 2 (0;T) such that @f(¯;t) @t j ¯=¯ ¤¡ ; t=T ¡ = 0. ¯ ¤¡ and T ¡ are a pair of ¯ ¤ and T values that satisfy (I), (II) and (III). Rewriting (III) for¯ ¤¡ andT ¡ , we have: 8¯ >¯ ¤¡ ; @f(¯;t) @t j t=T ¡ <0 (A.17) Since ¯ ¤ > ¯ ¤¡ and T > T ¡ , the equation above implies that @f(¯;t) @t j ¯=¯ ¤ ; t=T ¡ < 0, which contradicts with (II). Therefore, equation (A.16) must hold. Finally, we have: 80·¯ <¯ ¤ ; @f(¯;t) @t j 0<t·T >0 (A.18) 100

Abstract (if available)

Abstract Positron emission tomography (PET) with F-18 fluorodeoxyglucose (FDG) is improving the diagnosis, staging and treatment monitoring of a variety of cancers in comparison to alternative noninvasive imaging modalities. However, detection of small lesions is limited by image resolution and low signal to noise ratio. While computer-aided detection algorithms based on static PET may assist visual detection of lesions, the temporal domain information loss in the static image could impair the detection of lesions. In addition, spatial features such as the shape and contrast of a lesion in static PET images are highly variable and difficult to characterize. Since the dynamics of tracer uptake (tracer radioactivity) are different for normal and malignant tissue, the temporal domain provides potentially important additional features for use in lesion detection.

Linked assets

University of Southern California Dissertations and Theses

Asset Metadata

Creator Li, Zheng (author)

Core Title Statistical lesion detection in dynamic positron emission tomography

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 03/23/2008

Defense Date 12/13/2007

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

Tag computer-aided detection (CAD),dynamic PET,lesion detection,OAI-PMH Harvest,statistical

Language English

Advisor Leahy, Richard M. (committee chair), Mikulevicius, Remigijus (committee member), Narayanan, Shrikanth S. (committee member), Ortega, Antonio (committee member), Yu, Xiaoli (committee member)

Creator Email zhengli@usc.edu

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

Unique identifier UC193959

Identifier etd-Li-20080323 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-40619 (legacy record id),usctheses-m1053 (legacy record id)

Legacy Identifier etd-Li-20080323.pdf

Dmrecord 40619

Document Type Dissertation

Rights Li, Zheng

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