Sample Pages

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/262603426 Sample Pages DATASET · MAY 2014 READS 19 2 AUTHORS: Shuo Li The University of Western O… 116 PUBLICATIONS 507 CITATIONS SEE PROFILE Joao Manuel R. S. Tavares University of Porto 593 PUBLICATIONS 1,765 CITATIONS SEE PROFILE Available from: Shuo Li Retrieved on: 08 February 2016 Shape Analysis in Molecular Imaging Fei Gao and Pengcheng Shi Abstract Molecular imaging is a new research discipline enabling the visualization, characterization and quantification of biologic processes taking place at the cellular and subcellular levels within intact living subjects. Applications of molecular imaging techniques will benefit various clinical practices including classification and tracking of chemotherapy and treatment planning of radiotherapy, as well as drug discovery and development. Molecular imaging typically includes two or three dimensional imaging with quantification over time, and is often applied on molecular imaging modalities, such as Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT) etc. Image series acquired with spatiotemporal distribution of molecular biomarkers must be carefully analyzed to estimate the underlying physiology-related metabolic parameters. Shape analysis is one of the most powerful tools to analyze the geometrical properties from similar shapes or different groups, and can be applied to estimate both the concentration of biomarkers and interaction between biomarkers and tissue/organs. However, some limitations from molecular imaging modalities and clinical practices still hinder the quantitative accuracy of shape analysis, e.g. the low spatial and temporal resolution in PET scan, the inaccuracy of blood samplings from patients, the low Signal-toNoise (SNR) ratio of measurement data in dynamic PET/CT scan. In this chapter, firstly, we will introduce the definition of molecular imaging, the clinical advantages and limitations of various molecular imaging modalities, secondly, we will review the challenges in data analysis based on the data processing procedure, and explain how data corrections affect the accuracy of static and dynamic PET imaging, thirdly, the general frameworks of image processing in PET and SPECT are reviewed with focus on image reconstruction, at last, we will show some recent advancements and give examples of clinical applications. F. Gao (B) and P. Shi Rochester Institute of Technology, 102 Lomb Memorial Drive, Rochester, NY 14623-5608, USA e-mail: gaofei.rit@gmail.com P. Shi e-mail: pengcheng.shi@rit.edu S. Li and J. M. R. S. Tavares (eds.), Shape Analysis in Medical Image Analysis, Lecture Notes in Computational Vision and Biomechanics 14, DOI: 10.1007/978-3-319-03813-1_2, © Springer International Publishing Switzerland 2014 51 52 F. Gao and P. Shi 1 Introduction to Molecular Imaging and Molecular Imaging Modalities Molecular imaging provides the images of molecular and cellular level activities inside the body. Molecular imaging enables doctors to measure the biological processes quantitatively and reflects the functionality of organs and tissues inside patients. According to the definition from the Society of Nuclear Medicine and Molecular Imaging (SNMMI), molecular imaging is the visualization, characterization, and measurement of biological processes at the molecular and cellular levels in humans and other living systems [58]. Molecular imaging usually involves two- or three- dimensional images and conducts quantitative analysis over the time. Molecular imaging is a noninvasive procedure and can be used to study and diagnosis of cancer, brain diseases and disorders, cardiology, and various disorders in different organs and tissues. Molecular imaging modalities include Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT), Optical imaging, Magnetic Resonance Imaging (MRI), Computed Tomography (CT), and Ultrasound (US). Of all modalities, PET and SPECT are the full functional imaging modalities while others are only with limited abilities. Based on these imaging modalities, hybrid PET/CT [90], PET/MRI [74], PET/SPECT/CT [56] further enrich the ability of molecular imaging. Different imaging modalities require different considerations of data processing methods. The following sections will introduce the major molecular imaging modalities and their data processing methods. 1.1 Positron Emission Tomography PET as a biomedical research technique and clinical diagnostic procedure is one of the most important applications in nuclear medical imaging devices. In the past three decades, there have been significant advancements in PET scanners and image processing methods [1, 92, 94]. PET scan is a unique type of imaging test that helps doctors see how the organs and tissues inside your body are actually functioning. PET scan reveals the cellular level metabolic changes occurring in an organ or tissue. This is important and unique because disease processes often begin with functional changes at the cellular level. Currently, PET scans are most commonly used to detect cancer, heart problems, brain disorders and other central nervous system disorders. PET scan can be used to track the spread of disease inside body and patient response to drugs and therapies, which help to determine the more effective treatment plans for individual patient. PET scans can also be used to follow-up and manage ongoing cares. Quantitative dynamic PET imaging also offers good promise for personalized drug treatment by accurate pharmacokinetic analysis and will enable medicine to be tailored to each person’s needs, and improve the safety, quality and effectiveness of healthcare for every patient. Shape Analysis in Molecular Imaging (a) (c) 53 (b) (d) Fig. 1 PET images in different views: a Coronal view; b Sagittal view; c Horizontal view; d 3D view PET scans rely on the injected radiotracers which circulate inside the body. PET scanner detects the pairs of gamma photons emitted from the radiotracer, which are known as positron annihilation, and creates a projection data of radiotracer distribution. The reconstructed images from projection data are used for diagnosis. PET images of different views from a volunteer are shown in Fig. 1. One currently used PET radiotracer for daily clinical routines is 18 F-FDG (Fluoro-2-deoxy-d-glucose), which is a compound consist of glucose and radioactive fluorine-18. When disease occurs, the activities of cells begin to change, for example, the cancer cells need more glucose and more active than normal cells, so there will be more radiolabeled 18 FFDG accumulated in the cancer cells. With the PET images, it appears higher intense than surrounding tissues, which is called as a ‘hot spot’ indicating a high level of activity or metabolism is occurring there. Correspondingly, a ‘cold spot’ refers to the area of low metabolic activities indicated by lower intense in the PET images. With these PET images, doctors will be able to evaluate the working situation of organs 54 (a) F. Gao and P. Shi (b) (c) Fig. 2 PET/CT images: a CT image; b PET image; c Fused PET/CT image and tissues and determine the abnormals by analyzing the hot spots and cold spots. The ability of detecting cellular level change that occurs early in the disease makes PET imaging superior to the CT and MR images which show the structural change after accumulations in cellular level. Hybrid PET/CT or PET/MRI is the combination of PET and CT or MRI, the combination of multiple imaging modalities allows both anatomical and functional images in one image set. One example of PET/CT images is shown in Fig. 2. Different radiotracers will reveal different diseases. Besides 18 F-FDG mentioned before, which is widely used for cancer diagnosis, cardiology, neurology, there are many other radiotracers used in research and clinical applications. 18 F-FLT (3′ fluoro-3′ -deoxy-l-thymidine) is developed as a PET tracer to image tumor cell proliferation [12], 11 C-Acetate is developed to localize prostate cancer [62], 13 N-ammonia is developed to quantify the myocardial blood flow [48], 11 C-dihydrotetrabenazine (DTBZ) is developed for brain imaging, which can be used for differentiating Alzheimer’s disease from dementia and Parkinson’s disease [47], 11 C-WIN35,428 is a cocaine analogue and sensitive to the dopamine transporters [93]. Researchers are working on labeling different drugs with radioactive 11 F, 11 C, 13 N etc. for PET scan, and pharmaceutical companies are especially interested in applying quantitative PET analysis on radio-labeled new drugs, this has the potential to shorten the Phase I to Phase II studies from more than 5 years to be 2 years. Clinical PET studies include activity image reconstruction of radiotracer concentrations and parametric image reconstruction from dynamic PET scans. The quantitative accuracy of reconstructed images depends on the whole procedure: data corrections of measurement data, statistical modeling of acquisition process, and proper reconstruction methods. The data correction is the very first step, and directly determines the accuracy of following steps. The data correction consists of many specific corrections, including random correction, scatter correction, deadtime cor- Shape Analysis in Molecular Imaging 55 Fig. 3 One patient SPECT scan rection, attenuation correction etc. Among all the corrections, scatter correction is the most complicated and still a very active research topic. The percentage of scatter coincidence events in measurement data is around 30 % for PET scanners from major manufactures. And their distributions are more complicated and will vary with the dosage of injected radiotracer and the status of detector system, furthermore, the scatter coincidence events can either be corrected or modeled in the system probability matrix. Additionally, the development of new full 3D PET scanner with new scintillators and Photomultiplier Tubes (PMT) [3, 6, 28, 40, 52, 87], requires corresponding new correction methods. 1.2 Single Photon Emission Computed Tomography SPECT scan uses a gamma camera that rotates around the patient to detect the radiotracer inside body. SPECT will also produce a set of 3D images but generally have a lower resolution. The radiotracers commonly used for SPECT scan include 99m T c [59], 188 Re [39], 68 Ga [104], 82 Rb [23] etc. Electrocardiography (ECG)-Gated 82 Rb can also be used for myocardial perfusion PET [5]. One example image from SPECT scan is shown in Fig. 3. Hybrid SPECT/CT is also designed to provide more accurate anatomical and functional information [88]. SPECT scan differs from PET scan in that the tracer stays in your blood stream rather than being absorbed by surrounding tissues, therefore, SPECT scan can show how blood flows to the heart and brain are effective or not. SPECT scan is cheaper and more readily available than higher 56 F. Gao and P. Shi Fig. 4 One patient brain MRI scan resolution PET scan. Tests have shown that SPECT scan might be more sensitive to brain injury than either MRI or CT scanning because it can detect reduced blood flow to injured sites. SPECT scan is also useful for presurgical evaluation of medically uncontrolled seizures and diagnosing stress fractures in the spine (spondylolysis), blood deprived (ischemic) areas of brain following a stroke, and tumors [4, 9]. 1.3 Other Molecular Imaging Modalities Other molecular imaging modalities include MRI, CT, Optical imaging, Ultrasound, which have only limited molecular imaging abilities. 1.3.1 MRI MRI scanners produce strong magnetic fields where body tissue that contains hydrogen atoms is made to emit a radio signal which is detected by the scanner. MRI scanners produce detailed 3D images of the inside of the body. MRI scan is best for brain imaging, breast, heart and blood vessel, organs and soft tissues [13]. Contrast materials are often used to enhance MRI images. MRI has the highest spatial resolution and possibility to extract both physiological and anatomical information. However, MRI generally has low sensitivity and requires longer scan and data processing time. Samples of MRI images are shown in Fig. 4. Functional magnetic resonance imaging or functional MRI (fMRI) is a MRI procedure that measures brain activity by detecting associated changes in blood flow. The primary form of fMRI uses the blood-oxygen-level-dependent (BOLD) contrast, and measures the change in magnetization between oxygen-rich and oxygen-poor blood, which is mostly used in brain mapping research [37]. Shape Analysis in Molecular Imaging 57 Another MRI scan used in research is Dynamic Contrast Enhanced- MRI (DCEMRI), which is a quantitative method that allows for tumor vascular analysis, including blood volume, perfusion, vascular leakage space. DCE-MRI uses gadolinumbased contrast agents and standard MRI scanners to provide quantitative results [34]. 1.3.2 Optical Imaging Optical imaging has the highest sensitivity, but relatively low spatial resolution. Optical imaging only provides images in limited Field of View (FOV). Optical imaging reduces patient radiation exposed significantly by using non-ionizing radiation. Optical imaging can be used to differentiate between native soft tissues and tissues labeled with either endogenous or exogenous contrast media, using their different photon absorption or scattering profiles at different wavelengths. Optical imaging is also easy to be combined with other imaging modalities. In optical imaging, Diffusive Optical Imaging (DOI) is also known as Diffuse Optical Tomography (DOT) or Optical Diffusion Tomography (ODT). DOI is used to study the functions of brains, and can provide neural activities with its time courses [51]. Optical Coherence Tomography (OCT) produces 3D images from optical scattering media and penetrates biological tissues by using long wavelength light. OCT can provide higher SNR, faster signal acquisition [78]. 1.3.3 Ultrasound Ultrasound scanners emit high frequency sound and measure the reflected sound from the patients, which varies with the organs and tissues. Ultrasound imaging can be used to detect heart problems, examine liver, kidneys, abdomens, and guide a surgeon. Ultrasound can provide a real time imaging, but with limited spatial resolution. Most scans need targeted micro-bubbles to enhance the images [21]. 2 The Challenges in Molecular Imaging Data Analysis 2.1 Data Processing for PET Imaging The previous sections explain the advantages and limitations of molecular imaging modalities. Of all molecular imaging modalities, PET and SPECT are the full functional ones. Since PET and SPECT all belong to Emission Computed Tomography (ECT), many data correction and data analysis methods can be shared. In this section, we will focus on PET to demonstrate the challenges in data analysis. As shown in Fig. 5, the measurement data of PET is first stored in listmode. The listmode data 58 F. Gao and P. Shi Fig. 5 The whole procedure of PET data processing can be corrected and reconstructed directly [66, 72], or rebinned by Fourier Rebinning to sinogram, which is corrected and reconstructed to be images [18]. The data correction process generally includes random correction, normalization, deadtime correction, attenuation correction and scatter correction. The corrected data will be reconstructed by either analytical or statistical algorithms to generate both activity images which are the concentrations of radiotracers in body and parametric images which are the dynamic changes of radiotracers represented by kinetic parameters. The challenges in PET data analysis come from the change of statistical properties of measurement data after various data corrections. The quality of results from all image reconstruction algorithms depends on the accuracy of statistical models in each data correction and image reconstruction. However, due to the complexity of PET scan, it is nearly impossible to propose a perfect model. From next subsection, we will explain how the data corrections are applied and how they affect the data analysis. 2.1.1 Data Corrections 1. Random Correction. Random coincidence events are two gamma photons from different positron annihilations are recorded as a Line of Response (LOR). The rate of random coincidences on a particular LOR is defined as Ri j = 2τri r j , where Ri j is the random coincident rate on the LOR connecting detector i and j, τ is the coincidence timing window, ri and r j are the singles rate on detector Shape Analysis in Molecular Imaging 2. 3. 4. 5. 59 i and j. The random coincidence rate depends on the detector circuit response, the most common method for estimating the random coincidence events is the delayed timing window method. Normalization. All the image processing methods assume the detector response are uniform and identical, however, in real clinical data acquisition, the performance of each detector will be different, the PMT gains are not exactly the same, and angle of incident of each gamma photon will also affect the detector sensitivity. The process of correcting these effects is normalization. One common clinical solution is to perform a scan using a uniform cylindrical phantom and a line source. In this scan, every possible LOR is illuminated by the same coincidence source. Deadtime Correction. The detector deadtime will cause the loss of coincidence events. Especially in the high count rate cases, the detector response will be significantly delayed due to the pile-up of incident gamma photons. A common solution is to model both paralyzable and non-paralyzable components from measurements of sources of different radioactivities. Attenuation Correction. Some gamma photons will be absorbed when interacting with the patients. The widely-adopted method for PET only scanner calculates the attenuation correction factors by using a rotating rod source with the object in the FOV. For hybrid PET/CT and PET/MRI scanners, the attenuation factors can be obtained by assigning predefined attenuation coefficients to the anatomical information generated by CT and MRI, which generally have a better spatial resolution. Scatter Correction. The coincidence events include those that the annihilation photons have an interaction of Compton scattering and change their directional and energy information before they arrive at the detector system. They are called scatter coincidences. Scatter coincidences will decrease the contrast, resolution and SNR of reconstructed images and need to be corrected properly. Scatter correction is the most complicated correction, we will use the scatter correction as an example to show how corrections affect the quantitative accuracy in PET scan. 2.1.2 Scatter Correction Scatter Fraction. The scatter fraction is an important parameter included in both National Electrical Manufacturers Association (NEMA) and International Electrotechnical Commission (IEC) standards, and measured in every performance evaluation report of PET scanners. The scatter fractions of PET scanners from 3 major manufactures are listed in Table 1 [28, 40, 87]. Multiple Scatter Coincidence Events. The gamma photons can be scattered multiple times before recorded. Monte Carlo simulation is the best tool to study the Compton scattering during PET acquisition. We have done similar studies with our new PET scanner, HAMAMATSU SHR74000, we retrieve all the data, and analyze the composition of simulated measurement data according to the number of 60 F. Gao and P. Shi Table 1 Scatter fractions of PET scanners from three main manufactures in 3D mode Manufacture and model Scatter fraction (%) Siemens TruePoint PET/CT GE discovery PET/CT Philips Gemini PET/CT 32 33.9 27 Table 2 The percentage of coincidence events assorted by the number of scattering of each photon Photon2 \Photon1 0× (%) 1× (%) 2× (%) 3× (%) 4× (%) 0× 1× 2× 3× 4× 55.36 16.56 1.92 0.11 0.005 16.58 5.81 0.7 0.04 0.001 1.92 0.68 0.08 0.005 0 0.11 0.04 0 0 0 0 0 0 0 0 The data is from the simulation of a line source fixed in the transaxial center of standard phantom Table 3 The percentage of coincidence events assorted by the number of scattering of each photon Photon2 \Photon1 0× (%) 1× (%) 2× (%) 3× (%) 4× (%) 0× 1× 2× 3× 4× 37.66 18.82 3.75 0.38 0.02 18.87 10.98 2.18 0.23 0.02 3.7 2.18 0.44 0.05 0 0.39 0.23 0.05 0.01 0 0.03 0.01 0 0 0 The data is from the simulation of a line source fixed in the transaxial center of over-sized phantom scattering of each photon in the coincidence photon pairs. The first data set is analyzed with the standard NEMA phantom, a 70 cm long, 20 cm diameter cylinder, the composition of scatter coincidences is summarized in Table 2. Furthermore, we have a lot of over-sized and over-weighted patients, especially in the western countries, and these patients usually take higher risk of many diseases than other patients. For these patients, their size makes more photons scattered inside their body, we correspondingly analyze the composition of scatter coincidence photons with a 70 cm long, 35 cm diameter cylinder, the results are summarized in Table 3. The number of total scatter coincidence events increase by about 20 %. If the tumor is near the body surface, the scatter fractions shown in Table 4 turn to be similar as in Table 2. Scatter Coincidence from Activity Concentrations outside of FOV. A whole body PET scan generally needs 4–5 bed positions to complete, then gamma photons from outside of FOV will also have the possibility to be misrecorded by the detector system. Sossi et al. show the scatter coincidence events need to be accounted for when the amount of radioactivity outside of FOV is comparable to the radioactivity inside the FOV in [81]. In the case of line sources, the scatter fraction increases by about 15 % when the source is extended 4 cm outside of the FOV and 20 % when it is extended 8 cm outside of the FOV. Spinks et al. also show for a full 3D PET, Shape Analysis in Molecular Imaging 61 Table 4 The percentage of coincidence events assorted by the number of scattering of each photon Photon2 \Photon1 0× (%) 1× (%) 2× (%) 3× (%) 4× (%) 0× 1× 2× 3× 4× 59.23 16.29 2.94 0.38 0.04 16.28 1.18 0.1 0.01 0 2.94 0.1 0 0 0 0.38 0.01 0 0 0 0.03 0 0 0 0 The data is from the simulation of a line source fixed in the transaxial 169 mm offset of over-sized phantom the scatter fraction increased from 40 to 45 % with the scatter photons from outside of FOV, after introducing a fitted brain shielding, the scatter fraction decreases to be 41 % [82]. A recent paper by Ibaraki et al. show the use of the neck-shield to suppress scatters from outside FOV can improve the SNR by 8 and 19 % for H2 15 O and 15 O2 , respectively [38]. Scatter Correction Methods. Bailey et al. [2] and Bentourkia et al. [7] present a Convolution-Subtraction (CS) scatter correction technique for 3D PET data. The scatter distribution is estimated by iteratively convolving the photo-peak projections with a mono-exponential kernel. The method is based on measuring the scatter fraction and the scatter function at different positions in a cylinder. The method performs well on 2D measurement data and also accounts for the 3D acquisition geometry and nature of scatter by performing the scatter estimation on 2D projections. The method is easy to set up and still applies to a lot of animal studies, where the scatter correction are usually considered not necessary. Lubberink et al. propose a non-stationary CS scatter correction with a dual-exponential scatter kernel for scatter correction in both emission and transmission data of studies of conscious monkeys using Hamamatsu SHR7700 PET scanner [54] . Kitamura et al. implement a hybrid scatter correction method, which estimates scatter components with a dual energy acquisition using a CS to estimate the true coincidence events in the upper energy window for their four layer Depth of Interaction (DOI)—PET scanner [45]. NaidooVariawa et al. suggest that scatter correction methods based on spatially invariant scatter functions, such as CS, may be suitable for non-human primate brain imaging in [65]. Single Scatter Simulation (SSS) is one of the most important methods for PET scatter correction, and commercial PET scanners from Siemens are using SSS derived scatter correction methods. Ollinger [67] and Watson [98] first introduced the concept of model based scatter coincidence estimation. The single scatter approximation is defined with the well accepted formulation based on Klein-Nishina equation as S AB =   v σ AS σ B S 4π R 2AS R 2B S  µ dσc [I A + I B ], σc dΩ (1) 62 F. Gao and P. Shi Fig. 6 The procedure of a SSS based scatter correction for clinical PET scanner where I A = ε AS (E)ε B S (E ′ )e− ′ I B = ε AS (E )ε B S (E)e − A µ(E,x)ds − A µ(E ′ ,x)ds − S S e e B µ(E ′ ,x)ds B µ(E,x)ds S S × ×   A λds S B λds S v is the sampled scatter volume, S is the scatter point, A and B are sampled detectors, so S AB is the coincidence scatter rate in detectors A and B, σ AS and σ B S are the detector cross sections presented to the rays AS and B S, R AS and R B S are the distances from S to A and B. σc and dσc are the total and differential Compton scattering cross sections, Ω is the solid angle, E is the energy of the non-scattered photon, E ′ is the energy of the scattered photon, ε AS (E) and ε B S (E) are the approximated detector efficiencies for gamma rays which are incident along AS and B S, µ(E, x) is the linear attenuation coefficient at the energy E and position x, and λ is the estimated activity distribution from reconstructed images. The calculation is repeated through all the sampled scatter points from activity distribution and sampled detector blocks. The final result is the distribution of single scatter coincidence events, and will be scaled and then subtracted from the measurement data to perform scatter correction. Figure 6 demonstrates the procedure of SSS based scatter correction for clinical PET scanner. 2.1.3 Dynamic PET Imaging Dynamic PET imaging is a combination of short interval PET scans and reflects the dynamic metabolism of injected radiotracers. For example, a dynamic acquisition can Shape Analysis in Molecular Imaging (c) 1.2 ROI1 1.0 ROI2 ROI3 Concentration ( µCi/ml) (b) Concentration (µCi/ml) (a) 63 0.8 0.6 0.4 0.2 20 40 60 80 Time (min) 100 120 2.0 True False 1.5 1.0 0.5 20 40 60 80 100 120 Time (min) Fig. 7 a Standard zubal phantom; b TAC curves of three ROIs indicated in (a); c sample TAC curves of true lesion and false lesion consist of 85 frames in all: 15 × 0.2 min, 20 × 0.5 min, 40× 1 min, and 10 × 3 min. The series of acquisitions can be used to estimate the kinetic parameters which represent the metabolism of radiotracers in vivo. The difficulties in estimating kinetic parameters arise from the low count rates in the first several time frames and the low SNR. To obtain the kinetic parameters, a typical approach is first to reconstruct the activity distributions from the dynamic PET data, and then to fit the calculated time activity curve (TAC) to a predefined kinetic model. Samples of TAC curves in different organs are shown in Fig. 7. The accuracy of this kind of approaches relies on the reconstructed activity distributions. The complicated statistical noise properties, especially in the low-count dynamic PET imaging, and the uncertainties introduced by various PET data corrections will affect the activity reconstruction and lead to a suboptimal estimation of kinetic parameters [29]. There are also many efforts that try to estimate the kinetic parameters from PET projection data directly and achieve better bias and variance including both linear and nonlinear models [83, 91, 103]. One problem is that the optimization algorithms are very complicated. Kamasak et al. apply the coordinate descent algorithm for optimization but it is still limited to specific kinetic models [61]. Wang et al. apply a generalized algorithm for reconstruction of parametric images [96], however, it still lacks of estimation and analysis of individual kinetic parameter. Exactly, every kinetic parameter has its own physical meaning like radiotracer transport rate, phosphorylation rate and dephosphorylation rate (it is very sensitive to the system [35]) in FDG study, which will be critical to clinical research, drug discovery and drug development [15, 92]. 64 F. Gao and P. Shi 3 General Framework of Shape Analysis in Molecular Imaging 3.1 Image Segmentation In computer vision, image segmentation is the process of partitioning an image into multiple different segments (group of pixels). Especially in molecular imaging, the image segmentation is used to simplify the representation of an image and extract Region of Interest (ROI) that is more meaningful and followed by image analysis. Image segmentation is also important to find the boundaries of different regions and organs by applying different labels. Image segmentation can also be applied to 3D image stacks to help 3D image reconstruction [17, 106]. The aforementioned limitations of PET and SPECT also bring new challenges to image segmentation. Here we will introduce several basic image segmentation methods. 1. The simplest method of image segmentation is the thresholding method. Because of the simplicity and fast implementation, the thresholding method is still used in some clinical routines to identify ROI with high contrast, e.g. the lesions in lung [19, 20, 42]. 2. Cluster based method is a multivariate data analysis method that uses predefined criteria to partition a large number of objects into a smaller number of clusters, in which the objects are similar to each other. Cluster based method has been applied to fMRI imaging and then dynamic PET imaging, with limited spatial resolution and SNR [101]. 3. Gradient based method is to find the boundary of an object of interest with the gradient intensity observed in the image. This kind of method is fast and easy to apply, but generally works together with other method to achieve better results in molecular imaging [26]. 4. Level set based method derives from snake method (active contour model), which delineates an object outline from a noisy image by attempting to minimize an energy associated to the contour as a sum of internal and external energies [57, 68]. Many methods are derived from basic level set method including deformable level set models, which have the ability to automatically handle topology [33], 3D level set methods to compute 3D contours [105] etc. 5. Kinetic model guided segmentation assumes different ROIs have different tracer kinetic properties to seprate different functional regions [11]. This method can also be used to estimate the input functions for quantitative dynamic PET data analysis [71]. Shape Analysis in Molecular Imaging 65 3.2 Image Registration Image registration is the process to transform different sets of data into one coordinate system. Image registration is widely used in molecular imaging, e.g. patient radiotherapy follow-up by transforming PET images from a series of studies, diagnosis by images from multiple imaging modalities [16, 32, 36, 75]. With the development of hybrid PET/CT, PET/MRI, the image registration with multiple images in one study is made easier because the motions of patients are minimized by the simultaneous data acquisition. However, images from multiple studies still need good image registrations. Mainstream image registration methods for molecular imaging include 1. Intense-base image registration. Since PET and SPECT imaging reflects the concentrations of radiotracer, intense based methods compare intensity patterns in multiple images and register the reference image and target image by defining correlation metrics [46]. 2. Feature-based image registration. In a series of studies/images, common features can be extracted from the anatomical information of organs and tissues which do not change a lot and can be used as references [73]. This method can also be used for multiple imaging modalities [55]. 3. Multiple modality image registration. The hybrid PET/CT and PET/MRI make the image registration focus on the deformations from patients’ respirations and motions [27, 60]. The image registration can be improved by different patient preparation and pre-positioning [8], respiratory gating [10], various tracking devices etc. [76]. 3.3 Image Fusion Image fusion is the combination of relevant information from two or more images into one single image. The fused image will provide more information than any single input image. Accurate image fusion from combined PET, CT, MRI scans can significantly improve the diagnosis and provide better understandings of diseases. Image fusion generally works together and shares similar technologies with image registration [69, 89, 95]. 3.4 Image Reconstruction Image reconstruction is the process to reconstruct 2D and 3D images from acquisition data of molecular imaging modalities. Reconstruction algorithms include both analytical ones, e.g. Filtered Back Projection (FBP) and iterative ones, e.g. Maximum Likelihood (ML). Analytical algorithms are computationally fast, especially when applied to full 3D image reconstruction using 3D-FBP. FBP is also used in dynamic 66 F. Gao and P. Shi PET image reconstruction, where it is believed to provide better quantitative accuracy with the extremely low count data sets. Iterative algorithms are currently the mainstream reconstruction algorithms, which use statistical assumptions and provide images of overall better qualities. Image reconstruction is one of the most important processes in data processing, other image-based processes and clinical diagnosis all depend on the accuracy of reconstructed images. 3.5 Dynamic PET Analysis 3.5.1 Clinical Requirements The accuracy of quantitative dynamic PET studies depends on various factors including kinetic models, quantitative methods and the approximation of arterial input function from blood sampling. The most general kinetic models used are compartment model with assumptions that physiological process and molecular interactions are not influenced by injected radioligand. Current clinical adopted quantitative methods are actually semi-quantitative methods, which include methods using reference regions or calculating Standard Uptake Value (SUV). Methods using reference regions are easy to implement but have several drawbacks, e.g. the reference tissue is hard to define and has low SNR due to the low resolution of PET and SPECT scans, and the uptake of the reference tissue may change after the radiotherapy. SUV now is included in every clinical study, which is calculated as a ratio of tissue radioactivity concentration and injected dose divided by body weight, the advantage of SUV in clinical study is that the blood sampling is not required. However, the full quantitative analysis requires both dynamic PET scans and tracer concentration in the arterial blood plasma. The gold standard of blood sampling is serial arterial sampling of a superficial artery, and clinical alternative methods include venous blood sampling, image derived input function and population based input function. The drawback of the full quantitative method is only one FOV/bed position can be taken into consideration. For metastasized disease, not all lesions can be quantified simultaneously. 3.5.2 Compartment Models Compartment models are used in many fields including pharmacokinetics, biology, engineering etc. Compartment models are the type of mathematical models to describe the way materials (radiotracers and their metabolite in PET and SPECT scan) are transmitted among the compartments (different organs and tissues). Inside each compartment, the concentration of radiotracers is assumed to be uniformly equal. Due to their simplicity and plausibility, compartment models are widely used in the dynamic PET scans to describe the tracer/drug kinetics. Drug kinetic models include simple drug transport model, which generally contains equal or less than three compartments and can be solved directly, and compli- Shape Analysis in Molecular Imaging 67 Fig. 8 Illustration of six basic compartment models cated biological models, which can contain up to twenty compartments and generally require prior knowledge to solve [30, 31]. Most of the complicated models with many compartments can usually be decomposed into a combination of simple models with less than four compartments. The most basic compartment models used in kinetic analysis shown in Fig. 8 include two compartment blood flow model (Model 1), standard two tissue three compartment Phelps 4K model with reversible target tissue (Model 2) and Sokoloff 3K model with irreversible target tissue (Model 3), three tissue five parameter bertoldo model (Model 4), standard three tissue four compartment model (Model 5 and Model 6). More complicated models with more compartments and parallel model with multiple injection can be extended from aforementioned standard models [24, 44]. All six models can be represented by a set of differential equations with corresponding kinetic parameters K = {k1 , k2 · · · kn }, where n is the number of kinetic parameters. Here we utilize Model 2 as an example for demonstration. Model 2 can be represented by first-order differential equations dC F (t) = k1 (t)C P (t) + k4 (t)C B (t) − (k2 (t) + k3 (t))C F (t) dt dC B (t) = k3 (t)C F (t) − k4 (t)C B (t) dt (2) (3) The measurement of dPET is the combination of radiotracer in plasma C P , nonspecific binded radiotracer C F and specific binded radiotracer C B through C P E T = (1 − Vb ) · (C F + C B ) + Vb · C P Y = DC P E T + e (4) (5) where Vb is the blood volume fraction, Y is measured projection data , D is the system probability matrix, and e is the noises during acquisition. Equation (5) can be represented by a more general time-dependent form for all models as 68 F. Gao and P. Shi Y (t) = D X (K , t) + e(t) (6) Simple models with less than three compartments generally can be solved directly, while more complicated models need simplifications and various numerical approximations. 4 Review of Recent Advancements of Shape Analysis 4.1 Mathematical Modeling and Statistical Formulation of PET Image Reconstruction 4.1.1 Statistical Image Reconstruction Criterion The goal of mathematical modeling of data acquisition is to describe the transforms from spatial distribution of imaging objects to projection distributions on detector pairs in PET system. Denoting the spatial distribution of imaging object by a set of spatial variables x = {xi |i = 1 · · · n} ∈ Rn , where n is the total number of voxels, and the expected values (means) of projection bins by ȳ = { ȳ j | j = 1 · · · m} ∈ Rm , where m is the total number of bins, a mathematical expression of the transform can be obtained ȳ = Dx + ē (7) where D is the system response model giving the probability matrix of mapping the transform from x to ȳ, and ē is the means of background noises. A block diagram of the above procedure is shown in Fig. 9a. With the mathematical model specified above, the problem of PET image reconstruction is to find an estimation of the imaging object x̂ from the measurement data (projection bins) y. The image reconstruction is an ill-posed inverse problem, so an intuitive solution is to apply statistical assumptions of measurement data as regularizations. The statistical formulation tries to find an optimized relationship between measurement data y and imaging object x (or expected values of projection bins ȳ ) by defining different objective functions. If denoting yr ndsa is the measurement data y after all data corrections and assuming all the data correction are perfectly applied, ȳ is equal to the yr ndsa after all data corrections, however, in real clinical data acquisition, there will be some difference between ȳ and y, which is accumulated from residuals of each data correction. A block diagram of statistical image reconstruction framework is shown in Fig. 9b: the statistical properties of system noises or measurement data are first modeled based on certain statistical distributions (Gaussian, Poisson, their combination or other derivations) in block M, then inputted into system block D; the system output is generated and compared with corrected measurement data yr ndsa based on the predefined criteria, when the system goes convergent, estimations of the imaging Shape Analysis in Molecular Imaging 69 Fig. 9 Block diagrams. a PET data acquisition; b Statistical model based iterative PET reconstruction; c Designed system for PET image reconstruction; d Simplified block diagram (black box) of (c) object x̂ will be obtained. Three major criteria can be used are Least Square based (LS), Maximum Likelihood based (ML), Maximum A Posterior based (MAP). LS based methods try to obtain the estimation x̂ by minimizing the difference of fit between the predicted data by means of the modeling of the acquisition process and the measured data yr ndsa , the general objective function is x̂ = arg min ||yr ndsa − Dx||22 x (8) Extended algorithms are all based on the above equation and try to improve the performance by introducing different weights or penalization items. ML based methods try to obtain the estimation x̂ by maximizing the likelihood functions which represent the goodness of fitting statistical assumptions to measurement data. The general objective function is based on the conditional probability density of the image object with known measurement data yr ndsa as x̂ = arg max p(x|yr ndsa ) x (9) where p represents the probability density, and the above equation is a function of x. When a statistical model of measurement data or noises is defined, either Gaussian/Poisson or their combination, the objective likelihood function can be solved correspondingly. An improved method based on poisson assumption widely used is EM-SP (Shifted Possion) algorithm, which introduces two times the mean of randoms coincidences to the random precorrected data and models this sum as a Poisson random variable, and the objective function is still based on the above one [63]. 70 F. Gao and P. Shi Furthermore, a priori information can also be introduced in the form of statistical properties of the imaging object as constraints into the ML objective function. Bayesian formulation is such a probabilistic approach, statistical information of the unknown imaging object x is introduced by adopting the probability density of x, which is the prior p(x). Equation (9) can then be solved as follow: p(x, yr ndsa ) = p(x|yr ndsa ) = p(yr ndsa )  p(yr ndsa |x) p(yr ndsa ) p(yr ndsa |x) p(x) p(yr ndsa ) without prior with prior (10) where, p(x, yr ndsa ) is the joint probability density of x and yr ndsa . The expression with prior is just the objective function of MAP based algorithms, and different priors (e.g. image priors, independent priors) will lead to different implementations based on the above objective function. All above objective functions demonstrate the implied statistical knowledge assumptions on measurement data or noises in mathematical modelings, however, as discussed in the section of data corrections, the statistical distribution of corrected measurement data yr ndsa is neither single Gaussian/Poisson distribution nor their combination after various data corrections, and there will be further uncertainties during modelings. Due to the individual differences of patients during real PET scanning, the uncertainties during modeling will be more serious for over-weighted patients. With the development of new scintillators, PMTs and full 3D PET scanner, various system uncertainties and procedures of photo counting become more complicated. Additionally, there will be more scatter coincidences in measurement data. All above effects make the accurate modeling of PET measurement data an challenging issue, and it is almost impossible to establish a statistical model which can accommodate all the data errors. 4.1.2 Minimax Criterion To deal with the differences between ȳ and y, another solution is to adopt minimax criterion to reconstruct the measurement data yr ndsa from the point view of system energy. As shown in Fig. 9c, a new block F is designed to reconstruct measurement data yr ndsa , the output of the new system is designed to be the difference between expected values of imaging object x and estimation values x̂, a block P contains system response D and the input is all residuals and uncertainties e. Measurement data yr ndsa is constant during image reconstruction, so the system can be simplified with inner loop of block F and yr ndsa incorporated into block P and convert to a black box J as shown in Fig. 9d. When the output (difference between current value and estimated value) goes steady, the system J reaches convergent. With this system, the minimax criterion can be adopted intuitively, the PET image reconstruction procedure becomes to minimize the output (estimation errors x − x̂) with maximized input (disturbance e) through system J : x̂ = min max J x−x̂ e (11) Shape Analysis in Molecular Imaging 71 Fig. 10 Illustration of system gains and predefined upper bound From the point view of transition system in engineering, system J can be treated as a transformation from all residuals and uncertainties e to estimation errors x − x̂. In order to obtain a desired output, a proper upper bound for system gain of system J must be defined first, which means the transfer function of J shall also be required to comfort to the same predefined upper bound. Denoting the upper bound be γ 2 , the previous description is equal to max ||J ||2 < γ 2 (12) The ∞-norm of the system can be further interpreted as the peak system gain, so if the ∞-norm of the system satisfies ||J ||2∞ < γ 2 , then all the system gains will be less than γ 2 as illustrated in Fig. 10. Here J is an element of the Hardy space, whose members consist of all stable, causal, transfer functions [49]. The continuous form of the ∞-norm ||J ||∞ will be ||J ||2∞ = sup ||e||2 =0 ||x − x̂||22 ||e||22 = sup σ̄ 2 (J ) (13) ω where sup stands for supremum, and σ̄ (J ) is the maximum singular value of J . From Eq. (13), it is easy to obtain max ||J || ≤ ||J ||2∞ = sup ||e||2 =0 ||x − x̂||22 ||e||22 < γ2 (14) This criterion represents a family of solutions where the peak system gain of J is less than the predefined upper bound γ 2 . The problem is transformed from seeking the maximum system gain of J to calculating the ∞-norm of system gain with upper bound γ 2 . With both Eqs. (11) and (14), in order to calculate the supremum, intuitively we can minimize the numerator (estimation errors x − x̂) while maximizing denominator (system uncertainties e). Then the objective function for the problem will be min max ||x − x̂||22 − γ 2 ||e||22 < 0 x−x̂ e (15) 72 F. Gao and P. Shi Seeking solutions under disturbances is a difficult problem. By introducing the ∞-norm of system gain, solutions with peak system gain will perform well under any circumstance and make the problem globally optimized. Furthermore, the minimax criterion allows one to identify a robust solution that has the best worstcase performance. The robustness of the estimator arises from the fact that it yields an energy gain less than γ 2 for all bounded energy disturbances no matter what they are. 4.1.3 Static PET Image Reconstruction Under Minimax Criterion PET images reflect the concentration and metabolism of radiotracers in vivo, and the metabolic rates of radiotracers in organs or tissues may vary with time, which can be represented by differential equations with different constraint models as ẋ(t) = C(x, t), where C(x, t) models the concentrations of radiotracers with time. In this section, we focus on the static PET image reconstruction and analysis, which assumes that the distribution of the radiotracers in the body is temporally stationary corresponding to the autoradiographic (ARG) model. The ARG assumes the equilibrium of the metabolic ratio, which can be represented by the first-order differential equation of x as: ẋ(t) = 0 (16) Corresponding discrete-time form will be x ′ (m) = H (m)x(m) + v(m) (17) where m represents the iterations in one discretized time frame, H is the transition matrix representing the update of estimations, and v is the possible uncertainties in measurement. Similarly, the discrete-time form of PET measurement equation will be yr ndsa = Dx(m) + e(m) (18) Then the transition system J is extended to include uncertainties v as input together with e. The initialization of x can also affect the reconstruct accuracy and speed, so it is also included in the input as a source of uncertainties. Finally, the minimax criterion can be derived from Eq. (14) as ||J ||2∞ = sup x(0) − x̂(0)2 −1 p0 x(m) − x̂(m)2Z (m)  + m (v(m)2 −1 + e(m)2  m V (m) E(m)−1 ) (19) Z (m), p0−1 , V (m)−1 and E(m)−1 are the weighting matrices at iteration m to make the criterion more extensible. After setting the upper bound γ 2 , the final objective function of minimax criterion for PET image reconstruction can be derived based Shape Analysis in Molecular Imaging 73 on Eqs. (15) and (19), which minimizes the estimation errors z − ẑ with maximized initial disturbance x(0) − x̂(0), measurement disturbance e and v: min max z(m)−ẑ(m) v,e,x(0)−x̂(0) ||J ||2 =  m x(m) − x̂(m)2Z (m) − γ 2 x(0) − x̂(0)2p−1 o −γ 2  (v(m)2V (m)−1 + e(m)2E(m)−1(20) ) m 4.1.4 Solutions and Optimization The final objective function Eq. (20) of minimax criterion for static PET image reconstruction is established from the point view of transition system in engineering. The objective function calculates the system energy instead of building sophisticated statistical assumption of measurement data, and this makes the minimax criterion more robust. Of many methods that can be adopted to solve the above objective function, the well-validated H∞ filter is the best to optimize this ∞-norm problem [49, 79]. A specific H∞ filter has been modified based on Eq. (20). The H∞ filter represents a typical minimax problem where the worst situation is first induced by the disturbances and then the estimator is introduced for improvement, in other words, the H∞ filter is in fact a two-person game between the external disturbances and the estimator. This solution is just like optimization by using a game theoretic algorithm which can be implemented through recursive updating of the filter gain K (m), the Riccati difference equation solution P(m), and the state estimates x̂(m) as follow : K (m) = H (m)P(m)S(m)D T E(m)−1 P(m + 1) = H (m)P(m)S(m)H (m)T + V (m) x̂(m + 1) = H (m)x̂(m) + K (m)(yr ndsa − D x̂(m)) (21) (22) (23) S(m) = (I − γ −2 Z (m)P(m) + D T E(m)−1 D P(m))−1 P(0) = p0 where H∞ gain K (m) indicates the system gain (correspondingly shows convergence of the state estimation). Convergent estimations of the imaging object will be obtained when the gain K (m) goes steadily. m still represents the number of iterations. From the solution procedure, it can be noticed that the H2 norm filter for this objective function is just the Kalman filter widely used. Detailed proofs of above solution can be found in [79]. 74 F. Gao and P. Shi 4.1.5 Computation Issues Design of H∞ filter consists of choosing the weighting matrices Z , E, V , po and the performance bound γ 2 . When there are accurate modelings of some effects of PET acquisition, corresponding E, V , po can be initialized by the covariance matrices of measurement disturbance e, state transition disturbance v and initial condition x̂(0). If there is no modeling or one does not want to use current models, E, V , po can be simply initialized by identity matrix. Since we generally assume the estimation results after convergence are just what we desired, so Z is initialized by identity matrix, when scanning over-weighted patient and suffering obvious underestimation, weighting matrix Z can be scaled by a scaler less than 1. γ 2 is the predefined upper bound of performance, theoretically, the smaller the γ 2 value, the smaller the estimation error, however, the selection of γ 2 must make the Riccati equation have a positive definite solution. So firstly, we define and iteratively update a residual matrix R(m + 1)−1 through R(0) = (P(0)−1 − γ −2 Q(0))−1 (24) −1 R(m + 1)−1 = [H (m)(R(m)−1 + D T E(m)−1 D) H (m)T + V (m)]−1 − γ −2 Z (m) (25) As a result, the optimal γ value can be determined as: [H (m)(R(m)−1 + D T E(m)−1 D) −1 H (m)T + V (m)]−1 − γ −2 Z (m) > 0 → [H (m)(R(m)−1 + D T E(m)−1 D) −1 H (m)T + V (m)]−1 > γ −2 I → γ 2 I > H (m)(R(m)−1 + D T E(m)−1 D) −1 H (m)T + V (m) → γ 2 > max{eig[H (m)(R(m)−1 + D T E(m)−1 D) −1 H (m)T + V (m)]} 0.5 → γ = ξ max eig[H (m)(R(m)−1 + D T E(m)−1 D)−1 H (m)T + V (m)] (26) where max{eig(A)} denotes the maximum eigenvalue of the matrix A, and ξ is a constant larger than 1 to ensure that γ is always greater than certain optimal performance level. If the γ value is too close to the optimal performance level, i.e. ξ ≈ 1, it might lead to numerical errors because the matrix R(m) is now close to a singular matrix. The matrix inverse is required in every time step in the conventional H∞ filter in order to calculate the H∞ gain. Generally, inversion of small matrices is fairly easy, but the inversion of a large matrix will require more computational costs in a practical implementation. The steady state H∞ filter is much easier to be implemented in a system in which real-time computational effort or code size is a serious consideration [80]. The minimax objective function is given as Eq. (20), where the parameter matrices N , V and Q are symmetric positive definite matrices chosen by the designer based Shape Analysis in Molecular Imaging 75 on the specific problem. Since the designed parameters of the underlying system can be treated as fixed values for input, then the steady state solution to the minimax problem can be obtained. Referring to H∞ filter, the steady state solution will be K = P S D T N −1 P = H PSHT + V x̂(m + 1) = H x̂(m) + H K (m)(y(m) − D x̂(m)) (27) (28) (29) S = (I − γ −2 Q̄ P + D T N −1 D P)−1 Q̄ = F T QF In order to have a solution to the problem, the following condition must be hold: P −1 − γ −2 Q̄ + D T N −1 D > 0 (30) If γ −2 , F , N or Q is too large, or D is too small, the H∞ estimator will have no solution. After the conditions above satisfied, Eq. (28) can be written as P = H [P −1 − γ −2 Q̄ + D T N −1 D]−1 H T + V (31) Applying the matrix inversion lemma to the inverse of the above expression, we can get P = H P − P[P + (−γ −2 Q̄ + D T N −1 D)−1 ]−1 P H T + V = H P H T − H P[P + (−γ −2 Q̄ + D T N −1 D)−1 ]−1 P H T + V (32) Equation (32) is a discrete-time algebraic Riccati equation that can be solved by control system software or numerically iterating Eq. (28) until it converges to a steady state value. The disadvantage of the steady state H∞ filter is that theoretically it does not perform as well as the time-varying filter. However, the reduced performance that is seen in the steady state H∞ solution is often a small fraction of the optimal performance, whereas the computational saving can be significant [80]. 4.2 Dynamic PET Image Reconstruction 4.2.1 Radioisotope Decay Constrained Dynamic PET Imaging Both PET and SPECT use radiotracers, which decay with time. The natural decay property of the radioisotope, can be introduced into the objective function as the temporal guidance for multi-frame image sequence reconstruction, and can also be used to separate multiple radiotracers in dynamic imaging. The projection equation 76 F. Gao and P. Shi of dynamic PET imaging can be formulated through an affine transform between the projection data and emission object as: y(t) = Dx(t) + r (t) + s(t) (33) where the emission sinogram data is represented by a vector y, and the activity of emission object is represented by x. D is system probability matrix, which gives the probability of a photon emitted from ith voxel being detected in projection jth bin. t is the time frame. r and s are the contribution of random coincidence events and scatter coincidence events. After the conventional online delayed-window random correction, the Eq. (33) can be rewritten as: y(t) = Dx(t) + e(t) (34) here e is an error vector, which represents unknown measurement uncertainties including scatter coincidence events. In the section of static PET imaging reconstruction, the distribution of the radioisotopes in the body is assumed to be temporally stationary corresponding to the autoradiographic model, however, in the real situation, the radioisotope will decay with time, and its activity at time t should be x = X 0e ln(0.5) T t (35) here X 0 is the initial activity distribution, and T is the half life of the radioisotope. So the real-time change of radioisotope can be represented as dx ln(0.5) ln(0.5) t = X0 e T dt T (36) then the dynamic change of radioisotope from one frame to the next can be obtained from the integral of Eq. (36). A general representation of state transition will be x(t + 1) = H (t)x(t) + v(t) (37) where x(t) is the radioactivity concentration at time frame t, and H (t) is a coefficient matrix for state transition at time frame t. v(t) represents the uncertainties during state transition. With the introduction of decay model shown as Eq. (37), we are able to make use of the radioisotope’s own temporal properties as constraints to guide our reconstruction. 4.2.2 Dynamic PET Image Reconstruction with Minimax Criterion Minimax criterion, which allows one to identify a robust solution as one that has the best worstcase performance can also be applied to dynamic PET reconstruction. In Shape Analysis in Molecular Imaging 77 general, a robust discrete optimization problem can be formulated as follows. Let X be the set of all solutions, E be the set of uncertainties of measurement in single time frame, and M be the set of uncertainties for state transition among time frames, performance of a solution x ∈ X under uncertainties e ∈ E and v ∈ M is F(x, e, v). The problem is to find the solution that has the best worst-case performance, which is the same as minimizing (over all solutions) the maximum (over all uncertainties) performance: (38) min max F(x, e, v) x∈X e∈E,v∈M from the description of Eqs. (34) and (37), the estimation of activity distribution x(t) at time t is not only computed based on measurement y(t), but also affected by previous estimations, so we define a linear combination of x(t) as z(t) = F x(t) = g(x(k), H (k), v(k)) k = 1, 2...t (39) so the measure of performance F(x, e, v) is given by J= x(0) − x̂(0)2 −1 po z(t) − ẑ(t)2Q(t)  + (v(t)2 −1 + e(t)2  V (t) N (t)−1 ) (40) where the notation x2G is defined as the square of the weighted (by G) L 2 norm of x (i.e. x2G = x T Gx). N (t), V (t), Q(t) and po are the weighting matrices for the uncertainties of measurement in single time frame, the uncertainties of state transition among time frames, the estimation error, and the initial conditions respectively. x̂(0) is the initial estimate of the state. The optimal estimate z(t) among all possible ẑ(t) should satisfy: J ∞ = sup J < γ 2 (41) where γ 2 > 0 is a prescribed level of disturbances. It is assumed that the L 2 norms of e(t) and v(t) exist. Then the minimax performance criterion of Eq. (40) where the estimator strategy z(t) playing against the exogenous inputs e(t), v(t) and the initial state x(0) becomes min max J = z(t)−ẑ(t) v,e,x(0)  z(t) − ẑ(t)2Q(t) − γ 2 x(0) − x̂(0)2p−1 o  (v(t)2V (t)−1 + e(t)2N (t)−1 ) −γ 2 (42) Now the problem becomes to solve the above objective function, and the kinect model (e.g. the decay model in this section) is successfully incorporated in it. Same solution framework as static PET image reconstruction can be used here. 78 F. Gao and P. Shi 4.3 Kinetic Parameter Estimation From the section of compartment models, a two-tissue three-compartment model (Model2) is general enough to describe regional tracer kinetics as shown in Fig. 8, where C P (pmol/ml) is arterial concentration of radiotracer, C F and C B (pmol/ml) are the concentrations of non–specific binding and specific binding tracers in tissues. Parameters k1 , k2 , k3 and k4 (min −1 ) specify radiotracer transport rates. The time variation of kinetic model in voxel i can be denoted by first-order differential equations as: dC Fi (t) = k1i (t)C Pi (t) + k4i (t)C Bi (t) − (k2i (t) + k3i (t))C Fi (t) dt dC Bi (t) = k3i (t)C Fi (t) − k4i (t)C Bi (t) dt (43) (44) Here this model will be used to derive the solutions framework of kinetic parameter estimation. 4.3.1 Modeling of Dynamic PET Measurement with Tracer Kinetics Dynamic PET imaging involves a sequence of contiguous acquisition with different temporal resolutions, which can be formulated as a projection transform: y(t) = Dx(t) + e(t) (45) Here, y(t) is the projection data and x(t) = {xi (t)|i = 1, · · · , n}T is the activity concentration at time frame t. n is the total number of voxels. D is the system probability matrix. e(t) is the overall measurement uncertainties. Here we will transform Eq. (45) to accommodate kinetic models. Firstly, activity concentration x will be the combination of C F and C B , then Eq. (45) will be y(t) = D D  C F (t) + e(t) C B (t) (46) where C F (t) = {C Fi (t)|i = 1, · · · , n}T and C B (t) = [C Bi (t)|i = 1, · · · , n}T . After the dynamic change of measurement dydti (t) being deduced, we substitute the differential equations Eqs. (2) and (3) and do a simple transformation to arrange four kinetic parameters (k1 , k2 , k3 , k4 ) in a column vector will yield Shape Analysis in Molecular Imaging ⎡ .. . 79 ⎤ ⎢ dC Fi (t) ⎥ ⎥ ⎢ ⎢ dt ⎥ dyi (t) ′ ⎢ .. ⎥ = D D ⎢ . ⎥ + ei (t) ⎥ ⎢ dt ⎢ dC Bi (t) ⎥ ⎣ dt ⎦ .. . ⎡ ⎤ .. . ⎤ ⎢ ⎥⎡ ⎢ C Pi (t) −C Fi (t) −C Fi (t) C Bi (t) ⎥ k1i (t) ⎢ ⎥⎢ ⎢ ⎥ ⎢ k2i (t) ⎥ .. ⎥ +e′ (t)(47) = D D ⎢ ⎥⎣ . ⎢ ⎥ k3i (t) ⎦ i ⎢ 0 0 C Fi (t) −C Bi (t) ⎥ ⎣ ⎦ k4i (t) .. . ⎡ .. . ⎤ ⎡ ⎤ ⎥ ⎢ k1i (t) ⎢ C Pi (t) −C Fi (t) −C Fi (t) C Bi (t) ⎥ ⎥ ⎢ ⎢ k2i (t) ⎥ ⎥ ⎢ .. ⎥ By denoting Ri (t)= ⎢ ⎥ and Si (t)= ⎢ . ⎣ k3i (t) ⎦, ⎥ ⎢ ⎢ 0 0 C Fi (t) −C Bi (t) ⎥ k4i (t) ⎦ ⎣ .. . we can get the dynamic change of total measurement data from all voxels as n dy(t)  dyi (t) = = D D dt dt i=1 ⎡ ⎤ S1 (t) ⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎥ ′ ⎥ R1 (t) · · · Ri (t) · · · Rn (t) ⎢ ⎢ Si (t) ⎥ + e (t) ⎢ .. ⎥ ⎣ . ⎦ Sn (t) (48) ′ = D D R(t)S(t) + e (t) Now we have set up the relationship between the change of measurement data and kinetic parameters directly by Eq. (48). 4.3.2 Solution Under Minimax Criterion The minimax solution framework can also be transform to estimate kinetic parameter. No statistical assumptions needed makes the minimax criterion robust to the poor statistical properties in low count acquisition and system noises. Since kinetic parameters are generally assumed to be constant, we can set: 80 F. Gao and P. Shi S(t + 1) = S(t) + v(t) (49) Here, v is possible disturbances. With Eqs. (48) and (49), the corresponding minimax performance equation will be “min S∈L maxe∈E,v∈V F(S, e, v)”, where L, E and V are the sets of solutions, uncertainties of measurement and state transition. As an iterative solution, we also define a linear combination of S(t) as “z(t) = g(S(m), v(m)) where m = 1, 2...t”, then objective function J will be J= S(0) − Ŝ(0)2 −1 po z(t) − ẑ(t)2Q(t)  + (v(t)2 −1 + e(t)2  V (t) E(t)−1 ) (50) where the notation x2G is defined as the square of the weighted (by G) L 2 norm of x. po , E(t), V (t) and Q(t) are weighting matrices. Ŝ(0) is the initialization of x. More detailed settings and initialization of parameters can be found in [80]. 5 Clinical Practices and Future Directions 5.1 Personalized Drug Metabolism Analysis Dynamic PET is a molecular imaging technique that is used to monitor the spatiotemporal distribution of a radiotracer in vivo and enables cellular level metabolism analysis in clinical routine. dPET provides a good promise for quantitative lesion metabolism analysis to help identify lesions. However, due to poor statistical properties of the measurement data in low count dynamic PET acquisition and disturbances from surrounding tissues, identifying small lesions inside the human body is still a challenging issue. Furthermore, the mismatch between general purpose models and patient size/motions makes the situation even worse. Quantitative kinetic analysis of radiotracer uptakes requires the reconstruction of kinetic parameters [25, 77, 96]. The mainstream is statistical reconstruction algorithms, however, whose quality is determined by the accuracy of sophisticated system probability matrix (SM). Many efforts have been devoted to improve the accuracy of SM [64, 70, 85, 107]. However, the ideal SM is almost impossible to obtain under practical conditions. The general purpose SM also could not compensate different sizes of patients and the motions during acquisition, which will decrease the accuracy of reconstructions. Furthermore, the reconstruction of dynamic PET image sequences, whose poor temporal resolution, insufficient photon counts, more complicated data corrections and poor statistical properties of measurement data also requires a more accurate SM. To improve the personalized lesion metabolism analysis, we can generate a patient adaptive SM using machine learning techniques. Both patient size information and potential small lesion information are incorporated by simulations of phantoms of different sizes and individual point source responses [50, 53, 102]. Training data Shape Analysis in Molecular Imaging 81 set from simulations of 90 studies is conducted using 15 phantoms of different sizes based on Zubal thorax phantom. Each simulation has randomly generated motions and lesions in lung. The personalized SM should be able to differentiate true lesion and false lesion, and further deal with input functions of different accuracies. 5.1.1 System Matrix Derived from Supervised Learning Statistical reconstruction requires a well modeled SM, which directly determines the accuracy of reconstruction results. The SM D is extended to include 2 parts, D1 is a SM generated from geometry information and physical phenomena, and will account for sizes and motions of different patients, D2 is an additional SM generated from point source responses. D1 and D2 are full size SM, and combined together by weighting matrices w1 and w2 according to the anatomical information of patients. This effort makes the SM more patient adaptive. The measurement equation is extended to be  D1 X (κ, t) + e(t) (51) Y (t) = w1 w2 D2 w1 , w2 , D1 , D2 are updated by supervised learning. Training sets are provided by Monte Carlo simulations using GATE toolbox [41]. Correspondingly, two series of simulations are performed, one is performed with human thorax phantom of different sizes, and the other is done by point source response inside a thorax phantom of normal size. Denoting the activity concentrations as X = {x1 , x2 , · · · xn } and measurement datasets as Y = {y1 , y2 , · · · xn }. n is the number of training sets, and every dataset is a dynamic data sequence related to time t. For simplification of expression, the PET measurement equation is written as Y (t) = D ′ X (k, t) + e(t). Since ADALINE has been proved to be simple yet successful for updating SM in [85], we also adopt ADALINE for our SM training here. The initialization of D1 and D2 are the SMs generated with uniform cylindrical phantom. The update procedure by ADALINE using back-propagation and least mean square error is: ′ ŷm (t) = Dm X (k, t) + em (t) ′ Dm+1 (t) ′ Dm (t) + 2Lδm X T (t) = δk (t) = Y (t) − ŷm (t) em+1 (t) = em (t) + 2Lδm (t) (52) (53) where m is the iteration step of training, and L is the learning rate. After defining a precision level of learning ε, D ′ subject to  Y (t) − ŷm (t) < ε ŷm (t) − Y (t) < ε (54) the weighting matrices w1 and w2 will be obtained when convergence is achieved (Fig. 11). 82 F. Gao and P. Shi Fig. 11 Demonstration of five phantoms from the phantom library used to generate personalized SM 5.1.2 Parameter Reconstruction of Dynamic PET The kinetic model and image reconstruction are combined in one equation, the log likelihood function can be derived with measurement data y(t) as L (y|κ) =  y(t) log ȳ(κ, t) − ȳ(κ, t) (55) t where ȳ(κ, t) = D ′ x(κ, t) + e(t) (56) κ̂ = arg max Φ(κ) (57) Φ(κ) = L (y|k) − βU (k) (58) where U is the penalty regularization term with parameter β controlling resolution/noise tradeoff. Equation (57) can be solved by a paraboloidal surrogates algorithm in [22]. The above solution is global convergent, however, the kinetic parameter reconstruction has a higher data dimensionality/freedom, so we also define the evaluation of kinetic parameters through images by using student’s t-distribution hypothesis test to determine their statistical differences among iterations. By selecting Region of Interest (ROI), calculate |x̄m − x̄m+1 | (59) t= σ where σ =  and varm + varm+1 − 2covm,m+1 N 0.5 (60) N covm,m+1 = 1  (xm,i − x̄m )(xm+1,i − x̄m+1 ) N −1 (61) i=1 x̄m and x̄m+1 are the means in ROI at iteration m and m + 1, var is the corresponding variances across the image elements. cov is the covariance across the two iterations. t is calculated until less than t0.05 in the t-table to show a confidence level of 95 % that the difference between images is small enough. Shape Analysis in Molecular Imaging 83 (b) (a) Fig. 12 a Lesion in 40th slice; b 32nd slice as reference Table 5 Summary of experiments Lesion type Input function Successful estimation/Total studies Generic SM Personalized SM Group1 Group2 Group3 Group4 Group5 Group6 True 1 True 2 True 3 False 1 False 2 False 3 13/15 14/15 7/15 13/15 5/15 7/15 11/15 12/15 7/15 12/15 6/15 7/15 5.1.3 Clinical Results The clinical patient data in this study was a dynamic PET scan acquired from a 28-year-old, 75 kg male volunteer. 10 mCi 18 F-FDG was injected and a dynamic acquisition of the thoracic cavity started just after injection. The acquisition consists of 40 time frames: 20 × 0.5 min, 15 × 1 min, and 5 × 2 min. All corrections are performed properly with the software provided by the scanner. The input function is estimated by the image-derived method. Figure 12a shows a lesion region by a red arrow in the 40th slice. We calculate the influx rates of the lesion and compare them with the heart muscles in the 32nd slice. The lesion metabolism calculated by the new personalized SM is closer to the muscles than that by generic SM, and the lesion is confirmed by the doctor as a false lesion with temporarily increased metabolism than muscles. Results from the new personalized SM show potential improvements in diagnosis. Table 5 shows the statistical results from all 90 studies, which shows the ratio of successful identification of lesions using different input functions. Input Function 1 is perfect input function (equivalent to perfect blood sampling with less than 5 % error), Input Function 2 is good input function (equivalent to a disturbed blood sampling with about 20 % error), Input Function 3 is an Image Derived Input Function (IDIF). Both methods performs well by using Input Function 1. However, with Input Function 2 (with disturbances), the accuracy of generic SM decreases, but the personalized SM still provides good results. The results show the improvement in identifying lesions by personalized SM, and the reduction of requirement of accurate input function. 84 Table 6 Influx rate of patient study F. Gao and P. Shi Muscles Generic SM Personalized SM 0.0054 0.0085 0.0071 5.2 Model Selection In drug discovery and development, the procedure of drug selection is full of challenging issues. Kelloff et al. show that more than 90 % of all new oncology drugs fail in the late stages of development because of inadequate activity and difficulties in determining their efficacy [43]. Quantitative pharmacokinetic analysis with dynamic PET imaging now plays a promising role as determinants of in vivo drug action to help select drug candidate. Fast and accurate pharmacokinetic analysis with rapid information feedback in the early stage of drug discovery and development is critical to obtain the in vitro and in vivo drug properties [14, 100]. A typical procedure of pharmacokinetic analysis by dPET imaging includes, firstly, setting up a working hypothesis of the target enzyme or receptor for a particular disease, secondly, establishing suitable models (or surrogate markers) to test biological activities, and at last, screening the new drug molecules for biological activities. In this procedure, model selection by dPET has seldom been studied because of various scientific challenges, for example, (1) the kinetic models for drugs are generally very complicated, when facing a new biomarker (new drug), it is hard to determine which model will work best, (2) accurately solving these complicated models always needs special mathematical considerations, (3) although we can always use more complicated models to represent certain biological activity, the computational cost increase significantly due to the complex of the model, which cause a burden for the early drug discovery. (4) measurement data from dPET suffers from poor spatial and temporal resolutions, especially the first several time frames, (5) blood sampling is required in pharmacokinetic analysis but it is very hard to generate an accurate one [31]. 5.2.1 Temporal-Difference Reinforcement Learning Machine learning in image processing and analysis is growing rapidly [97]. Of various machine learning methods, reinforcement learning is meant to be a general approach to learn from interactions [86]. It is a control method which presents a robust mechanism for goal directed decision making. Unlike supervised learning methods, no examples of desired behaviors are provided during training, instead, behavior is guided through positive or negative reinforcements [99]. So this method do not require a large training dataset, and is especially suitable for preclinical drug selection and pharmacokinetic analysis with only limited data sets [84]. Additionally, as a control mechanism, the method can help solve the complicated kinetic models with noisy dPET acquisition data. At the same time, the method can inherently deal with disturbances during blood sampling. Therefore, in this section we will intro- Shape Analysis in Molecular Imaging 85 duce a reinforcement learning based method which combines model selection and parameter estimation for pharmacokinetic analysis by dPET. Temporal Difference reinforcement learning (TD Learning) is a combination of Monte Carlo ideas and dynamic programing [86, 99]. Therefore, like Monte Carlo methods, TD Learning can learn from raw experiences without pre-defined models, and like dynamic programing, TD learning can update its estimations based on a part of learning outcomes rather than the final outcome. These features are especially suitable for model selection and noisy dPET data. Regardless of model, when we have an initial K, we define an action set a, which contains 2n components, {a1+ , a1− , a2+ , a2− , · · · an+ , an− , }, with the subscript corresponds to the index of kinetic parameters, and the superscript represents increasing (+) or decreasing (−) that kinetic parameter by certain amount during estimation. We derive the TD Learning algorithm for model selection and parameter estimation from the classic one as shown in Algorithm 1. Details will be shown in next subsection. 5.2.2 Model Selection As shown in the algorithm, reinforcement learning acts according to the rewards, we define three rewards based on physical constraints for model selection. The combination of three rewards is able to exclude non-matching models fast, which can improve the computational efficiency, and reduce the disturbances from the noises in low count dPET data. By denoting K ′ as the estimated kinetic parameters after selecting an action from a, 1. Reward 1: We compare the measured total counts in each time frame of measurement data Y and estimated total counts with K ′ . M S E(T otalCounts(Y (t)), T otalCounts(D X (K ′ , t))) < T hr eshold1 (62) 2. Reward 2: We compare the first order difference of Time Activity Curve (TAC) from measurement data Y and TAC curve estimated with K ′ . M S E(Di f f er ence(T AC(Y )), Di f f er ence(T AC(D X (K ′ , t))) < T hr eshold2 (63) 3. Reward 3: This is an optional reward, if there is priori knowledge from clinical data available, they are learned together through a 2-hidden layer neural network (NN) as shown in Sect. 5.1.1, and then used as a reference for estimated K ′ . M S E(N N (Y ), K ′ ) < T hr eshold3 (64) 86 F. Gao and P. Shi Where M S E is the operation to calculate the mean squared error. When each reward criterion is met, we have a reward (r ew = +1), otherwise, (r ew = −1). Then we accumulate all rewards by eligible trace Q-Learning [86] as shown in the Algorithm1 Q(K , a) ← Q(K , a) + α[r ewm+1 + γ Q(K m+1 , am+1 ) − Q(K , a)] (65) Where α, γ are learning parameters, which control the width and depth of learning. Proofs in [86] show that α and γ mostly affect the convergence speed, and have only limited effect on learning accuracy after convergency. Q is the value function in the Q-Learning, which stores all the rewards, and m represents iteration steps. Then the algorithm is applied to every model in the model bank, with each estimated K ′ from the maximum in Q, we calculate the Bias between estimated TAC (T AC e ) and true TAC (T AC m ) from measurement data for each model by T ACm −T ACe  Bias = T1 , where T is number of time frames. The model with the T ACm lowest Bias in the model bank with be the selected model by the proposed method. 5.2.3 Parameter Estimation When using Algorithm 1 to choose model, simultaneously calculated K ′ will be the initial parameter for that model. And the kinetic parameter can also be further calculated with a refined action set ar e f containing smaller increasing or decreasing amount using Algorithm 1. Algorithm 1 Model Selection and Parameter Estimation by TD Learning Initialize Q(K,a) arbitrarily Initialize K Repeat until convergency Randomly choose one action from a Repeat for all steps Take action a, generate K′ , observe reward r ew using defined criteria Choose a′ from K′ derived from Q Accumulate all rewards using Eq. 65 K ← K′ End End Select the Maximum in Q, Corresponding K is the estimated K 5.2.4 Clinical Results We study three cases of real patient dPET scans. Figure 13 shows the three cases, the first case is the scan of patient thorax, a ROI is defined in the normal muscle region, Shape Analysis in Molecular Imaging (a) 87 (b) (c) Fig. 13 a Case 1; b Case 2; c Case 3 Table 7 Model selection results Case 1 Case 2 Case 3 Model1 Model2 Model3 Model4 Model5 Model6 1.5901 0.7963 1.3809 1.4248 NA 1.0721 1.1735 NA 1.6025 1.6549 0.9482 1.2169 1.5558 NA NA NA NA NA NA is “Not Applicable” Table 8 Estimated kinetic parameters for three cases Case 1 Case 2 Case 3 K1 K2 K3 K4 K5 K6 0.0960 0.1080 0.0400 0.1540 0.0838 0.0620 0.0375 NA 0.0145 NA NA 0.0015 NA NA NA NA NA NA NA is “Not Applicable” the second case studies a ROI in the heart region, and the third case is with a ROI in the liver. The dynamic PET scans are performed on our PET scanner, the dynamic data set consists of 40 time frames: 20 × 0.5 min, 15 × 1 min and 5 × 2 min. The input function is estimated by fitting the reconstructed dynamic images. This input function is equal to a disturbed one affected by noises in reconstructed images. The model bank used is the compartment models shown in Fig. 8. The model selection results is shown Table 7. The results of model selections are consistent with suggestions from clinical studies. ROI 1 is normal tissue and Model3 is mostly used. ROI 2 is near the left ventricular and highly affected by the blood flow, so the blood flow model (Model 1) is most suitable. For ROI 3, clinical results had shown the necessity and importance of estimation of K4 in liver cancer, and the proposed method correctly choose the right Model2 (Phelps 4K model). And the non-applicable models are excluded successfully. The estimated kinetic parameters are shown in Table 8. 5.3 Future Directions In this chapter, we review the shape analysis and application in molecular imaging. Molecular imaging is a relatively new but fast-developing area for both research and 88 F. Gao and P. Shi clinical applications. Although with some limitations, molecular imaging modalities especially PET show the superior ability to quantitatively measure the biologic processes in a functional way at the cellular and subcellular level within living subjects. All the data processing and analysis of molecular imaging depend on the physical natures of the molecular imaging modalities. The emerging new scanner systems with new detectors will further enhance their abilities, and bring new challenges in data correction and image analysis at the same time. The data corrections algorithms need to be adjusted with the properties of new system design, and new features in detector system correspondingly. Monte Carlo simulation is always a good way to study the new design and provide references for validation of new data correction methods. And image reconstruction with disease specified statistical model instead of the generic models can improve the image qualities of certain disease. With the advancement of data correction and image reconstruction methods, the image post-processing including image registration and image fusion must adopt related changes. Researchers are also actively using machine learning methods to extract applicable pathological models from a series of patient studies and apply the strategy to personalized treatment. Pharmaceutical companies are also interested in the accurate quantitative pharmacokinetic parameter estimation using PET to study the metabolism of new drugs, which has the potential to shorten the drug development cycle and save tons of money for the industry and patients. With the evolution of both image pre-processing and post-processing methods, molecular imaging is believed to be able to study more complicated diseases currently in the unknown area, for example, the origins of Alzheimer’s disease and dementia. References 1. Bailey D, Townsend D, Valk P, Maisey M (2005) Positron emission tomography: basic sciences. Springer, London 2. Bailey DL, Meikle SR (1994) A convolution-subtraction scatter correction method for 3D PET. Phys Med Biol 39(3):411 3. Bao Q, Newport D, Chen M, Stout DB, Chatziioannou AF (2009) Performance evaluation of the inveon dedicated PET preclinical tomograph based on the NEMA NU-4 standards. J Nucl Med 50:401–408 4. Bartenstein P, Minoshima S, Hirsch C, Buch K, Willoch F, Mösch D, Schad D, Schwaiger M, Kurz A et al (1997) Quantitative assessment of cerebral blood flow in patients with Alzheimer’s disease by SPECT. J Nucl Med 38(7):1095 5. Bateman TM, Heller GV, McGhie AI, Friedman JD, Case JA, Bryngelson JR, Hertenstein GK, Moutray KL, Reid K, Cullom SJ (2006) Diagnostic accuracy of rest/stress ECG-gated Rb-82 myocardial perfusion PET: comparison with ECG-gated Tc-99m sestamibi SPECT. J Nucl Cardiol 13(1):24–33 6. Bendriem B, Townsend DW (1998) Theor Pract 3D PET. Kluwer Academic/Plenum Publishers, New York 7. Bentourkia M, Msaki P, Cadorette J, Lecomte R (1995) Assessment of scatter components in high-resolution PET: correction by nonstationary convolution subtraction. J Nucl Med 36(1):121–130 Shape Analysis in Molecular Imaging 89 8. Beyer T, Tellmann L, Nickel I, Pietrzyk U (2005) On the use of positioning aids to reduce misregistration in the head and neck in whole-body PET/CT studies. J Nucl Med 46(4):596– 602 9. Bonte FJ, Weiner MF, Bigio EH, White C et al (1997) Brain blood flow in the dementias: SPECT with histopathologic correlation in 54 patients. Radiology 202(3):793–797 10. Boucher L, Rodrigue S, Lecomte R, Bénard F (2004) Respiratory gating for 3-dimensional PET of the thorax: feasibility and initial results. J Nucl Med 45(2):214–219 11. Brankov JG, Galatsanos NP, Yang Y, Wernick MN (2003) Segmentation of dynamic PET or fMRI images based on a similarity metric. IEEE Trans Nucl Sci 50(5):1410–1414 12. Buck AK, Halter G, Schirrmeister H, Kotzerke J, Wurziger I, Glatting G, Mattfeldt T, Neumaier B, Reske SN, Hetzel M (2003) Imaging proliferation in lung tumors with PET: 18F-FLT Versus 18F-FDG. J Nucl Med 44(9):1426–1431 13. Bushong SC (1988) Magnetic resonance imaging. CV Mosby Company, St. Louis 14. Catafau M, Bullich S (2013) Molecular imaging PET and SPECT approaches for improving productivity of antipsychotic drug discovery and development. Curr Med Chem 20(3):378– 388 15. Cobelli C, Foster D, Toffolo G (2000) Tracer kinetics in biomedical research: from data to model. Kluwer Academic/Plenum Publishers, New York 16. Crum W, Hartkens T, Hill D (2004) Non-rigid image registration: theory and practice. Br J Radiol 77(Suppl 2):S140–S153 17. Daisne JF, Sibomana M, Bol A, Doumont T, Lonneux M, Grégoire V (2003) Tri-dimensional automatic segmentation of PET volumes based on measured source-to-background ratios: influence of reconstruction algorithms. Radiother Oncol 69(3):247–250 18. Defrise M, Casey ME, Michel C, Conti M (2005) Fourier rebinning of time-of-flight PET data. Phys Med Biol 50(12):2749 19. Drever L, Roa W, McEwan A, Robinson D (2007) Iterative threshold segmentation for PET target volume delineation. Med Phys 34:1253 20. Erd, YE, Mawlawi O, Larson SM, Imbriaco M, Yeung H, Finn, R, Humm JL (1997), Segmentation of lung lesion volume by adaptive positron emission tomography image thresholding. Cancer 80(S12) 2505–2509. 21. Fenster A, Downey DB (1996) 3-D ultrasound imaging: a review. IEEE Eng Med Biol Mag 15(6):41–51 22. Fessler J Erdogan H (1998) A paraboloidal surrogates algorithm for convergent penalizedlikelihood emission image reconstruction. In: Nuclear science symposium, 1998. Conference record. 1998 IEEE, vol 2, pp 1132–1135. 23. Freedman N, Schechter D, Klein M, Marciano R, Rozenman Y, Chisin R (2000) SPECT attenuation artifacts in normal and overweight persons: insights from a retrospective comparison of Rb-82 positron emission tomography and Tl-201 SPECT myocardial perfusion imaging. Clin Nucl Med 25(12):1019–1023 24. Gao F, Liu H, Jian Y, Shi P (2009) Dynamic dual-tracer PET reconstruction. In: Information processing in medical imaging (IPMI 2009). LNCS, Vol 5636. Springer, Berlin, Heidelberg, pp 38–49 25. Gao F, Liu H, Shi P (2011) Robust estimation of kinetic parameters in dynamic PET imaging. In: Medical image computing and computer-assisted intervention–MICCAI 2011, vol 6891. pp 492–499. 26. Geets X, Lee JA, Bol A, Lonneux M, Grégoire V (2007) A gradient-based method for segmenting FDG-PET images: methodology and validation. Eur J Nucl Med Mol Imaging 34(9):1427– 1438 27. Goerres GW, Kamel E, Heidelberg TNH, Schwitter MR, Burger C, von Schulthess GK (2002) PET-CT image co-registration in the thorax: influence of respiration. Eur J Nucl Med Mol Imaging 29(3):351–360 28. Gregory R, Partridge M, Flower M (2006) Performance evaluation of the Philips Gemini PET/CT system. IEEE Trans Nucl Sci 53(1):93–101 90 F. Gao and P. Shi 29. Gunn R, Gunn S, Turkheimer F, Aston J, Cunningham V (2002) Tracer kinetic modeling via basis pursuit. In: Senda M, Kimura Y, Herscovitch P (eds) Brain Imaging using PET. Academic Press, San Diego 30. Gunn RN, Gunn SR, Cunningham VJ (2001) Positron emission tomography compartmental models. J Cereb Blood Flow Metab 21(6):635–652 31. Gunn RN, Gunn SR, Turkheimer FE, Aston JA, Cunningham VJ (2002) Positron emission tomography compartmental models: A basis pursuit strategy for kinetic modeling. J Cereb Blood Flow and Metab 22(12):1425–1439 32. Hajnal JV (2001) Medical image registration. CRC Press LLC, New York 33. Han X, Xu C, Prince JL (2003) A topology preserving level set method for geometric deformable models. IEEE Trans Pattern Anal Mach Intell 25(6):755–768 34. Hara N, Okuizumi M, Koike H, Kawaguchi M, Bilim V (2005) Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is a useful modality for the precise detection and staging of early prostate cancer. The Prostate 62(2):140–147 35. Hu Z, Shi P (2010) Sensitivity analysis for biomedical models. IEEE Trans Med Imaging 29(11):1870–1881 36. Huang X, Hill NA, Ren J, Guiraudon G, Boughner D, Peters TM (2005) Dynamic 3D ultrasound and MR image registration of the beating heart. In: Medical image computing and computer-assisted intervention-MICCAI 2005. Springer, pp 171–178. 37. Huettel SA, Song AW, McCarthy G (2004) Functional magnetic resonance imaging. Sinauer Associates, Sunderland 38. Ibaraki M, Sugawara S, Nakamura K, Kinoshita F, Kinoshita T (2011) The effect of activity outside the field-of-view on image signal-to-noise ratio for 3D PET with 15O. Phys Med Biol 56(10):3061 39. Iznaga-Escobar N (1998) 188 Re-direct labeling of monoclonal antibodies for radioimmunotherapy of solid tumors: biodistribution, normal organ dosimetry, and toxicology. Nuc Med Biol 25(5):441–447 40. Jakoby B, Bercier Y, Watson C, Bendriem B, Townsend D (2009) Performance characteristics of a new LSO PET/CT scanner with extended axial field-of-view and PSF reconstruction. IEEE Trans Nucl Sci 56(3):633–639 41. Jan S, Santin G, Strul D, Staelens S, Assie K, Autret D, Avner S, Barbier R, Bardies M, Bloomfield P et al (2004) GATE: a simulation toolkit for PET and SPECT. Phys Med Biol 49(19):4543 42. Jentzen W, Freudenberg L, Eising EG, Heinze M, Brandau W, Bockisch A (2007) Segmentation of PET volumes by iterative image thresholding. J Nucl Med 48(1):108–114 43. Kelloff GJ, Sigman CC (2012) Cancer biomarkers: selecting the right drug for the right patient. Nat Rev Drug Discov 11(3):201–214 44. Kelly CJ, Brady M (2006) A model to simulate tumour oxygenation and dynamic [18F]-Fmiso PET data. Phys Med Biol 51(22):5859 45. Kitamura K, Ishikawa A, Mizuta T, Yamaya T, Yoshida E, Murayama H (2007) Detector normalization and scatter correction for the jPET-D4: a 4-layer depth-of-interaction PET scanner. Nuc Instrum Methods Phys Res Sect A Accelerators Spectrom Detect Assoc Equip 571(1–2):231–234 46. Klein S, Staring M, Murphy K, Viergever MA, Pluim JP et al (2010) Elastix: a toolbox for intensity-based medical image registration. IEEE Trans Med Imaging 29(1):196–205 47. Koeppe RA, Gilman S, Junck L, Wernette K, Frey KA et al (2008) Differentiating Alzheimer’s disease from dementia with Lewy bodies and Parkinson’s disease with (+)-[11C] dihydrotetrabenazine positron emission tomography. Alzheimer’s and dementia J Alzheimer’s Assoc 4(Suppl 1):S67 48. Krivokapich J, Smith G, Huang SC, Hoffman E, Ratib O, Phelps M, Schelbert H (1989) 13N ammonia myocardial imaging at rest and with exercise in normal volunteers. quantification of absolute myocardial perfusion with dynamic positron emission tomography. Circulation 80(5):1328–1337 Shape Analysis in Molecular Imaging 91 49. Kwakernaak H (1993) Robust control and H∞ -optimization-tutorial paper. Automatica 29(2):255–273 50. Laffon E, de Clermont H, Vernejoux JM, Jougon J, Marthan R (2011) Feasibility of assessing [18F]FDG lung metabolism with late dynamic PET imaging. Mol Imaging Biol 13(2):378– 384 51. Leff DR, Warren OJ, Enfield LC, Gibson A, Athanasiou T, Patten DK, Hebden J, Yang GZ, Darzi A (2008) Diffuse optical imaging of the healthy and diseased breast: a systematic review. Breast Cancer Res Treat 108(1):9–22 52. Lewellen T (2008) Recent developments in PET detector technology. Phys Med Biol 53:R287 53. Li Z, Li Q, Yu X, Conti P, Leahy R (2009) Lesion detection in dynamic FDG-PET using matched subspace detection. IEEE Trans Med Imaging 28(2):230–240 54. Lubberink M, Kosugi T, Schneider H, Ohba H, Bergström M (2004) Non-stationary convolution subtraction scatter correction with a dual-exponential scatter kernel for the Hamamatsu SHR-7700 animal PET scanner. Phys Med Biol 49(5):833 55. Maes F, Collignon A, Vandermeulen D, Marchal G, Suetens P (1997) Multimodality image registration by maximization of mutual information. IEEE Trans Med Imaging 16(2):187–198 56. Magota K, Kubo N, Kuge Y, Nishijima KI, Zhao S, Tamaki N (2011) Performance characterization of the Inveon preclinical small-animal PET/SPECT/CT system for multimodality imaging. Eur J Nuc Med Mol Imaging 38(4):742–752 57. Malladi R, Sethian JA, Vemuri BC (1995) Shape modeling with front propagation: a level set approach. IEEE Trans Pattern Anal Mach Intell 17(2):158–175 58. Mankoff DA (2007) A definition of molecular imaging. J Nucl Med 48(6):18N–21N 59. Maruyama A, Hasegawa S, Paul AK, Xiuli M, Yoshioka J, Maruyama K, Hori M, Nishimura T (2002) Myocardial viability assessment with gated SPECT Tc-99m tetrofosmin wall thickening: comparison with F-18 FDG-PET. Ann Nucl Med 16(1):25–32 60. Mattes D, Haynor DR, Vesselle H, Lewellen TK, Eubank W (2003) PET-CT image registration in the chest using free-form deformations. IEEE Trans Med Imaging 22(1):120–128 61. Kamasak ME, Bouman CA, Morris ED, Sauer K (2005) Direct reconstruction of kinetic parameter images from dynamic PET data. IEEE Trans Med Imaging 25:636–650 62. Mena E, Turkbey B, Mani H, Adler S, Valera VA, Bernardo M, Shah V, Pohida T, McKinney Y, Kwarteng G, Daar D, Lindenberg ML, Eclarinal P, Wade R, Linehan WM, Merino MJ, Pinto PA, Choyke PL, Kurdziel KA (2012) 11C-Acetate PET/CT in localized prostate cancer: a study with MRI and histopathologic correlation. J Nucl Med 53(4):538–545 63. Michel C, Sibomana M, Boi A, Bernard X, Lonneux M, Defrise M, Comtat C, Kinahan P, Townsend D (1998) Preserving poisson characteristics of PET data with weighted OSEM reconstruction. In: Nuclear science symposium, 1998. Conference record. 1998 IEEE, vol. 2. pp. 1323–1329 64. Moehrs S, Defrise M, Belcari N, Guerra AD, Bartoli A, Fabbri S, Zanetti G (2008) Multiray-based system matrix generation for 3D PET reconstruction. Phys Med Biol 53(23):6925 65. Naidoo-Variawa S, Lehnert W, Banati RB, Meikle SR (2011) Scatter correction for large non-human primate brain imaging using microPET. Phys Med Biol 56(7):2131 66. Nichols TE, Qi J, Asma E, Leahy RM (2002) Spatiotemporal reconstruction of list-mode PET data. IEEE Trans Med Imaging 21(4):396–404 67. Ollinger JM (1996) Model-based scatter correction for fully 3D PET. Phys Med Biol 41:153– 176 68. Osher S, Fedkiw RP (2001) Level set methods: an overview and some recent results. J Comput Phys 169(2):463–502 69. Pajares G, Manuel de la Cruz J (2004) A wavelet-based image fusion tutorial. Pattern Recogn 37(9):1855–1872 70. Panin V, Kehren F, Rothfuss H, Hu D, Michel C, Casey M (2006) PET reconstruction with system matrix derived from point source measurements. IEEE Trans Nucl Sci 53(1):152–159 71. Parker BJ, Feng D (2005) Graph-based mumford-shah segmentation of dynamic PET with application to input function estimation. IEEE Trans Nucl Sci 52(1):79–89 92 F. Gao and P. Shi 72. Parra L, Barrett HH (1998) List-mode likelihood: EM algorithm and image quality estimation demonstrated on 2-D PET. IEEE Trans Med Imaging 17(2):228–235 73. Pellizzari C, Levin DN, Chen GT, Chen CT (1993) Image registration based on anatomic surface matching. In Maciunas RJ (ed) Interactive image-guided neurosurgery. AANS, Park Ridge, pp 47–62 74. Pichler BJ, Kolb A, Nagele T, Schlemmer HP (2010) PET/MRI: paving the way for the next generation of clinical multimodality imaging applications. J Nucl Med 51(3):333–336 75. Pietrzyk U, KarlHerholz G, AndreasJacobs R (1994) Image registration: Validation for PET, SPECT, MRI and CT brain studies. J Nucl Med 35(12):2011–2018 76. Rahmim A, Rousset O, Zaidi H (2007) Strategies for motion tracking and correction in PET. PET Clin 2(2):251–266 77. Rahmim A, Tang J, Zaidi H (2009) Four-dimensional (4D) image reconstruction strategies in dynamic PET: beyond conventional independent frame reconstruction. Med Phys 36(8):3654 78. Schmitt JM (1999) Optical coherence tomography (OCT): a review. IEEE J Sel Top Quantum Electron 5(4):1205–1215 79. Shen X, Deng L (1997) Game theory approach to discrete H∞ filter design. IEEE Trans Sig Process 45:1092–1095 80. Simon D (2006) Optimal state estimation: Kalman, H∞ and nonlinear approaches. Wiley, Hoboken, New Jersey 81. Sossi V, Barney J, Harrison R, Ruth T (1995) Effect of scatter from radioactivity outside of the field of view in 3D PET. IEEE Trans Nucl Sci 42(4):1157–1161 82. Spinks TJ, Miller MP, Bailey DL, Bloomfield PM, Livieratos L, Jones T (1998) The effect of activity outside the direct field of view in a 3D-only whole-body positron tomograph. Phys Med Biol 43(4):895 83. Meikle SR, Matthews JC, Cunningham VJ, Bailey DL, Livieratos L, Jones T, Price P (1998) Parametric image reconstruction using spectral analysis of PET projection data. Phyis Med Biol 43 651–666 84. Strauss LG, Pan L, Cheng C, Haberkorn U, Dimitrakopoulou-Strauss A (2011) Shortened acquisition protocols for the quantitative assessment of the 2-tissue-compartment model using dynamic PET/CT 18F-FDG studies. J Nucl Med 52(3):379–385 85. Su KH, Wu LC, Lee JS, Liu RS, Chen JC (2009) A novel method to improve image quality for 2-D small animal PET reconstruction by correcting a monte carlo-simulated system matrix using an artificial neural network. IEEE Trans Nucl Sci 56(3):704–714 86. Sutton RS, Barto AG (1998) Reinforcement learning: an introduction. Cambridge University Press, Cambridge 87. Teras M, Tolvanen T, Johansson J, Williams J, Knuuti J (2007) Performance of the new generation of whole-body PET/CT scanners: Discovery STE and Discovery VCT. Eur J Nucl Med Mol Imaging 34(10):1683–1692 88. Tharp K, Israel O, Hausmann J, Bettman L, Martin W, Daitzchman M, Sandler M, Delbeke D (2004) Impact of 131I-SPECT/CT images obtained with an integrated system in the follow-up of patients with thyroid carcinoma. Eur J Nucl Med Mol Imaging 31(10):1435–1442 89. Townsend DW, Beyer T (2002) A combined PET/CT scanner: the path to true image fusion. Br J Radiol 75(suppl 9):S24–S30 90. Townsend DW, Carney JP, Yap JT, Hall NC (2004) PET/CT today and tomorrow. J Nucl Med 45(1 suppl):4S–14S 91. Tsoumpas C, Turkheimer F, Thielemans K (2008) A survey of approaches for direct parametric image reconstruction in emission tomography. Med Phys 35:3963 92. Valk P, Delbeke D, Bailey D, Townsend D, Maisey M (2006) Positron emission tomography: clinical practice. Springer, London 93. Villemagne VL, Rothman RB, Yokoi F, Rice KC, Matecka D, Dannals RF, Wong DF (1999) Doses of GBR12909 that suppress cocaine self-administration in non-human primates substantially occupy dopamine transporters as measured by [11C] WIN35, 428 PET scans. Synapse 32(1):44–50 Shape Analysis in Molecular Imaging 93 94. Wahl RL, Buchanan JW (2002) Principles and practice of positron emission tomography. Lippincott Williams and Wilkins, Philadelphia 95. Wahl RL, Quint LE, Cieslak RD, Aisen AM, Koeppe RA, Meyer CR et al (1993) "Anatometabolic" tumor imaging: fusion of FDG PET with CT or MRI to localize foci of increased activity. J Nucl Med 34(7):1190 96. Wang G, Qi J (2009) Generalized algorithms for direct reconstruction of parametric images from dynamic PET data. IEEE Trans Med Imaging 28(11):1717–1726 97. Wang S, Summers R (2012) Machine learning and radiology. Med Image Anal 16:933–951 98. Watson C (2000) New, faster, image-based scatter correction for 3D PET. IEEE Trans Nucl Sci 47:1587–1594 99. Wiering M, van Otterlo M (2012) Reinforcement learning: state-of-the-art, vol 12. Springer, London 100. Willmann JK, Van Bruggen N, Dinkelborg LM, Gambhir SS (2008) Molecular imaging in drug development. Nat Rev Drug Discov 7(7):591–607 101. Wong KP, Feng D, Meikle SR, Fulham MJ (2002) Segmentation of dynamic PET images using cluster analysis. IEEE Trans Nucl Sci 49(1):200–207 102. Wu H, Pal DO, Sullivan J, Tai YC (2008) A feasibility study of a prototype PET insert device to convert a general-purpose animal PET scanner to higher resolution. J Nucl Med 49(1):79–87 103. Yan, J, Planeta-Wilson B, Gallezot J, Carson R (2010) Initial evaluation of direct 4D parametric reconstruction with human PET data. In: Nuclear science symposium conference record (NSS/MIC), 2009 IEEE. pp 2503–2506 104. Yang DJ, Azhdarinia A, Kim EE (2005) Tumor specific imaging using Tc-99m and Ga-68 labeled radiopharmaceuticals. Curr Med Imaging Rev 1(1):25–34 105. Yang J, Staib LH, Duncan JS (2004) Neighbor-constrained segmentation with level set based 3-D deformable models. IEEE Trans Med Imaging 23(8):940–948 106. Zaidi H, El Naqa I (2010) PET-guided delineation of radiation therapy treatment volumes: a survey of image segmentation techniques. Eur J Nucl Med Mol Imaging 37:2165–2187 107. Zhang L, Staelens S, Holen RV, Beenhouwer JD, Verhaeghe J, Kawrakow I, Vandenberghe S (2010) Fast and memory-efficient monte carlo-based image reconstruction for whole-body PET. Med Phys 37(7):3667 http://www.springer.com/978-3-319-03812-4