|Home | About | Journals | Submit | Contact Us | Français|
Parametric imaging of the cerebral metabolic rate for glucose (CMRGlc) using [18F]-fluorodeoxyglucose positron emission tomography is considered. Traditional imaging is hindered due to low signal to noise ratios at individual voxels. We propose to minimize the total variation of the tracer uptake rates while requiring good fit of traditional Patlak equations. This minimization guarantees spatial homogeneity within brain regions and good distinction between brain regions. Brain phantom simulations demonstrate significant improvement in quality of images by the proposed method as compared to Patlak images with post-filtering using Gaussian or median filters.
We focus on positron emission tomography (PET) parametric imaging for estimating the cerebral metabolic rate of glucose (CMRGlc) using the [18F]-fluorodeoxyglucose (FDG) tracer. The ability to derive accurate parameters depends upon the quality of data, the quantification method and the numerical algorithm. In this study, we refer to the time activity curve (TAC) from a given tissue location as the output, tissue TAC or TTAC, and the TAC from the blood pool (image-derived or arterial blood-sampled) as the input, plasma TAC or PTAC. Most existing quantification methods perform well for regions of interest (ROIs), but are not good for voxel level quantification due to the high level of noise. These include graphical methods, [20, 18], linear least squares, the weighted integration method, , generalized linear least squares, [6, 5], nonlinear least squares (NLS) and weighted NLS.
All the algorithms listed above perform the quantification at each voxel location separately; they do not consider the kinetic similarities among neighboring voxels within functionally-defined regions. Thus, voxel-by-voxel variation in a functionally-homogeneous region may be large because of noise in the data. But, by incorporating the spatial constraint that parameters in a functionally-homogeneous region should be similar, in any of the above methods, the quality of the resulting parametric image for the CMRGlc may be improved. Zhou et al, , for example, improved the parametric image quality by ridge regression with constraints on the rate constants. There, the estimation of parameters uses a linear components decomposition of the kinetics in which each component represents a functional kinetic curve generated by clustering the TTACs, and the problem is solved voxel by voxel. Here we propose the use of a total variation (TV) penalty term which imposes spatial consistency between neighboring voxels.
The total variation penalty was first introduced in the context of image deblurring by Rudin et al, . TV can significantly suppress noise while recovering sharp edges because it does not penalize discontinuities. It has received much theoretical research attention and been utilized in many signal and image processing applications. While it was introduced for PET image reconstruction by Jonsson et al, , and Kisilev et al, , it has apparently not been applied for parametric PET imaging. Instead of calculating the uptake rate for each voxel by Patlak’s method, we propose to minimize the TV of the uptake rate over the entire image while also maintaining a good least squares fit for the Patlak equations at all voxels. Thus the parameters of the whole image are spatially related by the TV and solved simultaneously. The resultant parametric image is expected to have spatial homogeneity over brain regions with similar kinetics and distinct edges between brain regions that have different kinetics. This is validated by phantom simulations.
In addition to proposing the new model with the TV penalty, we also pay careful attention to the computational complexity of the algorithm by taking advantage of implementations for large scale sparse matrix computations. In contrast to approximating the Hessian matrix, as is typical for quasi-Newton methods, our algorithm explicitly and accurately calculates both gradient and Hessian terms. The Hessian is efficiently recalculated at each iteration because of its sparsity. The Quasi-Minimal Residual (QMR) method, , is used to solve the resulting large-scale linear systems. With this efficient implementation, the procedure described here is computationally feasible.
The rest of the paper is organized as follows: The new algorithm with TV penalty is introduced in section 2, with relevant computational issues detailed in the appendix. The experimental data sets are described in section 3 and results reported in section 4. Issues relevant to the proposed TV-Patlak method and computational aspects are discussed in section 5. Conclusions are presented in section 6.
The Patlak plot has been developed for systems with irreversible trapping . Most often it is applied for the analysis of FDG. The measured TTAC undergoes a transformation and is plotted against a normalized time. It is given by the expression
where CT(t) is the measured TTAC, (in counts/min/g) and CP is the PTAC (in counts/min/ml), i.e. the FDG concentration in plasma. For systems with irreversible compartments this plot yields a straight line after sufficient equilibration time. For the FDG tracer, the slope K represents the uptake rate which, together with lumped constant (LC) and glucose concentration in plasma (Cpg) allows easy calculation of the CMRGlc=KCpg/LC, (in mg/min/100g). The intercept V is given by V0 + vB where V0 is the distribution volume of the reversible compartment and vB is the fractional blood volume.
The linear relationship (1) can be rewritten as
and, assuming m dynamic frames over the period of equilibration, its discretized version is
In matrix format,
where A is a m-by-2 matrix,
and vector b = (CT(t1), CT(t2), ···, CT(tm))T. If we were to solve (4) for each voxel independently we would obtain a parametric image lacking spatial homogeneity and with low signal-to-noise ratio (SNR). Image denoising techniques could then be applied as a separate task to improve the image quality. Instead, obtaining all voxel parameters as a result of a global optimization algorithm with a TV penalty for the entire image, the necessity for postprocessing should be eliminated.
Limiting the discussion here to 2D images (although our application of TV is 3D), we select the active voxels to be quantified by the application of a brain mask, yielding a total of N voxels. Equation (4) holds with common matrix A dependent on CP(t) for each voxel i, but with K, V and b replaced by K(i), V(i), and bi, respectively, where bi is obtained from the TTAC for voxel i. Collecting the unknowns of these N voxels in vectors of uptake rates and intercepts
and requiring (4) in the least squares sense over all voxels, while maintaining minimal TV of the uptake rate for the selected image voxels, yields the global minimization problem
Here the total variation norm is given by with , and xir and xib are the values associated with voxels to the right and below voxel i. Theoretically the TV norm is ||x||TV,0, which is a seminorm on a space of bounded variation, . The small constant β is used to avoid the numerical difficulty due to the lack of differentiability at the origin of i for β = 0. The diagonal weight matrix for voxel i is given by
where Δtj is the scan duration of the frame at time tj, λ is the tracer’s decay constant and bij is the value of the ith TTAC at frame j. This weighting is consistent with using a simulation with variance
Further details on the calculation of gradient vector Φ(z) and Hessian matrix 2Φ(z) are provided in the appendix. Some other aspects of the algorithm, including discussion about constants α and β, here chosen to be 0.2 and 10−8, respectively, as well as other approaches for the solution of the TV problem are discussed in section 5.
To validate the proposed parametric imaging method we performed experiments with simulated data. The MRI-based high-resolution Zubal head phantom is used to define the brain structures, . Each voxel in the slice of 256 ×256 voxels is of size 1.5 ×1.5mm2. There are 128 slices for the entire head and 62 defined anatomical, neurological, and taxonomical structures. The 11 regions for slice 64, representing a total 13708 voxels, are given in table 1. For the purposes of the simulations, well-accepted values of the kinetic rate parameters of these structures for the two-tissue compartmental model of the FDG tracer , are assigned. Notice that in order to better model the true biochemical process we assume k4 > 0. Because k4 is, however, relatively small, Patlak’s method, for which it is assumed k4 = 0, is still a suitable graphical method for quantifying the uptake rate.
The noise-free input function, with values given in kBq/ml, is given by the formulation introduced in ,
Although any reasonable input function, including clinical plasma samples, could be used for the simulation this formulation was validated as providing a good approximation to the plasma samples of a healthy subject. Given the input function, the exact phantom can then be generated using the rate constants for the structures detailed in table 1, . The output time frames were generated assuming time frames with durations Δtj, j = 1, ···, m, m = 22, given in minutes, 0.2, 8 × 0.0333, 2 × 0.1667, 0.2, 0.5, 2 × 1, 2 × 1.5, 3.5, 2 × 5, 10 and 30.
In our experiments, in order to control the computational overhead of the reconstructions, we reduced the image size of the phantom from 256 ×256 to 128 × 128. To do this we needed to relabel the structure assigned to a given voxel. This was done by averaging the quantity K1k3/(k2 + k3) over the 2-by-2 neighboring voxels at the finer resolution. Then the voxel at the coarser level was labeled as belonging to the structure for which this average is closest to the structure value. The kinetic values were then assigned, the TTAC output values calculated, and then projected to the sinogram space yielding projected noise-free sinogram data P. For simplification, the instrumental and physical effects, including attenuation, Compton scattering, decay and random coincidences, were not simulated. Poisson noise was then added to the projected data using S = poissrnd(P) where the Matlab®, function poissrnd uses vector P as the means of Poisson densities to generate the noisy sinogram data S. Based on this noisy sinogram data, concentration images are reconstructed using the Expectation-Maximization (EM) algorithm, .
Several data sets were generated: To investigate errors introduced by violation of the irreversibility assumption, k4 = 0, we tested both k4 = 0 and k4 > 0. Similarly, we tested both with and without noise in the sinograms to investigate the effects of the proposed method due to noise in the sinograms. Gaussian noise with noise levels 0%, 5%, 10%, 15% and 20% was added to the input function, i.e.
where ηj is selected from a standard normal distribution (G(0, 1)), and CV = 0, 0.05, 0.10, 0.15 and 0.20. 100 random realizations are tested in each case.
We summarize the test data sets in table 2, using a character triple to classify each test. The first character of the triple indicates whether k4 = 0 or k4 > 0, 0 or +, respectively. The second character indicates the noise level on the input function, 0 or +, for noise-free, or with added noise, respectively. To indicate the noise level + is replaced by 1, 2, 3 and 4 to indicate noise levels 5%, 10%, 15% and 20% as necessary. The third character indicates whether noise is added to the sinogram, again 0 or +, respectively. Therefore, for example, +3+ represents the test case for k4 > 0, 15% noise in u(t) and Poisson noise added to the sinogram.
To evaluate the simulations quantitatively we define the relative error of the kth realization for voxel i:
where is the estimated value of the true value of the uptake rate at the kth realization for voxel i. Over 100 random simulations, and N voxels, we calculate the bias, i.e. the mean of the relative errors, and the deviation from the mean:
Absolute relative errors |rik| and associated mean , and deviation, , are also calculated. Note the deviations are l1 measurements which do not overweight outlier and large error samples, as is the case for the l2-based measurements such as the root mean squared error.
In the images shown in the figures we illustrate the calculated uptake rates K of the FDG. Images for the CMRGlc can be obtained by directly scaling K. In figure 1 we compare the result of using Patlak and TV-Patlak for estimating the uptake rates with respect to no noise, 20% noise in the input function, Poisson noise in the sinogram, and finally with respect to the case in which the irreversibility assumption is violated but without noise in the sinogram or input data. In each case the histogram of the relative errors is given on the left, the Patlak image in the middle and the TV-Patlak on the right. The different scales in the histograms are due to the total number of results illustrated. When there is no noise (triples 000 and +00) the histogram illustrates results over all voxels but only one simulation, while for the noisy simulations the results are for all voxels over all 100 realizations of the noise. The TV-Patlak images are more homogeneous in all cases and the relative errors are smaller. The figures clearly show the improvements of employing the TV-Patlak method as compared to using Patlak independently for each voxel. This is confirmed in figure 2 in which images with noise in the sinogram, positive k4 and different noise levels in the input function are shown.
Quantitative measurements, confirming the illustrations, are presented in table 3. There we also present the results for conventional Patlak’s method with post-smoothing by two standard filters:
Consistent with the observation in [21, 16], we find that violation of the Patlak assumption, k4 = 0, introduces about 10% bias; ≈ 0 when k4 = 0 but ≈ 10% for k4 > 0. The rows (std) and “# 10% (#15%)” provide complementary supporting information, indicating that the TV is minimized by TV-Patlak; as compared to Patlak, Patlak-GF and Patlak-MF the number of voxels with larger error is reduced. In particular, we emphasize that TV-Patlak provides a better noise removal mechanism than popular post-filtering approaches.
In figures 3 and and44 we illustrate the uptake rates and relative error in the uptake rates, respectively, calculated by Patlak, TV-Patlak, Patlak-GF and Patlak-MF for one simulated data case +4+, i.e. k4 > 0, 20% noise in the input function and Poisson noise in the sinograms. The uptake rate image generated by Patlak-MF is visually smoother than that by TV-Patlak, but the equivalent histograms show that the relative error is higher for Patlak-MF than for TV-Patlak; the Patlak-MF image is over-smoothed.
In figure 5 we illustrate the relative error of noisy case +4+ for gray matter regions in the phantom, including frontal lobes, occipital lobes, insula cortex, temporal lobes and globus pallidus, which are of interest for Alzheimer’s disease research. The error bars show that all estimation methods for gray matter regions have negative bias and there are fewer cases with high error by TV-Patlak.
Finally, we note that the computational cost for the TV-Patlak Algorithm 1 is about 6.5 seconds while for Patlak, Patlak-GF and Patlak-MF they are 0.48, 0.55 and 0.58 seconds respectively for a 2D image on a PC with 1GHz CPU and Matlab code.
In this section we discuss issues relevant to the proposed TV-Patlak method.
A qualitative improvement in imaging of PET uptake can be achieved by using a global model, with the total variation as a penalty term, to obtain the voxel uptake rates. The resultant uptake images have spatial homogeneity over brain regions with similar kinetics and distinct edges between brain regions that have different kinetics. It is statistically validated that the TV-Patlak significantly reduces the relative errors of the calculated parameters as compared with those generated by Patlak’s graphical method, and post-smoothing by Gaussian and median filters.
This work was supported by grants from the state of Arizona, the NIH, R01 MH057899 and P30 AG19610 and the NSF DMS 0652833 and DMS 0513214. The authors thank Chi-Chuan Chen for providing the transition matrix, generated by the Monte Carlo method, which is used in the PET image reconstructions and Drs Bouman and Musfata for offering their MAP reconstruction code, which provides comparisons between reconstruction methods. We also acknowledge Dr. Laszlo Balkay and Dr. Jeffrey Fessler for providing PET simulator and PET reconstruction software, respectively.
Hongbin Guo received his Ph.D. degree in computational mathematics from the Department of Mathematics, Fudan University, in 2000. He is currently an Assistant Research Professor with the Department of Mathematics and Statistics at Arizona State University. His main research interest is computational mathematics with a focus on the numerical linear algebra and its applications in medical imaging, particularly for the quantification of brain function with positron emission tomography and magnetic resonance imaging, co-registration, functional clustering, image restoration and classification.
Rosemary A. Renaut received the Ph.D. degree in applied mathematics from the University of Cambridge, U.K., in 1985. Since 1987, she has been with the Department of Mathematics, Arizona State University, where she is now a Full Professor. Dr. Renaut is a Fellow of the Institute for Mathematics and its Applications, and a Chartered Mathematician. Her research interests are broad and include the design and evaluation of computational methods for the solution of partial differential equations, with specific emphasis on high-order and spectral methods, as well as the development of novel algorithms for solving inverse problems, specifically as applied to medical image reconstruction and restoration.
Kewei Chen got his master degree from Beijing Normal University in 1986 and his Ph.D. Degree from UCLA in 1993. He is currently the director and a senior biomathematician of the Computational Image Analysis Program, Banner Alzheimer’s Institute. He has been using and developing neuroimaging analytic techniques in PET, magnetic resonance imaging (MRI) and functional MRI. His main methodological research interests include the tracer kinetic modeling in PET, measuring local and global volume changes using various MRI techniques, brain functional connectivity, and multi-modal data integration. One of his primary research interests is the use of neuroimaging techniques in the study of Alzheimer’s disease.
Eric M Reiman is Executive Director of the Banner Alzheimer’s Institute, Clinical Director of the Neurogenomics Division at the Translational Genomics Research Institute (TGen), Professor and Associate Head of Psychiatry at the University of Arizona, and Director of the Arizona Alzheimer’s Consortium. He received his undergraduate and medical degrees and most of psychiatry residency training at Duke University. He completed his residency training, became an Assistant Professor of Psychiatry and developed a leadership role in positron emission tomography research at Washington University in St. Louis, before moving to Arizona. His research interests include brain imaging, genomics, and their use in the unusually early detection and tracking of Alzheimer’s disease and the rigorous and rapid evaluation of promising Alzheimer’s disease-slowing, risk-reducing and prevention therapies.
To simplify notation, we assume that the weighting matrix in (5) has been absorbed into A and bi. To simplify the expressions in the objective function (6) we introduce the vector function f(z) = (f1(z), f2(z), ···, f2N(z))T, where fi(z) = k(x)2 for i = 1 … N and fi(z) = ||A[xi−N; yi−N] − bi−N||2 for i = N + 1 … 2N.
To find the gradient of Φ(z), we first derive the Jacobian of f(z).
To find the Hessian matrix for Φ(z), we first derive the Hessian matrix for each function fi(z).